[交通運(yùn)輸]21ABAQUS用戶材料子程序

1、.21 ABAQUS用戶材料子程序(UMAT)雖然ABAQUS為用戶提供了大量的單元庫和求解模型，使用戶能夠利用這些模型處理絕大多數(shù)的問題；但是實(shí)際問題畢竟非常復(fù)雜，ABAQUS不可能直接求解所有可能出現(xiàn)的問題。所以ABAQUS提供了大量的用戶自定義子程序（User Subroutine），允許用戶在找不到合適模型的情況下自行定義符合自己問題的模型。這些用戶子程序涵蓋了建模、載荷到單元的幾乎各個部分。用戶子程序具有以下的功能和特點(diǎn)：（1）如果ABAQUS的一些固有選項(xiàng)模型功能有限，用戶子程序可以提高ABAQUS中這些選項(xiàng)的功能；（2）通常用戶子程序是用FORTRAN語言的代碼寫成；（3）它可以

2、以幾種不同的方式包含在模型中；（4）由于它們沒有存儲在restart文件中，如果需要的話，可以在重新開始運(yùn)行時修改它；（5）在某些情況下它可以利用ABAQUS允許的已有程序。要在模型中包含用戶子程序，可以利用ABAQUS執(zhí)行程序，在執(zhí)行程序中應(yīng)用user選項(xiàng)指明包含這些子程序的FORTRAN源程序或者目標(biāo)程序的名字。提示：ABAQUS的輸入文件除了可以通過ABAQUS/CAE的作業(yè)模塊中提交運(yùn)行外，還可以在ABAQUS Command窗口中輸入ABAQUS執(zhí)行程序直接運(yùn)行：ABAQUS job=輸入文件名 user=用戶子程序的Fortran文件名ABAQUS/Standard和ABAQUS/

3、Explicit都支持用戶子程序功能，但是他們所支持的用戶子程序種類不盡相同，讀者在需要使用時請注意查詢手冊。在接下來的兩章里，我們將討論兩種常用的用戶子程序用戶材料子程序和用戶單元子程序。本章將通過在ABAQUS/Standard中創(chuàng)建Johnson-Cook的材料模型，介紹編寫ABAQUS/Standard的用戶材料子程序UMAT。在ABAQUS/Explicit中編寫用戶材料子程序VUMAT與之相似，但是由于隱式和顯式兩種方法本身的差異，它們之間也有一些不同，請讀者在具體使用前仔細(xì)查閱ABAQUS手冊中的相關(guān)內(nèi)容。21.1 引言用戶材料子程序是ABAQUS提供給用戶自定義材料屬性的For

4、tran程序接口，它使用戶能使用ABAQUS材料庫中沒有的材料模型。ABAQUS中自有的Johnson-Cook模型只能應(yīng)用于顯式ABAQUS/Explicit程序中，而我們希望能在隱式ABAQUS/Standard程序中更精確地實(shí)現(xiàn)本構(gòu)積分，而且應(yīng)用Johnson-Cook模型的修正形式。這就需要通過ABAQUS/Standard的用戶材料子程序UMAT編程實(shí)現(xiàn)。在UMAT編程中使用了率相關(guān)塑性理論以及完全隱式的應(yīng)力更新算法。21.2 模型的數(shù)學(xué)描述21.2.1 Johnson-Cook強(qiáng)化模型簡介Johnson-Cook模型用來模擬在沖擊載荷作用下的變形。Johnson-Cook強(qiáng)化模型（

5、JC）表示為三項(xiàng)的乘積，分別反映了應(yīng)變硬化、應(yīng)變率硬化和溫度軟化。這里使用JC模型的修正形式：(21-1)模型中包含五個參數(shù)，需要通過實(shí)驗(yàn)來確定。使參考應(yīng)變率，這樣公式中的即為材料的靜態(tài)屈服應(yīng)力。公式中的為無量綱化的溫度其中為室溫，為材料的熔點(diǎn)。Johnson-Cook模型在溫度從室溫到材料熔點(diǎn)溫度的范圍內(nèi)都是有效的。高應(yīng)變率的變形經(jīng)常伴有溫升現(xiàn)象，這是因?yàn)椴牧献冃芜^程中塑性功轉(zhuǎn)化為熱量。對于大多數(shù)金屬，90-100的塑性變形將耗散為熱量。所以JC模型中溫度的變化可以用如下的公式計算：(21-2)其中，為溫度的增量，為塑性耗散比，表示塑性功轉(zhuǎn)化為熱量的比例，為材料的比熱，為材料的密度。公式（2

