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

材料本構(gòu)模型及編程-ABAQUS-UMAT材料本構(gòu)模型及編程實(shí)現(xiàn)：簡(jiǎn)介1、什么時(shí)候用用戶定義材料(User-definedmaterial,UMAT)?很簡(jiǎn)單，當(dāng)ABAQUS沒有提供我們需要的材料模型時(shí)。所以，在決定自己定義一種新的材料模型之前，最好對(duì)ABAQUS已經(jīng)提供的模型心中有數(shù)，并且盡量使用現(xiàn)有的模型，因?yàn)檫@些模型已經(jīng)經(jīng)過詳細(xì)的驗(yàn)證，并被廣泛接受。2、好學(xué)嗎？需要哪些基礎(chǔ)知識(shí)？先看一下ABAQUS手冊(cè)(ABAQUSAnalysisUser'sManual)里的一段話：Warning:Theuseofthisoptiongenerallyrequiresconsiderableexpertise.Theuseriscautionedthattheimplementationofanyrealisticconstitutivemodelrequiresextensivedevelopmentandtesting.Initialtestingonasingleelementmodelwithprescribedtractionloadingisstronglyrecommended.但這并不意味著非力學(xué)專業(yè)，或者力學(xué)基礎(chǔ)知識(shí)不很豐富者就只能望洋興嘆，因?yàn)槲覀兊娜蝿?wù)不是開發(fā)一套完整的有限元軟件，而只是提供一個(gè)描述材料力學(xué)性能的本構(gòu)方程(Constitutiveequation)而已。當(dāng)然，最基本的一些概念和知識(shí)還是要具備的，比如應(yīng)力(stress),應(yīng)變;strain及其分量；volumetricpart和deviatoricpart模量modulus)、泊松比(Poisson'sratio)、拉美常數(shù)(Lameconstant)；矩陣的加減乘除甚至求逆；還有一些高等數(shù)學(xué)知識(shí)如積分、微分等。3、UMAT的基本任務(wù)？我們知道，有限元計(jì)算(增量方法)的基本問題是：已知第n步的結(jié)果(應(yīng)力應(yīng)變等),；然后給出一個(gè)應(yīng)變?cè)隽?計(jì)算新的應(yīng)力。UMAT要完成這一計(jì)算并要計(jì)算Jacobian矩陣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)變同步線性增長(zhǎng)，所以用一個(gè)簡(jiǎn)單的數(shù)學(xué)公式描述。為了解釋彈塑性材料的實(shí)驗(yàn)現(xiàn)象，又提出了一些彈塑性模型，并用數(shù)學(xué)公式表示出來(lái)。對(duì)各向同性材料(Isotropicmaterial)，經(jīng)常采用的辦法是先研究材料單向應(yīng)力-應(yīng)變規(guī)律(如單向拉伸、壓縮試驗(yàn))，并用一數(shù)學(xué)公式加以描述，然后把講該規(guī)律推廣到各應(yīng)力分量。這叫做“泛化”(generalization)。5、一個(gè)完整的例子及解釋下面這個(gè)UMAT取自ABAQUS手冊(cè)，是一個(gè)用于大變形下的彈塑性材料模型。希望我的注釋能幫助初學(xué)者理解。需要了解J2理論。SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,1DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,2CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,3PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)STRESS--應(yīng)力矩陣，在增量步的開始，保存并作為已知量傳入U(xiǎn)MAT；在增量步的結(jié)束應(yīng)該保存更新的應(yīng)力；STRAN--當(dāng)前應(yīng)變，已知。DSTRAN一應(yīng)變?cè)隽?，已知。STATEV--狀態(tài)變量矩陣，用來(lái)保存用戶自己定義的一些變量，如累計(jì)塑性應(yīng)變，粘彈性應(yīng)變等等。增量步開始時(shí)作為已知量傳入，增量步結(jié)束應(yīng)該更新；DDSDDE=。需要更新DTIME—時(shí)間增量dt。已知。NDI—正應(yīng)力、應(yīng)變個(gè)數(shù)，對(duì)三維問題、軸對(duì)稱問題自然是3(11,22,33)，平面問題是2(11,22)；已知。NSHR—剪應(yīng)力、應(yīng)變個(gè)數(shù)，三維問題時(shí)3(12,13,23)，軸對(duì)稱問題是1(12)；已知。NTENS=NTENSNSHR,已知。PROPS材料常數(shù)矩陣，如模量啊，粘度系數(shù)啊等等；作為已知量傳入，已知。DROT一對(duì)finitestrain問題，應(yīng)變應(yīng)該排除旋轉(zhuǎn)部分，該矩陣提供了旋轉(zhuǎn)矩陣，詳見下面的解釋。已知。PNEWDT一可用來(lái)控制時(shí)間步的變化。如果設(shè)置為小于1的數(shù)，則程序放棄當(dāng)前計(jì)算，并用新的時(shí)間增量DTIMEXPNEWDT作為新的時(shí)間增量計(jì)算；這對(duì)時(shí)間相關(guān)的材料如聚合物等有用；如果設(shè)為大余1的數(shù)，則下一個(gè)增量步加大DTIME為DTIMEXPNEWDT?？梢愿?。其他變量含義可參看手冊(cè)，暫時(shí)用不到。CINCLUDE'ABA_PARAM.INC'定義了一些參數(shù)，變量什么的，不用管CCHARACTER*8CMNAMECDIMENSIONSTRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),1DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),2PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),3DFGRD0(3,3),DFGRD1(3,3)矩陣的尺寸聲明CCLOCALARRAYSC CEELAS-ELASTICSTRAINSCEPLAS-PLASTICSTRAINSCFLOW-DIRECTIONOFPLASTICFLOWC C局部變量，用來(lái)暫時(shí)保存彈性應(yīng)變、塑性應(yīng)變分量以及流動(dòng)方向DIMENSIONEELAS(6),EPLAS(6),FLOW(6)CPARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,1ENUMAX=.4999D0,NEWTON=10,TOLER=1.0D-6)CC CUMATFORISOTROPICELASTICITYANDISOTROPICMISESPLASTICITYCCANNOTBEUSEDFORPLANESTRESSC CPROPS⑴-ECPROPS(2)-NUCPROPS(3..)-SYIELDANHARDENINGDATACCALLSHARDSUBFORCURVEOFYIELDSTRESSVS.PLASTICSTRAINC

