1、材料本構模型及編程-ABAQUS-UMAT材料本構模型及編程實現：簡介1、什么時候用用戶定義材料（User-defined material, UMAT）？很簡單，當ABAQUS沒有提供我們需要的材料模型時。所以，在決定自己定義一種新的材料模型之前，最好對ABAQUS已經提供的模型心中有數，并且盡量使用現有的模型，因為這些模型已經經過詳細的驗證，并被廣泛接受。2、好學嗎？需要哪些基礎知識？先看一下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.但這并不意味著非力學專業，或者力學基礎知識不很豐富者就只能望洋興嘆，因為我們的任務不是開發一套完整的有限元軟件，而只是提供一個描述材料力學性能的本構方程（Constitutive equation）而已。當然，最基本的一些概念和知識還是要具備

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

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

6、160;對各向同性材料（Isotropic material）,經常采用的辦法是先研究材料單向應力-應變規律（如單向拉伸、壓縮試驗），并用一數學公式加以描述，然后把講該規律推廣到各應力分量。這叫做“泛化“(generalization)。5、一個完整的例子及解釋  下面這個UMAT取自ABAQUS手冊，是一個用于大變形下的彈塑性材料模型。希望我的注釋能幫助初學者理解。需要了解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-應力矩陣，在增量步的開始，保存并作為已知量傳入UMAT ；在增量步的結束應該保存更新的應力；STRAN-當前應變，已知 。 DSTRAN應變增量，已知。STATEV-狀

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

9、旋轉矩陣，詳見下面的解釋。已知。PNEWDT可用來控制時間步的變化。如果設置為小于1的數，則程序放棄當前計算，并用新的時間增量DTIME X PNEWDT作為新的時間增量計算；這對時間相關的材料如聚合物等有用；如果設為大余1的數，則下一個增量步加大DTIME為DTIME X PNEWDT。可以更新。其他變量含義可參看手冊，暫時用不到。C INCLUDE 'ABA_PARAM.INC'定義了一些參數，變量什么的，不用管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局部變量，用來暫時保存彈性應變、塑性應變分量以及流動方向 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中，剪切應變采用工程剪切應變的定義，所以剪切部分模量是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讀取彈性應變分量，塑性應變分量，并旋轉（調用了ROTSIG），分別保存在EELAS和EPLAS中；  CALL ROTSIG(STATEV( 1),DROT,EELAS,2,NDI,NSHR) CALL ROTSIG(STATEV(NTENS 1),DROT,EPLAS,2,NDI,NSHR)讀取等效塑性應變  EQPLAS=STATEV(1 2*NTENS)先假設沒

17、有發生塑性流動，按完全彈性變形計算試算應力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應力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 根據當前等效塑性應變，調用HARDSUB得到對應的屈服應力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應力大余屈服應力，屈服發生，計算流動方向 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根據J2理論并應用Newton-Rampson方法求得等效塑性應變增量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更新應力，應變分量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將彈性應變，塑性應變分量保存到狀態變量中，并傳到下一個增量步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.子程序，根據等效塑性應變，利用插值的方法得到對應的屈服應力 SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)C INCLUDE 'ABA_PARAM.INC'C DIMENSION TABLE(2,NV