4.3.2　用戶子程序UMAT和應用程序（Utilities）

1）用戶子程序（User Subroutines）和應用程序（Utilities）

雖然ABAQUS為用戶提供了大量的單元庫和材料模型（如金屬、橡膠、塑料、混凝土、巖土等），使用戶能夠利用這些模型處理絕大多數(shù)的問題，但上述道路工程（包含部分巖土工程）中常見材料模型，并沒有包含在ABAQUS軟件中。為了彌補這一不足，ABAQUS提供了大量的用戶子程序（User Subroutines）和應用程序（Utilities），用戶可以自行定義符合特定問題的模型。

（1）用戶材料子程序UMAT

用戶材料子程序UMAT（User－Defined Material Mechanical Behavior）是ABAQUS提供給用戶進行材料本構(gòu)模型二次開發(fā)的一個用戶子程序接口，可以定義用戶所需要的各類材料本構(gòu)模型，這大大增強了ABAQUS的應用面和靈活性。

用戶子程序UMAT具有如下特點：

·用來定義材料的本構(gòu)關系；

·當材料的定義包含用戶自定義材料模型時，每一個計算單元的材料積分點都可以調(diào)用UMAT；

·可以用于力學行為分析的任何分析過程；

·可以使用狀態(tài)變量；

·對于力學本構(gòu)關系，必須在UMAT中提供材料本構(gòu)模型的雅可比矩陣

·可以和用戶子程序USDFLD聯(lián)合使用，通過USDFLD重新定義任何常變量值并傳遞到UMAT。

UMAT子程序的核心內(nèi)容就是給出定義材料本構(gòu)模型的雅可比矩陣（Jacobian矩陣，即應力增量對應變增量的變化率aΔσ

aΔε，Δσ是應力增量，Δε是應變增量），并更新應力提供給ABAQUS主程序。例如，已知第n步的結(jié)果σ_n及ε_n等，然后ABAQUS主程序給出一個應變增量dε_n＋1，UMAT根據(jù)提供的雅可比矩陣DDSDDE計算出新的應力σ_n＋1。

UMAT子程序采用Fortran語言編寫，從主程序獲取數(shù)據(jù)，計算單元的材料積分點的雅可比矩陣，并更新應力張量和狀態(tài)變量，UMAT的接口格式如下：

SUBROUTINE UMAT（STRESS，STATEV，DDSDDE，SSE，SPD，SCD，RPL，DDSDDT，

1DRPLDE，DRPLDT，STRAN，DSTRAN，TIME，DTIME，TEMP，DTEMP，PREDEF，DPRED，

2CMNAME，NDI，NSHR，NTENS，NSTATV，PROPS，NPROPS，COORDS，DROT，

3PNEWDT，CELENT，DFGRD0，DFGRD1，NOEL，NPT，LAYER，KSPT，KSTEP，KINC）

C

INCLUDE′ABA＿PARAM.INC′

C

CHARACTER＊80CMNAME

C

DIMENSIONSTRESS（NTENS），STATEV（NSTATV），DDSDDE（NTENS，NTENS），

1DDSDDT（NTENS），DRPLDE（NTENS），STRAN（NTENS），DSTRAN（NTENS），

2TIME（2），PREDEF（1），DPRED（1），PROPS（NPROPS），COORDS（3），DROT（3，3），

3DFGRD0（3，3），DFGRD1（3，3）

User coding to define DDSDDE，STRESS，STATEV，SSE，SSD，SCD

and，if necessary，RPL，DDSDDT，DRPLDE，DRPLDT，PNEWDT

RETURN

END

UMAT子程序主要變量說明：

·DDSDDE（NTENS，NTENS）

是一個NTENS×NTENS維的方陣，稱為本構(gòu)關系的雅可比矩陣，DDSDDE（I，J）表示增量步結(jié)束時第J個應變分量的改變引起的第I個應力分量的變化。通常雅可比矩陣是一個對稱矩陣。

·STRESS（NTENS）