先假設(shè)沒有發(fā)生塑性流動(dòng)，按完全彈性變形計(jì)算試算應(yīng)力CCCALCULATEPREDICTORSTRESSANDELASTICSTRAINCDOK1=1,NTENSDOK2=1,NTENSSTRESS(K2)=STRESS(K2)DDSDDE(K2,K1)*DSTRAN(K1)ENDDOEELAS(K1)=EELAS(K1)DSTRAN(K1)ENDDOC計(jì)算乂1$6$應(yīng)力CCALCULATEEQUIVALENTVONMISESSTRESSCSMISES=(STRESS(1)-STRESS(2))**2(STRESS(2)-STRESS(3))**21(STRESS(3)-STRESS(1))**2DOK1=NDI1,NTENSSMISES=SMISESSIX*STRESS(K1)**2ENDDOSMISES=SQRT(SMISES/TWO)C根據(jù)當(dāng)前等效塑性應(yīng)變，調(diào)用HARDSUB得到對(duì)應(yīng)的屈服應(yīng)力CGETYIELDSTRESSFROMTHESPECIFIEDHARDENINGCURVECNVALUE=NPROPS/2-1CALLHARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)CCDETERMINEIFACTIVELYYIELDINGC如果Mises應(yīng)力大余屈服應(yīng)力，屈服發(fā)生，計(jì)算流動(dòng)方向IF(SMISES.GT.(ONETOLER)*SYIEL0)THENCCACTIVELYYIELDINGCSEPARATETHEHYDROSTATICFROMTHEDEVIATORICSTRESSCCALCULATETHEFLOWDIRECTIONCSHYDRO=(STRESS(1)STRESS(2)STRESS(3))/THREEDOK1=1,NDIFLOW(K1)=(STRESS(K1)-SHYDRO)/SMISESENDDODOK1=NDI1,NTENSFLOW(K1)=STRESS(K1)/SMISESENDDOC根據(jù)J2理論并應(yīng)用Newton-Rampson方法求得等效塑性應(yīng)變?cè)隽緾SOLVEFOREQUIVALENTVONMISESSTRESSCANDEQUIVALENTPLASTICSTRAININCREMENTUSINGNEWTONITERATIONCSYIELD=SYIEL0DEQPL=ZERODOKEWTON=1,NEWTONRHS=SMISES-EG3*DEQPL-SYIELDDEQPL=DEQPLRHS/(EG3HARD)CALLHARDSUB(SYIELD,HARD,EQPLASDEQPL,PROPS(3),NVALUE)IF(ABS(RHS).LT.TOLER*SYIEL0)GOTO10ENDDOCCWRITEWARNINGMESSAGETOTHE.MSGFILECWRITE(7,2)NEWTON2FORMAT(//,30X,'***WARNING-PLASTICITYALGORITHMDIDNOT',1'CONVERGEAFTER',I3,'ITERATIONS')10CONTINUEC更新應(yīng)力，應(yīng)變分量CUPDATESTRESS,ELASTICANDPLASTICSTRAINSANDCEQUIVALENTPLASTICSTRAINCDOK1=1,NDISTRESS(K1)=FLOW(K1)*SYIELDSHYDROEPLAS(K1)=EPLAS(K1)THREE/TWO*FLOW(K1)*DEQPLEELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPLENDDODOK1=NDI1,NTENSSTRESS(K1)=FLOW(K1)*SYIELDEPLAS(K1)=EPLAS(K1)THREE*FLOW(K1)*DEQPLEELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPLENDDOEQPLAS=EQPLASDEQPLCCCALCULATEPLASTICDISSIPATIONCSPD=DEQPL*(SYIEL0SYIELD)/TWO

STATEV(12*NTENS)=EQPLASRETURNENDc...子程序，根據(jù)等效塑性應(yīng)變，利用插值的方法得到對(duì)應(yīng)的屈服應(yīng)力SUBROUTINEHARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE)INCLUDE'ABA_PARAM.INC'DIMENSIONTABLE(2,NVALUE)PARAMETER(ZERO=0.D0)CSETYIELDSTRESSTOLASTVALUEOFTABLE,HARDENINGTOZEROSYIELD=TABLE(LNVALUE)HARD=ZEROCIFMORETHANONEENTRY,SEARCHTABLEIF(NVALUE.GT.1)THENDOK1=1,NVALUE-1EQPL1=TABLE(2,K11)IF(EQPLAS.LT.EQPL1)THENEQPL