1、材料本構(gòu)模型及編程-ABAQUS-UMAT材料本構(gòu)模型及編程實現(xiàn)：簡介1、什么時候用用戶定義材料（User-defined material, UMAT）？很簡單，當(dāng)ABAQUS沒有提供我們需要的材料模型時。所以，在決定自己定義一種新的材料模型之前，最好對ABAQUS已經(jīng)提供的模型心中有數(shù)，并且盡量使用現(xiàn)有的模型，因為這些模型已經(jīng)經(jīng)過詳細(xì)的驗證，并被廣泛接受。2、好學(xué)嗎？需要哪些基礎(chǔ)知識？先看一下ABAQUS手冊（ABAQUS Analysis User's Manual）里的一段話：Warning: The use

2、 of this option generally requires considerable expertise. The user is cautioned that the implementation of any realistic constitutive model requires extensive development and test

3、ing. Initial testing on a single element model with prescribed traction loading is strongly recommended.但這并不意味著非力學(xué)專業(yè)，或者力學(xué)基礎(chǔ)知識不很豐富者就只能望洋興嘆，因為我們的任務(wù)不是開發(fā)一套完整的有限元軟件，而只是提供一個描述材料力學(xué)性能的本構(gòu)方程（Constitutive equation）而已。當(dāng)然，最基本的一些概念和知識還是要具備

4、的，比如應(yīng)力(stress),應(yīng)變（strain）及其分量； volumetric part和deviatoric part；模量（modulus）、泊松比(Poissons ratio)、拉美常數(shù)(Lame constant)；矩陣的加減乘除甚至求逆；還有一些高等數(shù)學(xué)知識如積分、微分等。3、UMAT的基本任務(wù)？ 我們知道，有限元計算（增量方法）的基本問題是： 已知第n步的結(jié)果（應(yīng)力，應(yīng)變等） ，； 然后給出一個應(yīng)變增量, 計算新的應(yīng)力 。 UMAT要完成這一計算，并要計算J

5、acobian矩陣DDSDDE(I,J) =。是應(yīng)力增量矩陣（張量或許更合適）， 是應(yīng)變增量矩陣。DDSDDE(I,J) 定義了第J個應(yīng)變分量的微小變化對第I 個應(yīng)力分量帶來的變化。該矩陣只影響收斂速度，不影響計算結(jié)果的準(zhǔn)確性（當(dāng)然，不收斂自然得不到結(jié)果）。4、怎樣建立自己的材料模型？ 本構(gòu)方程就是描述材料應(yīng)力應(yīng)變（增量）關(guān)系的數(shù)學(xué)公式，不是憑空想象出來的，而是根據(jù)實驗結(jié)果作出的合理歸納。比如對彈性材料，實驗發(fā)現(xiàn)應(yīng)力和應(yīng)變同步線性增長，所以用一個簡單的數(shù)學(xué)公式描述。為了解釋彈塑性材料的實驗現(xiàn)象，又提出了一些彈塑性模型，并用數(shù)學(xué)公式表示出來。&#

6、160;對各向同性材料（Isotropic material）,經(jīng)常采用的辦法是先研究材料單向應(yīng)力-應(yīng)變規(guī)律（如單向拉伸、壓縮試驗），并用一數(shù)學(xué)公式加以描述，然后把講該規(guī)律推廣到各應(yīng)力分量。這叫做“泛化“(generalization)。5、一個完整的例子及解釋  下面這個UMAT取自ABAQUS手冊，是一個用于大變形下的彈塑性材料模型。希望我的注釋能幫助初學(xué)者理解。需要了解J2理論。 SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT, 1 DRPLDE,

7、DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED, 2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT, 3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)STRESS-應(yīng)力矩陣，在增量步的開始，保存并作為已知量傳入UMAT ；在增量步的結(jié)束應(yīng)該保存更新的應(yīng)力；STRAN-當(dāng)前應(yīng)變，已知 。 DSTRAN應(yīng)變增量，已知。STATEV-狀

