abaqus_standard用戶材料子程序?qū)嵗齙W

1、ABAQUS 軟件 2003 年度用戶年會論ABAQUS/Standard 用戶材料子程序?qū)嵗?Johnson-Cook 金屬本構(gòu)模型盧劍鋒 莊茁* 張帆 清華大學(xué)工程力學(xué)系 北京 100084摘要：用戶材料子程序是 ABAQUS 提供給用戶定義自己的材料屬性的 Fortran 程序接口，使用戶能使用 ABAQUS 材料庫中沒有定義的材料模型。 ABAQUS 中自有的 Johnson-Cook 模型只能應(yīng)用于顯式 ABAQUS/Explicit 程序中，而我們希望能在隱式 ABAQUS/Standard 程序中更精確的實現(xiàn)本構(gòu)積分，而且應(yīng)用 Johnson-Cook 模型的修正形式。這就需要通

2、過 ABAQUS/Standard 的用戶材料子程序 UMAT 編程實現(xiàn)。在 UMAT 編程中使用了率相關(guān)塑性理論以及完全隱式的應(yīng)力更新算法。 1 Johnson-Cook 強化模型簡介 Johnson-Cook（JC）模型用來模擬高應(yīng)變率下的金屬材料。JC 強化模型表示為三項的乘積， 分別反映了應(yīng)變硬化，應(yīng)變率硬化和溫度軟化。這里使用 JC 模型的修正形式： o = ( A + Be n ) + + e& (1 - T *m )1 C ln 1 e& 0 并使參考應(yīng)變率e&0 = 1 ，這樣公式中的 A 即為材料的靜態(tài)屈服應(yīng)力。公式中包含 A, B, n, C, m 五個參數(shù)，需要通過實驗來

3、確定。 2 ABAQUS 用戶材料子程序 用戶材料子程序（ User-defined Material Mechanical Behavior，簡稱 UMAT）通過與ABAQUS 主求解程序的接口實現(xiàn)與 ABAQUS 的數(shù)據(jù)交流。在輸入文件中，使用關(guān)鍵字“*USER MATERIAL”表示定義用戶材料屬性。 子程序概況與接口 UMAT 子程序具有強大的功能，使用 UMAT 子程序： (1) 可以定義材料的本構(gòu)關(guān)系，使用 ABAQUS 材料庫中沒有包含的材料進行計算，擴充程序功能。 - 18 -(2) 幾乎可以用于力學(xué)行為分析的任何分析過程，幾乎可以把用戶材料屬性賦予 ABAQUS 中的任何單元

4、； (3) 必須在 UMAT 中提供材料本構(gòu)模型的雅可比（Jacobian）矩陣，即應(yīng)力增量對應(yīng)變增量的變化率。 (4) 可以和用戶子程序“USDFLD”聯(lián)合使用，通過“USDFLD”重新定義單元每一物質(zhì)點上傳遞到 UMAT 中場變量的數(shù)值。 由于主程序與 UMAT 之間存在數(shù)據(jù)傳遞，甚至共用一些變量，因此必須遵守有關(guān) UMAT 的書寫格式，UMAT 中常用的變量在文件開頭予以定義，通常格式為： SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,2 STRAN,DSTRAN,TIME,DTI

5、ME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)CINCLUDE ABA_PARAM.INCCCHARACTER*80 CMNAMEDIMENSION STRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2 STRAN(NTENS),DSTRAN(NTE

6、NS),TIME(2),PREDEF(1),DPRED(1),3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDTRETURNENDUMAT 中的應(yīng)力矩陣、應(yīng)變矩陣以及矩陣 DDSDDE ，DDSDDT ，DRPLDE 等，都是直接分 量存儲在前，剪切分量存儲在后。直接分量有 NDI 個，剪切分量有 NS

7、HR 個。各分量之間的順序根據(jù)單元自由度的不同有一些差異，所以編寫 UMAT 時要考慮到所使用單元的類別。下面對UMAT 中用到的一些變量進行說明： DDSDDE ( NTENS, NTENS )是一個 NTENS 維的方陣，稱作雅可比矩陣，D / D ，D 是應(yīng)力的增量，D 是應(yīng)變的增量， DDSDDE ( I , J ) 表示增量步結(jié)束時第 J 個應(yīng)變分量的改變引起的第 I 個應(yīng)力分量的變化。通常雅可比是一個對稱矩陣，除非在“*USER MATERIAL”語句中加入了“UNSYMM”參數(shù)。 STRESS ( NTENS )應(yīng)力張量矩陣，對應(yīng) NDI 個直接分量和 NSHR 個剪切分量。在增

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