6、1-2）考慮的是一個絕熱過程，即認(rèn)為溫度的升高完全起因于塑性耗散。21.2.2 率相關(guān)塑性的基本公式Johnson-Cook本構(gòu)模型考慮率相關(guān)塑性，塑性變形是關(guān)聯(lián)的，即塑性流動沿著屈服面的法線方向，并采用Mises屈服面。將應(yīng)變的增量分解為彈性部分和塑性部分：(21-3)將上式兩端同時對時間的增量微分得到率形式：(21-4)在率相關(guān)塑性中，材料的塑性反應(yīng)取決于加載率，以率的形式給出材料的彈性反應(yīng)：(21-5)為了發(fā)生塑性變形，率相關(guān)塑性必須滿足或者超過屈服條件，塑性應(yīng)變率為：，(21-6)上式即為流動法則，其中為塑性流動勢能，為塑性率參數(shù)，為運(yùn)動硬化時的背應(yīng)力。對于各項(xiàng)同性硬化，不存在背應(yīng)力，

7、因此有，此時有：(21-7) (21-8)與率無關(guān)塑性不同的是，率相關(guān)塑性中等效塑性應(yīng)變率不能通過一致性條件獲得，而是直接通過經(jīng)驗(yàn)定律給出，成為過應(yīng)力模型： (21-9)式中是過應(yīng)力，為粘性。在過應(yīng)力模型中，等效塑性應(yīng)變率取決于超過了多少屈服應(yīng)力。上面一維的率相關(guān)塑性公式可以很方便地推廣到三維情況。對于小應(yīng)變的情形，應(yīng)力度量之間無需區(qū)分，這里采用Cauchy應(yīng)力，塑性率參數(shù)由應(yīng)力和內(nèi)變量的經(jīng)驗(yàn)函數(shù)給出。對照一維情況，三維情況下分解應(yīng)變率為彈性和塑性部分：(21-10)應(yīng)力率和彈性應(yīng)變率之間的關(guān)系為：(21-11)塑性流動法則和內(nèi)變量的演化方程為：，(21-12)塑性率參數(shù)為：(21-13)對于

8、J2流動理論，Perzyna（1971）中提出了典型的過應(yīng)力模型為： (21-14)式中為Macualay括號，如果，則；如果，則。為Mises等效應(yīng)力，為等效應(yīng)變，為率敏感系數(shù)。對于Johnson-Cook模型，可以得到等效塑性應(yīng)變率的表達(dá)式為：(21-15)其中為靜態(tài)屈服應(yīng)力。(21-16)21.2.3 完全隱式的應(yīng)力更新算法對率形式的本構(gòu)方程進(jìn)行積分的算法稱為應(yīng)力更新算法。在完全隱式的算法中，在步驟結(jié)束時計算塑性應(yīng)變和內(nèi)變量的增量，同時強(qiáng)化屈服條件。積分算法寫為： (21-17a) (21-17b)(21-17c)(21-17d)(21-17e)在時刻給出一組和應(yīng)變增量，公式（21-17

9、）是一組關(guān)于求解的非線性代數(shù)方程。將公式（21-17b）代入（21-17d）得到：(21-18)式中是彈性預(yù)測的試應(yīng)力，而數(shù)值是塑性修正量，它沿著結(jié)束點(diǎn)塑性流動的方向。對于J2流動理論，塑性流動的方向?yàn)椋?21-19)它是屈服面的法向，即。在偏應(yīng)力空間，von Mises屈服面為環(huán)狀，所以屈服面的法向通過圓心，如圖21-1所示。圖21-1 徑向返回算法從圖21-1可以看出，在彈性預(yù)測階段，塑性應(yīng)變和內(nèi)變量保持固定；而在塑性修正階段，總體應(yīng)變保持不變。在迭代開始時，程序?qū)?yīng)力和應(yīng)變設(shè)初始值為：，(21-20)應(yīng)力在第次迭代時為：(21-21)定義屈服面的單位法向矢量為：，(21-22)且在整個算

10、法的塑性修正狀態(tài)過程中始終保持不變，因此塑性應(yīng)變的更新是的線性函數(shù)。在次迭代時將檢查屈服條件：(21-23)若收斂，則迭代完畢，增量步結(jié)束。否則將計算塑性參數(shù)的增量： (21-24)并對塑性應(yīng)變和內(nèi)變量進(jìn)一步更新：， (21-25a)(21-25b) (21-25c)(21-25d)(21-25e)然后將更新的變量返回屈服條件進(jìn)行檢查，整個過程將重復(fù)直至收斂為止。這就是增量步中應(yīng)力更新的過程。21.3 ABAQUS用戶材料子程序用戶材料子程序（User-defined Material Mechanical Behavior，簡稱UMAT）通過與ABAQUS主求解程序的接口實(shí)現(xiàn)與ABAQUS的