8、態(tài)變量矩陣，用來保存用戶自己定義的一些變量，如累計塑性應(yīng)變，粘彈性應(yīng)變等等。增量步開始時作為已知量傳入，增量步結(jié)束應(yīng)該更新；DDSDDE=。需要更新DTIME時間增量dt。已知。NDI正應(yīng)力、應(yīng)變個數(shù),對三維問題、軸對稱問題自然是3（11,22,33），平面問題是2(11,22)；已知。NSHR 剪應(yīng)力、應(yīng)變個數(shù)，三維問題時3(12,13,23)，軸對稱問題是1(12)；已知。NTENS=NTENS  NSHR，已知。PROPS材料常數(shù)矩陣，如模量啊，粘度系數(shù)啊等等；作為已知量傳入，已知。DROT對finite strain問題，應(yīng)變應(yīng)該排除旋轉(zhuǎn)部分，該矩陣提供了

9、旋轉(zhuǎn)矩陣，詳見下面的解釋。已知。PNEWDT可用來控制時間步的變化。如果設(shè)置為小于1的數(shù)，則程序放棄當(dāng)前計算，并用新的時間增量DTIME X PNEWDT作為新的時間增量計算；這對時間相關(guān)的材料如聚合物等有用；如果設(shè)為大余1的數(shù)，則下一個增量步加大DTIME為DTIME X PNEWDT。可以更新。其他變量含義可參看手冊，暫時用不到。C INCLUDE 'ABA_PARAM.INC'定義了一些參數(shù)，變量什么的，不用管C CHARACTER*8 CMNAMEC DIMENSION 

10、STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS), 1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS), 2 PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3), 3 DFGRD0(3,3),DFGRD1(3,3)矩陣的尺寸聲明CC LOCAL ARRAYSC -C EELAS - ELASTIC STR

11、AINSC EPLAS - PLASTIC STRAINSC FLOW - DIRECTION OF PLASTIC FLOWC -C局部變量，用來暫時保存彈性應(yīng)變、塑性應(yīng)變分量以及流動方向 DIMENSION EELAS(6),EPLAS(6),FLOW(6)C PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0, 1 ENUMAX=.4999D0,NEWTON=10,T

12、OLER=1.0D-6)CC -C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITYC CANNOT BE USED FOR PLANE STRESSC -C PROPS(1) - EC PROPS(2) - NUC PROPS(3.) - SYIELD AN

13、0;HARDENING DATAC CALLS HARDSUB FOR CURVE OF YIELD STRESS VS. PLASTIC STRAINC -CC ELASTIC PROPERTIESC 獲取楊氏模量，泊松比，作為已知量由PROPS向量傳入 EMOD=PROPS(1) E ENU=PROPS(2)   EBULK3=EMOD/(ONE-TWO*ENU) 3K&#

14、160;EG2=EMOD/(ONE ENU) 2G EG=EG2/TWO G EG3=THREE*EG 3G ELAM=(EBULK3-EG2)/THREE  DO K1=1,NTENS DO K2=1,NTENS DDSDDE(K1,K2)=ZERO END DO END DO彈性部分，Jacobian矩陣很容易計算注意，在ABAQUS中，剪切應(yīng)變采用工程剪切應(yīng)變的定義，所以剪切部分模量是G而不是2G!CC ELASTIC&

15、#160;STIFFNESSC DO K1=1,NDI DO K2=1,NDI DDSDDE(K2,K1)=ELAM END DO DDSDDE(K1,K1)=EG2 ELAM END DO DO K1=NDI 1,NTENS DDSDDE(K1,K1)=EG END DOCC RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE

16、0;FORWARDC ALSO RECOVER EQUIVALENT PLASTIC STRAINC讀取彈性應(yīng)變分量，塑性應(yīng)變分量，并旋轉(zhuǎn)（調(diào)用了ROTSIG），分別保存在EELAS和EPLAS中；  CALL ROTSIG(STATEV( 1),DROT,EELAS,2,NDI,NSHR) CALL ROTSIG(STATEV(NTENS 1),DROT,EPLAS,2,NDI,NSHR)讀取等效塑性應(yīng)變  EQPLAS=STATEV(1 2*NTENS)先假設(shè)沒

17、有發(fā)生塑性流動，按完全彈性變形計算試算應(yīng)力CC CALCULATE PREDICTOR STRESS AND ELASTIC STRAINC DO K1=1,NTENS DO K2=1,NTENS STRESS(K2)=STRESS(K2) DDSDDE(K2,K1)*DSTRAN(K1) END DO EELAS(K1)=EELAS(K1) DSTRAN(K1) END DOC計算Mises應(yīng)力C CALCULATE