9、陣不同的是：對于有限應(yīng)變問題，除了材料本構(gòu)行為引起的數(shù)據(jù)更新以外，狀態(tài)變量矩陣中的任何矢量或者張量都必須通過旋轉(zhuǎ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è)定的值。 PROPS ( NPROPS )材料常數(shù)矩陣，矩陣中元素的數(shù)值對應(yīng)于關(guān)鍵字“*USER MATERIAL”下面的數(shù)據(jù)行。 SSE ， SPD ， SCD分別定義每一增量步的彈性應(yīng)變

10、能，塑性耗散和蠕變耗散。它們對計算結(jié)果沒有影響，僅僅 作為能量輸出。其他變量： STRAN ( NTENS ) ： 應(yīng) 變 矩 陣 ； DSTRAN ( NTENS )：應(yīng)變增量矩陣； DTIME ：增量步的時間增量； NDI ：直接應(yīng)力分量的個數(shù)； NSHR ：剪切應(yīng)力分量的個數(shù)； NTENS ：總應(yīng)力分量的個數(shù)， NTENS = NDI + NSHR 。 使用 UMAT 時需要注意單元的沙漏控制剛度和橫向剪切剛度。通常減縮積分單元的沙漏控制剛度和板、殼、梁單元的橫向剪切剛度是通過材料屬性中的彈性性質(zhì)定義的。這些剛度基于材料初始剪切模量的值，通常在材料定義中通過“*ELASTIC”選項定義。

11、但是使用 UMAT 的時候， ABAQUS 對程序輸入文件進行預(yù)處理的時候得不到剪切模量的數(shù)值。所以這時候用戶必須使用“*HOURGLASS STIFFNESS” 選項來定義具有沙漏模式的單元的沙漏控制剛度，使用“ * T R A N S V E R S E SHEAR STIFFNESS”選項來定義板、殼、梁單元的橫向剪切剛度。 編程 基于上面所述的率相關(guān)材料公式和應(yīng)力更新算法，參照 ABAQUS 用戶材料子程序的接口規(guī)范，進行 UMAT 的編程。有限元模擬結(jié)果將在下一節(jié)給出，在最后一節(jié)中還給出了相應(yīng)的程序源代碼。 由于 UMAT 在單元的積分點上調(diào)用，增量步開始時，主程序路徑將通過 UMA

12、T 的接口進入UMAT，單元當(dāng)前積分點必要變量的初始值將隨之傳遞給 UMAT 的相應(yīng)變量。在 UMAT 結(jié)束時， 變量的更新值將通過接口返回主程序。整個 UMAT 的流程如圖 1 所示。 一共有 8 個材料常數(shù)需要給定，并申請一個 13 維的狀態(tài)變量矩陣，它們表示的物理含義如表 1 所示。 下一步將使用建立的UMAT 結(jié)合ABAQUS/Standard 進行霍布金森沖擊桿(SHPB)實驗的有限元模擬，并對結(jié)果進行比較。 圖 1 UMAT 流程圖 PROPS12345678物理性質(zhì) 楊氏模量 泊松比 塑性耗散比 ABnCMSTATEV1-67-1213變量意義 彈性應(yīng)變 塑性應(yīng)變 等效塑性應(yīng)變

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

14、以使用軸對稱單元進行二維分析。 二維軸對稱模型如圖 2 所示。在模型中，對試件以及入力桿，出力桿和試件接觸的部分進行了局部網(wǎng)格加密，這樣的網(wǎng)格劃分可以取得比較經(jīng)濟的結(jié)果。 模 型 尺寸mm( F L ) 單元類型 單元個數(shù) 總節(jié)點 數(shù) 總單元數(shù) 二維模型 入力桿 251000CAX4530 1475 1220試件 1822CAX4160出力桿 251000CAX4530圖 2二維軸對稱有限元模型表 2 模型信息 材料定義 入力桿和出力桿使用線彈性材料，彈性模量和泊松比分別為 200GPa 和 0.3，密度為 7.85103 kg/m3。試件采用用戶在 UMAT 中自定義材料，材料參數(shù)如表 3