11、數(shù)據(jù)交流。在輸入文件中，使用關(guān)鍵字“*USER MATERIAL”表示定義用戶材料屬性。21.3.1 子程序概況與接口UMAT子程序具有強(qiáng)大的功能，使用UMAT子程序：(1) 可以定義材料的本構(gòu)關(guān)系，使用ABAQUS材料庫中沒有包含的材料進(jìn)行計算，擴(kuò)充程序功能；(2) 幾乎可以用于力學(xué)行為分析的任何分析過程，幾乎可以把用戶材料屬性賦予ABAQUS中的任何單元；(3) 必須在UMAT中提供材料本構(gòu)模型的雅可比（Jacobian）矩陣，即應(yīng)力增量對應(yīng)變增量的變化率；(4) 可以和用戶子程序“USDFLD”聯(lián)合使用，通過“USDFLD”重新定義單元每一物質(zhì)點(diǎn)上傳遞到UMAT中場變量的數(shù)值。由于主程序

12、與UMAT之間存在數(shù)據(jù)傳遞，甚至共用一些變量，因此必須遵守有關(guān)UMAT的書寫格式，UMAT中常用的變量在文件開頭予以定義，通常格式為：SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KS

13、PT,KSTEP,KINC)C INCLUDE 'ABA_PARAM.INC'C CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) user coding to define D

14、DSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT RETURN ENDUMAT中的應(yīng)力矩陣、應(yīng)變矩陣以及矩陣，等，都是直接分量存儲在前，剪切分量存儲在后。直接分量有個，剪切分量有個。各分量之間的順序根據(jù)單元自由度的不同有一些差異，所以編寫UMAT時要考慮到所使用單元的類別。下面對UMAT中用到的一些變量進(jìn)行說明：是一個維的方陣，稱作雅可比矩陣，即，是應(yīng)力的增量，是應(yīng)變的增量，表示增量步結(jié)束時第個應(yīng)變分量的改變引起的第個應(yīng)力分量的變化。通常雅可比是一個對稱矩陣，除

15、非在“*USER MATERIAL”語句中加入了“UNSYMM”參數(shù)。應(yīng)力張量矩陣，對應(yīng)個直接分量和個剪切分量。在增量步的開始，應(yīng)力張量矩陣中的數(shù)值通過UMAT和主程序之間的接口傳遞到UMAT中；在增量步的結(jié)束，UMAT將對應(yīng)力張量矩陣更新。對于包含剛體轉(zhuǎn)動的有限應(yīng)變問題，一個增量步調(diào)用UMAT之前就已經(jīng)對應(yīng)力張量進(jìn)行了剛體轉(zhuǎn)動，因此在UMAT中只需處理應(yīng)力張量的共旋部分。UMAT中應(yīng)力張量的度量為柯西（真實(shí)）應(yīng)力。用于存儲狀態(tài)變量的矩陣，在增量步開始時將數(shù)值傳遞到UMAT中。也可在子程序USDFLD或UEXPAN中先更新數(shù)據(jù)，然后在增量步開始時將更新后的數(shù)據(jù)傳遞到UMAT中。在增量步結(jié)束時必

16、須更新狀態(tài)變量矩陣中的數(shù)據(jù)。和應(yīng)力張量矩陣不同的是：對于有限應(yīng)變問題，除了材料本構(gòu)行為引起的數(shù)據(jù)更新以外，狀態(tài)變量矩陣中的任何矢量或者張量都必須通過旋轉(zhuǎn)來考慮材料的剛體運(yùn)動。狀態(tài)變量矩陣的維數(shù)，等于關(guān)鍵字“*DEPVAR”中定義的數(shù)值。狀態(tài)變量矩陣的維數(shù)通過ABAQUS輸入文件中的關(guān)鍵字“*DEPVAR”定義，關(guān)鍵字下面數(shù)據(jù)行的數(shù)值即為狀態(tài)變量矩陣的維數(shù)。材料常數(shù)的個數(shù)，等于關(guān)鍵字“*USER MATERIAL”中“CONSTANTS”常數(shù)設(shè)定的值。材料常數(shù)矩陣，矩陣中元素的數(shù)值對應(yīng)于關(guān)鍵字“*USER MATERIAL”下面的數(shù)據(jù)行。，分別定義每一增量步的彈性應(yīng)變能，塑性耗散和蠕變耗散。它們

17、對計算結(jié)果沒有影響，僅僅作為能量輸出。其他變量：應(yīng)變矩陣：應(yīng)變增量矩陣：增量步的時間增量：直接應(yīng)力分量的個數(shù)：剪切應(yīng)力分量的個數(shù)：總應(yīng)力分量的個數(shù)，使用UMAT時需要注意單元的沙漏控制剛度和橫向剪切剛度。通常減縮積分單元的沙漏控制剛度和板、殼、梁單元的橫向剪切剛度是通過材料屬性中的彈性性質(zhì)定義的。這些剛度基于材料初始剪切模量的值，通常在材料定義中通過“*ELASTIC”選項(xiàng)定義。但是使用UMAT時，ABAQUS對程序輸入文件進(jìn)行預(yù)處理的時候得不到剪切模量的數(shù)值。所以這時候用戶必須使用“*HOURGLASS STIFFNESS”選項(xiàng)來定義具有沙漏模式的單元的沙漏控制剛度，使用“*TRANSVER