應力張量數(shù)組，對應NDI個直接分量和NSHR個剪切分量。增量步開始時，該數(shù)組從ABAQUS主程序獲取數(shù)據(jù)σ_n并作為已知量傳入UMAT；增量步結(jié)束時，在子程序中更新應力張量為σ_n＋1。如果定義了初始應力，那么分析開始時該應力張量的數(shù)值即為初始應力。對于包含剛體轉(zhuǎn)動的有限應變問題，一個增量步調(diào)用UMAT之前就已經(jīng)對應力張量進行了剛體轉(zhuǎn)動，因此在UMAT中只需處理應力張量的共旋部分。UMAT中應力張量的度量為柯西（真實）應力。

·STATEV（NTENS）

狀態(tài)變量數(shù)組，用于存儲狀態(tài)變量的數(shù)組，數(shù)值在增量步開始時從主程序傳遞到UMAT，也可以在USDFLD或UEXPAN中先更新數(shù)據(jù)，然后增量步開始時將更新的數(shù)據(jù)傳遞到UMAT子程序中。在增量步結(jié)束時，必須更新STATEV中的數(shù)據(jù)返回給主程序。狀態(tài)變量數(shù)組的維數(shù)通過ABAQUS輸入文件中的關鍵字“＊DEPVAR”定義，關鍵字下面數(shù)據(jù)行的數(shù)值即為狀態(tài)變量數(shù)組的維數(shù)。狀態(tài)變量數(shù)組是用來保存用戶自己定義的一些變量，如累計塑性應變、黏彈性應變等。

·STRAN（NTENS）

總應變數(shù)組，存儲增量步開始時的總應變ε_n，由ABAQUS主程序自動更新。

·DSTRAN（NTENS）

總應變增量dε_n＋1。

·PROPS（NPROPS）

材料常數(shù)矩陣，用來保存用戶自定義材料參數(shù)，即用戶自定義材料本構(gòu)關系的模型參數(shù)。材料常數(shù)的個數(shù)NPROPS等于關鍵字“＊USER MATERIAL”中的“CONSTANTS”參數(shù)設定的值，矩陣中元素的數(shù)值對應關鍵字“＊USER MATERIAL”下面的數(shù)據(jù)行。

·SSE，SPD，SCD

分別定義每一增量步的彈性應變能、塑性耗散和蠕變耗散。它們對計算結(jié)果沒有影響，僅僅作為能量輸出。

·DTIME

時間增量dt。

·NDI

法向應力、應變個數(shù)，三維問題、軸對稱問題是3（11，22，33），平面問題是2（11，22）。

·NSHR

剪切應力、應變個數(shù)，三維問題是3（12，13，23），軸對稱問題是1（12）。

·NTENS

應力和應變分量數(shù)組的大?。∟TENS＝NDI＋NSHR）。

·DROT

對Finite strain問題，應變應該排除旋轉(zhuǎn)部分，該矩陣提供了旋轉(zhuǎn)矩陣。

·PNEWDT

可用來控制時間步的變化。如果設置為小于1的數(shù)，則程序放棄當前計算，并用新的時間增量DTIME×PNEWDT作為新的時間增量計算；對時間相關的材料如聚合物等有用；如果設為大于1的數(shù)，則下一個增量步加大DTIME為DTIME×PNEWDT。

·KSTEP，KINC

ABAQUS傳到UMAT的當前分析步和增量步數(shù)值。

·TIME（1），TIME（2）

當前分析步時間（STEP TIME）和增量步時間（INCREMENT TIME）值。

（2）應用程序（Utilities）

ABAQUS提供了一些實用的應用程序（Utility Routines），供用戶在開發(fā)用戶子程序UMAT時調(diào)用，這些程序包括應力不變量的計算，主應力的計算等。

·SINV：用于計算應力不變量

CALLSINV（STRESS，SINV1，SINV2，NDI，NSHR）

STRESS：應力張量；

NDI：法向應力分量的數(shù)量；

NSHR：剪切應力分量的數(shù)量；

SINV1：第一應力不變量，SINV1＝traceσ，其中σ為應力張量；