15、所示，其中 Johnson-Cook 模型中參數(shù)的數(shù)值來源于前面的數(shù)值擬合程序。 表 3 試件的材料定義 楊氏模量 MPa 泊松比 Johnson-Cook 模型參數(shù) 性質(zhì) 密度 Kg/m3A MPa B MPa n C M數(shù) 值 2.7103 68.0103 0.3366.562108.8530.238 0.029 0.5邊界條件 為了保證 SHPB 實驗的要求，在二維模型中施加了必要的邊界條件。在對稱軸上施加了對稱性邊界條件，同時保證壓桿和試件可以沿軸線方向自由無約束的運動。壓桿和試件之間的接觸為硬接觸，光滑無摩擦。 為了確定輸入應(yīng)力脈沖的時間，進行了簡單的計算。彈性材料中縱波波速的計算公

16、式為： ErCd =其中 E 為材料彈性模量， r 為材料密度。由此可以計算輸入應(yīng)力波在壓桿中的傳播速度為 Cd = 5048 m/s。 要求在入力桿應(yīng)力波的輸入端不能出現(xiàn)入射波和反射波的重疊，也就是說在輸入應(yīng)力脈沖的時間內(nèi)，應(yīng)力波的傳播距離不應(yīng)超過兩倍的桿長，即： sT 2L =2 4.0 10-4 (s)Cd5048sr根據(jù)這一估計，選擇輸入應(yīng)力脈沖的持續(xù)時間T = 2.0 10-4 s，上升時間t = 3.0 10-5 s。 基 準(zhǔn) 壓 力 (1.7 10-4,200)(3.0 10-5,170)應(yīng) 變 率(s-1)系數(shù)2500.902000.771000.52700.46經(jīng)過若干次試算

17、，對輸入應(yīng)力脈沖的波形進行了適當(dāng)?shù)恼{(diào)整，使試件中產(chǎn)生較均勻的應(yīng)變率。最后輸入應(yīng)力脈沖的波形如圖 3 所示： 300250應(yīng) 力 (MPa)2001501005000.00000.0001圖 30.00020.00030.00040.0005時 間 (s) 輸入應(yīng)力脈沖 為了確定增量步的最大時間步長，需要先簡單計算一下單元的穩(wěn)定極限?；谝粋€單元的估 d算，穩(wěn)定極限可以用單元特征長度 Le 和材料波速C 定義如下： CDLetstable =d壓桿單元的特征單元長度 Le 10 mm，由此可以計算出應(yīng)力波在壓桿傳遞的穩(wěn)定極限為 stableDt 2.010-6 (s)將它作為 ABAQUS 自動

18、增量控制里面的最大時間步長。二維動態(tài)分析 我們所對照的 SHPB 實驗正是屬于這一情況，所以可以將 ABAQUS/Standard 結(jié)合 UMAT 進行有限元模擬的結(jié)果和實驗數(shù)據(jù)進行對比。下面是應(yīng)變率 250 s-1 下的動態(tài)模擬過程。 在時間t = 1.9810-4 (s) 左右，應(yīng)力波前沿到達試件，這一時間和前面使用彈性波波速計算的傳播時間是相同的，此前試件上的 Mises 應(yīng)力幾乎為零，如圖 4 所示。 圖 4應(yīng)力波前沿到達試件時的 Mises 應(yīng)力 (t=1.9810-4 s)在時間t = 3.0 10-4 (s) ，試件經(jīng)過應(yīng)力波的上升時間后達到穩(wěn)定變形的狀態(tài)，一部分入射波 反射回入

19、力桿，一部分應(yīng)力波經(jīng)過試件進入出力桿，試件各點的變形都很均勻，如圖 5 (a)所示。在圖 5 (b)試件的放大圖上可以看出，各點 Mises 應(yīng)力相差不超過 1MPa，這個精度是相當(dāng)可靠的。 (a) 全局視圖 (b) 試件的放大視圖 圖 5 試件經(jīng)歷均勻變形時的 Mises 應(yīng)力 (t=3.010-4 s)經(jīng)過穩(wěn)定變形階段后，反射波和傳遞波分別向入力桿和出力桿擴散，試件上 Mises 應(yīng)力逐漸減小到較低的水平，試件開始經(jīng)歷卸載，如圖 6 所示。圖中 Mises 應(yīng)力云圖的單位為 KPa。 圖 6應(yīng)力波消退后試件時的 Mises 應(yīng)力 (t=4.210-4 s)實際上有限元模擬的應(yīng)力應(yīng)變曲線和恒

20、定應(yīng)變率下實驗的結(jié)果也能夠很好的吻合。取出試件表面中間的一點，將應(yīng)變率 250 s-1 和 200 s-1 下 ABAQUS 有限元模擬的結(jié)果與實驗的結(jié)果對比 見圖 7 和圖 8。 160140120100500實驗值 ABAQUS模擬應(yīng) 變 率450400350應(yīng) 變 率 (s-1)300應(yīng) 力 (MPa)8025060200150401002050000.0000.0050.0100.0150.0200.0250.0300.0350.0400.045應(yīng) 變 圖 7 應(yīng)力應(yīng)變曲線的對比及模擬過程中真實應(yīng)變率變化 (250 s-1)160140120力 (MPa)100500實驗數(shù)據(jù) ABAQ

21、US模擬 應(yīng) 變 率450400350應(yīng) 變 率 (s-1)30080250應(yīng)60200150401002050000.0000.0050.0100.0150.0200.0250.0300.035應(yīng) 變 圖 8 應(yīng)力應(yīng)變曲線的對比及模擬過程中真實應(yīng)變率變化 (200 s-1) 4 J-C 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,NTEN

22、S,3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)CINCLUDE ABA_PARAM.INCCCHARACTER*80 MATERLDIMENSION 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),