18、SE SHEAR STIFFNESS”選項(xiàng)來定義板、殼、梁單元的橫向剪切剛度。21.3.2 編程基于上面所述的率相關(guān)材料公式和應(yīng)力更新算法，參照ABAQUS用戶材料子程序的接口規(guī)范，進(jìn)行UMAT的編程。有限元模擬結(jié)果將在下一節(jié)給出，在最后一節(jié)中還給出了相應(yīng)的程序源代碼。由于UMAT在單元的積分點(diǎn)上調(diào)用，增量步開始時，主程序路徑將通過UMAT的接口進(jìn)入UMAT，單元當(dāng)前積分點(diǎn)必要變量的初始值將隨之傳遞給UMAT的相應(yīng)變量。在UMAT結(jié)束時，變量的更新值將通過接口返回主程序。整個UMAT的流程如圖21-2所示。圖21-2 UMAT流程圖一共有8個材料常數(shù)需要給定，并申請了一個13維的狀態(tài)變量矩陣，

19、它們表示的物理含義如表21-1所示。表21-1 UMAT材料常數(shù)PROPS12345678物理性質(zhì)楊氏模量泊松比塑性耗散比ABnCMSTATEV1-67-1213變量意義彈性應(yīng)變塑性應(yīng)變等效塑性應(yīng)變下一步將使用建立的UMAT結(jié)合ABAQUS/Standard進(jìn)行分離式Hopkinson壓桿（SHPB）實(shí)驗(yàn)的有限元模擬，并對結(jié)果進(jìn)行比較。21.4 SHPB實(shí)驗(yàn)的有限元模擬下面將建立SHPB實(shí)驗(yàn)的有限元模型，并把前面所建立的UMAT接入ABAQUS/Standard進(jìn)行有限元模擬。進(jìn)行有限元模擬的目的不僅是為了再現(xiàn)SHPB實(shí)驗(yàn)的過程，同時也是為了對選擇的數(shù)值模型和建立的UMAT進(jìn)行評價。21.4.

20、1 分離式Hopkinson壓桿（SHPB）實(shí)驗(yàn)分離式Hopkinson壓桿（Split Hopkinson Pressure Bar，簡稱SHPB）實(shí)驗(yàn)是從經(jīng)典Hopkinson實(shí)驗(yàn)基礎(chǔ)之上發(fā)展而來的一種實(shí)驗(yàn)技術(shù)，用來測量材料的動態(tài)應(yīng)力應(yīng)變行為。該實(shí)驗(yàn)技術(shù)的理論基礎(chǔ)是一維應(yīng)力波理論，通過測量兩根壓桿上的應(yīng)變來推導(dǎo)試件上的應(yīng)力應(yīng)變關(guān)系。分離式Hopkinson壓桿實(shí)驗(yàn)示意圖見圖21-3。圖21-3 分離式Hopkinson壓桿裝置21.4.2 有限元建模有限元模型主要是參照前面介紹的SHPB實(shí)驗(yàn)裝置，通過載荷、邊界條件等的定義在有限元中模擬SHPB實(shí)驗(yàn)的環(huán)境，盡量在較少的機(jī)時耗費(fèi)下達(dá)到更高的精

21、度。模型的簡化與有限元網(wǎng)格為了不使模型過于龐大，對模型進(jìn)行了一些簡化。首先，改變?nèi)肓U和出力桿的尺寸，長度由原來的3040mm減小為1000mm，直徑增加到25mm，試件的長度和直徑也分別變化為22mm和18mm。這樣不僅優(yōu)化了網(wǎng)格的質(zhì)量，還成倍地減小了模型的規(guī)模，其帶來的負(fù)面影響就是試件能達(dá)到的應(yīng)變將降低。另外，由于撞擊桿僅僅起到產(chǎn)生應(yīng)力脈沖的作用，在數(shù)值模型中沒必要考慮撞擊桿，取代的方法是直接在入力桿的輸入端施加均布的應(yīng)力脈沖?？紤]到實(shí)驗(yàn)裝置的對稱性，也做了一些簡化。整個實(shí)驗(yàn)裝置以及載荷等都是關(guān)于桿的中心線軸對稱的，所以可以使用軸對稱單元進(jìn)行二維分析。另外也建立了四分之一橫截面的三維模型作