18、0;EQUIVALENT VON MISES STRESSC SMISES=(STRESS(1)-STRESS(2)*2 (STRESS(2)-STRESS(3)*2 1  (STRESS(3)-STRESS(1)*2 DO K1=NDI 1,NTENS SMISES=SMISES SIX*STRESS(K1)*2 END DO SMISES=SQRT(SMISES/TWO)C 根據(jù)當(dāng)前等效塑性應(yīng)變，調(diào)用HARDSUB得到對應(yīng)的屈服應(yīng)力C GET

19、 YIELD STRESS FROM THE SPECIFIED HARDENING CURVEC NVALUE=NPROPS/2-1 CALL HARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)CC DETERMINE IF ACTIVELY YIELDINGC 如果Mises應(yīng)力大余屈服應(yīng)力，屈服發(fā)生，計算流動方向 IF (SMISES.GT.(ONE TOLER)*SYIEL0)

20、60;THENCC ACTIVELY YIELDINGC SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESSC CALCULATE THE FLOW DIRECTIONC SHYDRO=(STRESS(1) STRESS(2) STRESS(3)/THREE DO K1=1,NDI FLOW(K1)=(STRESS(K1)-SHYDRO)/SMISES END

21、60;DO DO K1=NDI 1,NTENS FLOW(K1)=STRESS(K1)/SMISES END DOC根據(jù)J2理論并應(yīng)用Newton-Rampson方法求得等效塑性應(yīng)變增量C SOLVE FOR EQUIVALENT VON MISES STRESSC AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATIONC SY

22、IELD=SYIEL0 DEQPL=ZERO DO KEWTON=1,NEWTON RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL RHS/(EG3 HARD) CALL HARDSUB(SYIELD,HARD,EQPLAS DEQPL,PROPS(3),NVALUE) IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10 END DOCC WRITE WARNING MESSAGE 

23、TO THE .MSG FILEC WRITE(7,2) NEWTON 2 FORMAT(/,30X,'*WARNING - PLASTICITY ALGORITHM DID NOT ', 1 'CONVERGE AFTER ',I3,' ITERATIONS') 10 CONTINUEC更新應(yīng)力，應(yīng)變分量C UPDATE STR

24、ESS, ELASTIC AND PLASTIC STRAINS AND C EQUIVALENT PLASTIC STRAINC DO K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD SHYDRO EPLAS(K1)=EPLAS(K1) THREE/TWO*FLOW(K1)*DEQPL EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL END DO DO

25、0;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 END DO EQPLAS=EQPLAS DEQPLCC CALCULATE PLASTIC DISSIPATIONC SPD=DEQPL*(SYIEL0 SYIELD)/TWOCC 計算塑性變形下的Jacobian矩陣 

26、60;FORMULATE THE JACOBIAN (MATERIAL TANGENT)C FIRST CALCULATE EFFECTIVE MODULIC EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG EFFG3=THREE/TWO*EFFG2 EFFLAM=(EBULK3-EFFG2)/THREE EFFHRD=EG3*HARD/(EG3 HARD)-EFFG3c. if (props(7).lt.001)

27、60;go to 99c. DO K1=1,NDI DO K2=1,NDI DDSDDE(K2,K1)=EFFLAM END DO DDSDDE(K1,K1)=EFFG2 EFFLAM END DO DO K1=NDI 1,NTENS DDSDDE(K1,K1)=EFFG END DO DO K1=1,NTENS DO K2=1,NTENS DDSDDE(K2,K1)=DDSD

28、DE(K2,K1) EFFHRD*FLOW(K2)*FLOW(K1) END DO END DOc. 99 continuec. ENDIFC將彈性應(yīng)變，塑性應(yīng)變分量保存到狀態(tài)變量中，并傳到下一個增量步C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS C IN STATE VARIABLE ARRAYC DO K1=1,NTENS STATEV(K1)=EELAS(K1) STATEV(K1 NTENS)=EPLAS(K1) END DO STATEV(1 2*NTENS)=EQPLASC RETURN ENDc.c.子程序，根據(jù)等效塑性應(yīng)變，利用插值的方法得到對應(yīng)的屈服應(yīng)力 SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)C INCLUDE 'ABA_PARAM.INC'C DIMENSION TABLE(2,NV