23、COORDS(3),DROT(3,3),4 DFGRD0(3,3),DFGRD1(3,3)CDIMENSION 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 CUMAT FOR JOHNSON-COOK MODELC CPROPS(1) - YANGS MODULUSCPROPS(2) - POISSON RATIOCPROPS(3) - INELASTIC HEAT FRACTIONCPARAME

24、TERS OF JOHNSON-COOK MODEL: CPROPS(4) - ACPROPS(5) - BCPROPS(6) - nCPROPS(7) - CCPROPS(8) - mC CIF (NDI.NE.3) THEN WRITE(6,1)1FORMAT(/,30X,*ERROR - THIS UMAT MAY ONLY BE USED FOR ,1ELEMENTS WITH THREE DIRECT STRESS COMPONENTS) ENDIFCCELASTIC PROPERTIES CEMOD=PROPS(1) ENU=PROPS(2)IF(ENU.GT.0.4999.AND

25、.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)/THREECCELASTIC STIFFNESS CDO 20 K1=1,NTENSDO 10 K2=1,NTENS DDSDDE(K2,K1)=0.010CONTINUE20CONTINUE CDO 40 K1=1,NDIDO 30 K2=1,NDI DDSDDE(K2,K1)=ELAM30CONTINUE DDSDDE(K1,K1)=EG2+ELAM40CONTINUE

26、DO 50 K1=NDI+1,NTENSDDSDDE(K1,K1)=EG50CONTINUE CCCALCULATE STRESS FROM ELASTIC STRAINS CDO 70 K1=1,NTENSDO 60 K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)60CONTINUE70CONTINUE CCRECOVER ELASTIC AND PLASTIC STRAINS CDO 80 K1=1,NTENS EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTEN

27、S)80CONTINUEEQPLAS=STATEV(1+2*NTENS)CCCALCULATE MISES STRESS CIF(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)*(STRESS(3)-STRESS(1) DO 90 K1=NDI+1,NTENSSMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)9

28、0CONTINUE SMISES=SQRT(SMISES/TWO)CCCALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESS CCCALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4) CDETERMINE IF ACTIVELY YIELDINGCIF (SMISES.GT.(1.0+TOLER)*SYIEL0) THENCCMATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTION CSHYDRO=(STRESS(1)+STRESS(2)+STR

29、ESS(3)/THREE ONESY=ONE/SMISESDO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)110CONTINUEDO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY120CONTINUE CCREAD PARAMETERS OF JOHNSON-COOK MODEL CA=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8)CCNEWTON ITERATION CSYIELD=SYIEL0DEQPL=(SMISES-SYIELD)

30、/EG3 DSTRES=TOLER*SYIEL0/EG3 DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTONDEQPL=MAX(DEQPL,DEQMIN)CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4) TVP=C*LOG(DEQPL/DTIME)TVP1=TVP+ONE HARD1=HARD*TVP1+SYIELD*C/DEQPL SYIELD=SYIELD*TVP1RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1)IF(ABS(

31、RHS/EG3) .LE. DSTRES ) GOTO 140130CONTINUE WRITE(6,2) NEWTON2FORMAT(/,30X,*WARNING - PLASTICITY ALGORITHM DID NOT ,1CONVERGE AFTER ,I3, ITERATIONS)140CONTINUE EFFHRD=EG3*HARD1/(EG3+HARD1)CCCALCULATE STRESS AND UPDATE STRAINS CDO 150 K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD+SHYDROEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO150CONTINUEDO 160 K1=NDI+1,NTENS