22、為補(bǔ)充。二維軸對稱模型和三維模型分別如圖21-4和圖21-5所示。在模型中，對試件以及入力桿，出力桿和試件接觸的部分進(jìn)行了局部網(wǎng)格加密，這樣的網(wǎng)格劃分可以取得比較經(jīng)濟(jì)的結(jié)果。圖21-4 二維軸對稱有限元模型圖21-5 三維有限元模型單元類型上，選擇一階常規(guī)單元，由于沒有使用減縮積分單元，所以使用UMAT時無需指定單元的沙漏控制剛度。最后的模型中，二維網(wǎng)格單元總數(shù)為1220，三維模型網(wǎng)格的單元總數(shù)為17160。關(guān)于單元的詳細(xì)信息如表21-2所示。表21-2 模型信息模 型尺寸mm(×) 單元類型單元個數(shù)總節(jié)點(diǎn)數(shù)總單元數(shù)二維模型入力桿25×1000CAX453014751220

23、試件18×22CAX4160出力桿25×1000CAX4530三維模型入力桿25×1000C3D875002104917160試件18×22C3D82160出力桿25×1000C3D87500在二維模型局部網(wǎng)格存在疏密連接的部位，一個單元邊要同時和兩個單元連接，這在通常的有限元網(wǎng)格中是不好實(shí)現(xiàn)的。本例在這里使用了ABAQUS的多點(diǎn)約束技術(shù)（Multi-Points Constrain，簡稱MPC）來解決，線性多點(diǎn)約束方式如圖21-6所示。圖21-6 線性多點(diǎn)約束技術(shù)圖中點(diǎn)使用線性多點(diǎn)約束后，其節(jié)點(diǎn)自由度均由旁邊的節(jié)點(diǎn)和線性插值得到。所以使用多點(diǎn)

24、約束方式可以很好連接模型中網(wǎng)格疏密不同的部位，劃分出比較精練的網(wǎng)格。材料定義入力桿和出力桿使用線彈性材料，彈性模量和泊松比分別為200GPa和0.3，密度為7.85×103 kg/m3。試件采用用戶在UMAT中的自定義材料，材料參數(shù)如表21-3所示，其中Johnson-Cook模型中參數(shù)的數(shù)值來源于對實(shí)驗(yàn)數(shù)據(jù)的擬合。表21-3 試件的材料定義性質(zhì)密度Kg/m3楊氏模量MPa泊松比Johnson-Cook模型參數(shù)AMPaBMPanCM數(shù)值2.7×10368.0×1030.3366.562108.8530.2380.0290.5邊界條件為了保證SHPB實(shí)驗(yàn)的要求，在二