SINV2：第二應力不變量，SINV2＝，其中S為偏應力張量，S＝σ－traceσI；

·SPRINC：用于計算主應力值

CALLSPRINC（S，PS，LSTR，NDI，NSHR）

S：應力或應變張量；

PS（I），I＝1，2，3：三個主應力值；

LSTR：標識，LSTR＝1表示S為應力張量，LSTR＝2表示S為應變張量。

·SPRIND：用于計算主應力值和主應力方向

CALLSPRIND（S，PS，AN，LSTR，NDI，NSHR）

S：應力或應變張量；

AN（K1，I），I＝1，2，3：PS（K1）法向應力的方向余弦。

2）ABAQUS主程序與UMAT子程序協(xié)同工作過程

ABAQUS主程序與UMAT子程序之間是一個動態(tài)交互傳遞數(shù)據(jù)、協(xié)同工作的過程。UMAT子程序作為ABAQUS主程序的一個接口，在單元的積分點上調(diào)用，增量步開始時，主程序通過UMAT的接口進入UMAT子程序，單元當前積分點必要變量的初始值將隨之傳遞給UMAT子程序的相應變量，UMAT子程序計算單元材料積分點的雅可比矩陣，并更新應力張量和狀態(tài)變量，最后將這些變量的更新值通過接口返回主程序。

ABAQUS主程序與UMAT子程序的交互計算過程如下：從第t_n時刻開始，ABAQUS在Δt時間增量內(nèi)產(chǎn)生一個由外載荷產(chǎn)生的應變增量Δε，UMAT子程序通過給定的本構(gòu)方程為主程序提供新的柯西應力張量σ（t_n＋Δt），如果計算的應力－應變結(jié)果收斂，那么ABAQUS主程序繼續(xù)計算第t_n＋1步，并根據(jù)上一步的收斂情況來選取下一步增量步長。雅可比矩陣即DDSDDE的精度影響程序的收斂速度，但是并不影響計算結(jié)果的準確性。

以下節(jié)研究的Burgers蠕變模型為例，某一材料單元積分點上ABAQUS主程序與UMAT子程序協(xié)同工作過程，如圖4.12所示。

圖4.12　ABAQUS主程序與UMAT子程序協(xié)同工作過程

第一步：平衡時刻t_n，ABAQUS主程序提供給UMAT子程序總應變ε（t_n）和應力σ（t_n），同時ABAQUS主程序自動產(chǎn)生一個時間增量Δt。

第二步：ABAQUS主程序產(chǎn)生一個總應變增量Δε（t_n），調(diào)用UMAT子程序計算Δε（t_n）對應的蠕應變增量Δε_c（t_n）及應力增量Δσ（t_n），此過程是一個非線性迭代求解的過程。迭代收斂后更新應力張量σ（t_n＋1）＝σ（t_n）＋Δσ（t_n）。

第三步：UMAT子程序?qū)⒏潞蟮膽Ζ遥╰_n＋1）返回給ABAQUS主程序，同時ABAQUS主程序自動更新總應變ε（t_n＋1）＝ε（t_n）＋Δε（t_n）。

第四步：ABAQUS主程序進行最大迭代次數(shù)的檢查（控制最大迭代次數(shù)是出于計算成本的考慮），如果迭代次數(shù)超過了限定的最大迭代次數(shù)（ABAQUS默認的最大迭代次數(shù)為16次），那么ABAQUS主程序就會自動減小時間增量，然后返回到第一步重新進行迭代平衡；然后ABAQUS主程序?qū)ⅵ遥╰_n＋1）代入系統(tǒng)平衡方程進行平衡判斷，若滿足系統(tǒng)平衡標準，結(jié)束本增量步的迭代，進入下一增量步。若不滿足系統(tǒng)平衡，此時ABAQUS主程序?qū)⒎艞壉敬蔚膽?、應變更新，回復到增量步初始時刻的值，并繼續(xù)進行本增量步的下一次平衡迭代，直到滿足系統(tǒng)平衡標準為止。