ABAQUS-材料本構(gòu)模型及編程

1、材料本構(gòu)模型及編程-ABAQUS-UMAT材料本構(gòu)模型及編程實(shí)現(xiàn)：簡(jiǎn)介1、什么時(shí)候用用戶(hù)定義材料（User-defined material, UMAT）？很簡(jiǎn)單，當(dāng)ABAQUS沒(méi)有提供我們需要的材料模型時(shí)。所以，在決定自己定義一種新的材料模型之前，最好對(duì)ABAQUS已經(jīng)提供的模型心中有數(shù)，并且盡量使用現(xiàn)有的模型，因?yàn)檫@些模型已經(jīng)經(jīng)過(guò)詳細(xì)的驗(yàn)證，并被廣泛接受。2、好學(xué)嗎？需要哪些基礎(chǔ)知識(shí)？先看一下ABAQUS手冊(cè)（ABAQUS Analysis User's Manual）里的一段話(huà)：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é)專(zhuān)業(yè)，或者力學(xué)基礎(chǔ)知識(shí)不很豐富者就只能望洋興嘆，因?yàn)槲覀兊娜蝿?wù)不是開(kāi)發(fā)一套完整的有限元軟件，而只是提供一個(gè)描述材料力學(xué)性能的本構(gòu)方程（Constitutive equation）而已。當(dāng)然，最基本的一些概念和知識(shí)還是要具備

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

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

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

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

9、旋轉(zhuǎn)矩陣，詳見(jiàn)下面的解釋。已知。PNEWDT可用來(lái)控制時(shí)間步的變化。如果設(shè)置為小于1的數(shù)，則程序放棄當(dāng)前計(jì)算，并用新的時(shí)間增量DTIME X PNEWDT作為新的時(shí)間增量計(jì)算；這對(duì)時(shí)間相關(guān)的材料如聚合物等有用；如果設(shè)為大余1的數(shù)，則下一個(gè)增量步加大DTIME為DTIME X PNEWDT?？梢愿?。其他變量含義可參看手冊(cè)，暫時(shí)用不到。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局部變量，用來(lái)暫時(shí)保存彈性應(yīng)變、塑性應(yīng)變分量以及流動(dò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矩陣很容易計(jì)算注意，在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è)沒(méi)

17、有發(fā)生塑性流動(dòng)，按完全彈性變形計(jì)算試算應(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計(jì)算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得到對(duì)應(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ā)生，計(jì)算流動(dòng)方向 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 計(jì)算塑性變形下的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)變量中，并傳到下一個(gè)增量步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)變，利用插值的方法得到對(duì)應(yīng)的屈服應(yīng)力 SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)C INCLUDE 'ABA_PARAM.INC'C DIMENSION TABLE(2,NV