25、維模型和三維模型中均施加了必要的邊界條件。在對稱軸或?qū)ΨQ面上施加了對稱性邊界條件，同時保證壓桿和試件可以沿軸線方向自由無約束的運(yùn)動。壓桿和試件之間的接觸為硬接觸，光滑無摩擦。為了確定輸入應(yīng)力脈沖的時間，進(jìn)行了簡單的計算。彈性材料中縱波波速的計算公式為：(21-26)其中為材料彈性模量，為材料密度。由此可以計算輸入應(yīng)力波在壓桿中的傳播速度為m/s。要求在入力桿應(yīng)力波的輸入端不能出現(xiàn)入射波和反射波的重疊，也就是說在輸入應(yīng)力脈沖的時間內(nèi)，應(yīng)力波的傳播距離不應(yīng)超過兩倍的桿長，即：(21-27)根據(jù)這一估計，選擇輸入應(yīng)力脈沖的持續(xù)時間s，上升時間s。經(jīng)過若干次試算，對輸入應(yīng)力脈沖的波形進(jìn)行了適當(dāng)?shù)恼{(diào)整，

26、使試件中產(chǎn)生較均勻的應(yīng)變率。最后輸入應(yīng)力脈沖的波形如圖21-7所示：為了確定增量步的最大時間步長，需要先簡單計算一下單元的穩(wěn)定極限。基于一個單元的估算，穩(wěn)定極限可以用單元特征長度和材料波速定義如下：(21-28)圖21-7 輸入應(yīng)力脈沖壓桿單元的特征單元長度Le = 10 mm，由此可以計算出應(yīng)力波在壓桿傳遞的穩(wěn)定極限為 (s)(21-29)將它作為ABAQUS自動增量控制里面的最大時間步長。21.4.3 二維動態(tài)分析下面將討論二維動態(tài)分析的結(jié)果。為了便于比較，進(jìn)行了三類的分析，首先是無溫度影響時的強(qiáng)化模型，然后是考慮溫度影響的強(qiáng)化模型，最后是考慮單元失效的模型。無溫度影響強(qiáng)化模型分析所進(jìn)行的

27、SHPB實(shí)驗(yàn)正是屬于這一情況，可以將ABAQUS/Standard結(jié)合UMAT進(jìn)行有限元模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比。下面是應(yīng)變率250 s-1下的動態(tài)模擬過程。在時間左右，應(yīng)力波前沿到達(dá)試件，這一時間和前面使用彈性波波速計算的傳播時間是相同的，此前試件上的Mises應(yīng)力幾乎為零，如圖21-8所示。圖21-8 應(yīng)力波前沿到達(dá)試件時的Mises應(yīng)力 (t=1.98×10-4 s)在時間，試件經(jīng)過應(yīng)力波的上升時間后達(dá)到穩(wěn)定變形的狀態(tài)，一部分入射波反射回入力桿，一部分應(yīng)力波經(jīng)過試件進(jìn)入出力桿，試件各點(diǎn)的變形都很均勻，如圖21-9 (a) 所示。在圖21-9 (b) 試件的放大圖上可以看出，各

28、點(diǎn)Mises應(yīng)力相差不超過1MPa，這個精度是相當(dāng)可靠的。(a) 全局視圖(b) 試件的放大視圖圖21-9 試件經(jīng)歷均勻變形時的Mises應(yīng)力 (t=3.0×10-4 s)經(jīng)過穩(wěn)定變形階段后，反射波和傳遞波分別向入力桿和出力桿擴(kuò)散，試件上Mises應(yīng)力逐漸減小到較低的水平，試件開始經(jīng)歷卸載，如圖21-10所示。圖中Mises應(yīng)力云圖的單位為kPa，如不作特別說明，下面各種應(yīng)力云圖中應(yīng)力單位都為kPa。圖21-10 應(yīng)力波消退后試件時的Mises應(yīng)力 (t=4.2×10-4 s)在離壓桿兩端0.2m處各取一個單元作為輸出，其應(yīng)力歷史如圖21-11所示，從中可直觀地看出壓桿上應(yīng)

29、力的傳播過程。圖21-11 壓桿上的應(yīng)力輸出（實(shí)際輸出）取出試件同一橫截面的三個單元以及試件表面長度方向的三個單元，將不同點(diǎn)的應(yīng)力應(yīng)變歷史比較如下：圖21-12 試件橫截面應(yīng)變、應(yīng)力及應(yīng)變率歷史比較圖21-13 試件長度方向應(yīng)變、應(yīng)力及應(yīng)變率歷史比較從試件各點(diǎn)的應(yīng)力應(yīng)變分布上看，圖21-12中應(yīng)變、應(yīng)力及應(yīng)變率歷史曲線基本重合，同一橫截面內(nèi)各點(diǎn)的變化歷史基本一致。圖21-13中試件長度方向各點(diǎn)的變形歷史也基本一致，除了初始階段和最后卸載階段曲線出現(xiàn)較小波動外，其他階段曲線也基本重合。出現(xiàn)小幅度波動的原因可能是應(yīng)力波剛剛到達(dá)試件不能立即達(dá)到平衡狀態(tài)，另外應(yīng)力波經(jīng)歷試件的彈塑性變形后沿其他方向也會

30、出現(xiàn)擴(kuò)散?？梢哉J(rèn)為總體上試件的變形是均勻的，試件各點(diǎn)的應(yīng)變率基本保持在250 s-1。其他三種應(yīng)變率下（分別為70、100、200 s-1）的有限元模擬結(jié)果與250 s-1的情況大體相同，試件的變形均勻，各點(diǎn)的應(yīng)變率基本保持在設(shè)計的水平。這給研究恒定應(yīng)變率下的應(yīng)力應(yīng)變曲線提供了有力的前提條件。實(shí)際上有限元模擬的應(yīng)力應(yīng)變曲線和恒定應(yīng)變率下實(shí)驗(yàn)的結(jié)果也能夠很好的吻合。取出試件表面中間的一點(diǎn)，將應(yīng)變率250 s-1和200 s-1下ABAQUS有限元模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比見圖21-14和圖21-15。四種應(yīng)變率（分別為70、100、200、250 s-1）下使用ABAQUS進(jìn)行有限元模擬結(jié)果的對比見

31、圖21-16到圖21-20。圖21-14 應(yīng)力應(yīng)變曲線的對比及模擬過程中真實(shí)應(yīng)變率變化 (250 s-1)圖21-15 應(yīng)力應(yīng)變曲線的對比及模擬過程中真實(shí)應(yīng)變率變化 (200 s-1)圖21-16 四種應(yīng)變率下的應(yīng)變時間曲線圖21-17 四種應(yīng)變率下的應(yīng)力時間曲線圖21-18 四種應(yīng)變率下的應(yīng)力應(yīng)變曲線圖21-19 四種設(shè)計應(yīng)變率下的真實(shí)應(yīng)變率時間曲線圖21-20 四種應(yīng)變率下的應(yīng)變率應(yīng)變曲線在不同的應(yīng)變率下，材料表現(xiàn)出不同的應(yīng)力應(yīng)變關(guān)系，體現(xiàn)出本構(gòu)模型中應(yīng)變率的作用。而且應(yīng)力應(yīng)變水平相對于應(yīng)變率變化較大，在較高的應(yīng)變率作用下，試件發(fā)生了較大的變形，承受的應(yīng)力也較大。以上結(jié)果說明，選用John

32、son-Cook強(qiáng)化模型能很好地反映不同應(yīng)變率作用下材料的力學(xué)行為，對于金屬材料的沖擊問題是適用的。在此基礎(chǔ)上建立的用戶材料子程序是可靠的?？紤]溫度影響的強(qiáng)化模型分析下面將在UMAT中考慮塑性耗散向熱量轉(zhuǎn)化所引起的溫度的改變，從而影響本構(gòu)模型中的應(yīng)力應(yīng)變關(guān)系。這一過程是非耦合的，看作是絕熱過程，即熱量在單元中產(chǎn)生，引起單元溫度的變化，但是相鄰單元之間沒有熱量的傳遞。在應(yīng)變率250 s1下試件最終的溫度場如圖21-21所示。(a) t=1.98×10-4 s時的溫度場(b) t=3.0×10-4 s時的溫度場(c) t=3.0×10-4 s時的溫度場圖21-21 試

33、件上溫度場的變化取塑性耗散比為0.95，即95的塑性耗散將轉(zhuǎn)化為熱量。這一比例對于大多數(shù)金屬是合適的。從圖21-21中看出，早期應(yīng)力波沒有到達(dá)試件時，試件上的溫度沒有變化，經(jīng)歷塑性變形后，塑性耗散開始向熱量轉(zhuǎn)化，引起單元溫度的變化。試件上的溫升較為均勻，變形結(jié)束時試件上的最高溫升約為1.9。為了使溫度對屈服應(yīng)力的影響明顯，取Johnson-Cook模型中參數(shù)（通常對于大多數(shù)金屬值在1.0左右），這相當(dāng)于夸大了模型中的溫度軟化效應(yīng)。將考慮溫度效應(yīng)時單元的應(yīng)力應(yīng)變與前面沒有考慮溫度效應(yīng)的情況做一下比較，試件表面中間單元的溫度曲線如圖21-22所示。曲線中發(fā)生塑性變形后，單元的溫度逐漸上升。圖21-

34、22 單元溫度的變化在溫度的影響下單元的應(yīng)力、應(yīng)變歷史曲線如圖21-23所示，并同前面不考慮溫度作用的情況進(jìn)行了對比。圖中隨著單元溫度的升高，應(yīng)力較前面不考慮溫度的情況要略低，而應(yīng)變則略大。這說明單元表現(xiàn)出溫度軟化的特征。應(yīng)力應(yīng)變曲線表示如下（圖21-24），在單元進(jìn)入塑性后，溫度升高的同時，屈服強(qiáng)度和硬化率要比前面不考慮溫度的情況要小。Johnson-Cook模型中最后一項(xiàng)關(guān)于溫度軟化的情況在有限元模擬中得到了驗(yàn)證，同時也說明UMAT能正確地考慮溫度對材料本構(gòu)關(guān)系的影響。(a) 應(yīng)變時間曲線(b) 應(yīng)力時間曲線圖21-23 單元應(yīng)力、應(yīng)變歷史圖21-24 考慮溫度效應(yīng)與不考慮溫度效應(yīng)的比較2

35、1.4.4 三維動態(tài)分析本例題進(jìn)行的三維動態(tài)分析也相當(dāng)成功，試件在模擬過程中各處應(yīng)力都比較均勻，最大應(yīng)力和最小應(yīng)力的差別不超過1MPa。圖21-25為試件的Mises應(yīng)力分布。圖21-25 三維有限元模擬中試件的Mises應(yīng)力分布變形過程中單元的應(yīng)力、應(yīng)變以及應(yīng)變率歷史都與二維的情況非常吻合，如圖21-26所示。圖21-26 三維有限元模擬的應(yīng)力、應(yīng)變率歷史下面是應(yīng)力應(yīng)變曲線：圖21-27 三維有限元模擬中的應(yīng)力應(yīng)變曲線這是對編寫的UMAT用于三維實(shí)體單元的一個驗(yàn)證。UMAT雖然是從一維Johnson-Cook模型中建立起來，但是在UMAT中率相關(guān)塑性公式將它擴(kuò)展到了三維情況，所以能用于三維實(shí)

36、體單元也是意料之中。雖然取得的結(jié)果是一致的，但是計算時間上卻有很大差別，在個人計算機(jī)上，使用軸對稱二維網(wǎng)格完成一次計算只需要10分鐘，然后完成一次三維網(wǎng)格的計算需要3個小時。相比而言使用軸對稱二維模型要經(jīng)濟(jì)得多。21.5 UMAT的Fortran程序21.5.1 UMAT SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS, 3 NSTATV,PRO

37、PS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)C INCLUDE 'ABA_PARAM.INC'C CHARACTER*80 MATERL DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROP

38、S),COORDS(3),DROT(3,3), 4 DFGRD0(3,3),DFGRD1(3,3)C DIMENSION EELAS(6),EPLAS(6),FLOW(6) PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0, HALF =0.5d0) DATA NEWTON,TOLER/40,1.D-6/CC -C UMAT FOR JOHNSON-COOK MODELC -C PROPS(1) - YANG'S MODULUSC PROPS(2) - POISSON RATIOC PROPS(3) - INELASTIC HEA

39、T FRACTIONC PARAMETERS OF JOHNSON-COOK MODEL:C PROPS(4) - AC PROPS(5) - B C PROPS(6) - nC PROPS(7) - CC PROPS(8) - mC -C IF (NDI.NE.3) THEN WRITE(6,1) 1 FORMAT(/,30X,'*ERROR - THIS UMAT MAY ONLY BE USED FOR ', 1 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS') ENDIFCC ELASTIC PROPERTI

40、ESC EMOD=PROPS(1) ENU=PROPS(2) IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499 EBULK3=EMOD/(ONE-TWO*ENU) EG2=EMOD/(ONE+ENU) EG=EG2/TWO EG3=THREE*EG ELAM=(EBULK3-EG2)/THREECC ELASTIC STIFFNESSC DO 20 K1=1,NTENS DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.0 10 CONTINUE 20 CONTINUEC DO 40 K1=1,NDI DO 30 K2=1,NDI DD

41、SDDE(K2,K1)=ELAM 30 CONTINUE DDSDDE(K1,K1)=EG2+ELAM 40 CONTINUE DO 50 K1=NDI+1,NTENS DDSDDE(K1,K1)=EG 50 CONTINUECC CALCULATE STRESS FROM ELASTIC STRAINSC DO 70 K1=1,NTENS DO 60 K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 60 CONTINUE 70 CONTINUECC RECOVER ELASTIC AND PLASTIC STRAINSC D

42、O 80 K1=1,NTENS EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS) 80 CONTINUE EQPLAS=STATEV(1+2*NTENS)CC CALCULATE MISES STRESSC IF(NPROPS.GT.5.AND.PROPS(4).GT.0.0) THEN SMISES=(STRESS(1)-STRESS(2)*(STRESS(1)-STRESS(2) + 1 (STRESS(2)-STRESS(3)*(STRESS(2)-STRESS(3) + 1 (STRESS(3)-STRESS(1)*(

43、STRESS(3)-STRESS(1) DO 90 K1=NDI+1,NTENS SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 90 CONTINUE SMISES=SQRT(SMISES/TWO)CC CALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESSCC CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4)C DETERMINE IF ACTIVELY YIELDINGC IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THENCC M

44、ATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTIONC SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3)/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) 110 CONTINUE DO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUECCREAD PARAMETERS OF JOHNSON-COOK MODELC A=PROPS(4) B=PROPS(

45、5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8)CC NEWTON ITERATIONC SYIELD=SYIEL0 DEQPL=(SMISES-SYIELD)/EG3 DSTRES=TOLER*SYIEL0/EG3 DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTON DEQPL=MAX(DEQPL,DEQMIN) CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4) TVP=C*LOG(DEQPL/DTIME) TVP1=TVP+ONE HARD1=HARD*TVP1

46、+SYIELD*C/DEQPL SYIELD=SYIELD*TVP1 RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1) IF(ABS(RHS/EG3) .LE. DSTRES ) GOTO 140 130 CONTINUE WRITE(6,2) NEWTON 2 FORMAT(/,30X,'*WARNING - PLASTICITY ALGORITHM DID NOT ', 1 'CONVERGE AFTER ',I3,' ITERATIONS') 140 CONTINUE EFFHR

47、D=EG3*HARD1/(EG3+HARD1)CC CALCULATE STRESS AND UPDATE STRAINSC DO 150 K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO 150 CONTINUE DO 160 K1=NDI+1,NTENS STRESS(K1)=FLOW(K1)*SYIELD EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL 160 CONTINUE EQPLAS=EQPLAS+DEQPL SPD=DEQPL*(SYIEL0+SYIELD)/TWO RPL = PROPS(3)*SPD/DTIMECC JAC