一起學(xué)習(xí)戶定義材料User

一起學(xué)習(xí)UMATZHANGchunyu1、 什么時(shí)候用用戶定義材料(User-definedmaterial,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è)(ABAQUSAnalysisUser'sManual)里的一段話:Warning: Theuseofthisoptiongenerallyrequiresconsiderableexpertise.Theuseriscautionedthattheimplementationofanyrealisticconstitutivemodelrequiresextensivedevelopmentandtesting.Initialtestingonasingleelementmodelwithprescribedtractionloadingisstronglyrecommended.但這并不意味著非力學(xué)專業(yè)，或者力學(xué)基礎(chǔ)知識(shí)不很豐富者就只能望洋興嘆，因?yàn)槲覀兊娜蝿?wù)不是開(kāi)發(fā)一套完整的有限元軟件，而只是提供一個(gè)描述材料力學(xué)性能的本構(gòu)方程(Constitutiveequation)而已。當(dāng)然，最基本的一些概念和知識(shí)還是要具備的，比如應(yīng)力(stress)，應(yīng)變(strain)及其分量；volumetriepart和deviatoricpart；模量(modulus)、泊松比(Poisson'sratio)、拉美常數(shù)(Lameconstant)；矩陣的加減乘除甚至求逆；還有一些高等數(shù)學(xué)知識(shí)如積分、微分等。3、 UMAT的基本任務(wù)？我們知道，有限元計(jì)算(增量方法)的基本問(wèn)題是：已知第n步的結(jié)果(應(yīng)力，應(yīng)變等)g,s；然后給出一個(gè)應(yīng)變?cè)隽縟s,n n n+1計(jì)算新的應(yīng)力。UMAT要完成這一計(jì)算，并要計(jì)算Jacobian矩陣DDSDDE(I,J)= 。 是應(yīng)力增量矩陣(張量或許更合適)， 一J是應(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,DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRDO,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)STRESS-應(yīng)力矩陣，在增量步的開(kāi)始，保為并作為已知量傳入U(xiǎn)MAT;在增量步的結(jié)束應(yīng)該保存更新的應(yīng)o；STRAW-當(dāng)前應(yīng)變e，已知。DSTRAN—應(yīng)變?cè)隽浚阎?。STATEV-狀態(tài)變量矩陣，用來(lái)保存用戶自己定義的一些變量，如累計(jì)塑性應(yīng)變，粘彈性應(yīng)變等等。增量步開(kāi)始時(shí)作為已知量傳入，增量步結(jié)束應(yīng)該更新；DDSDD曲△刃&士。需要更新DTIME—時(shí)間增量dt。已知。NDI—正應(yīng)力、應(yīng)變個(gè)數(shù)對(duì)三維問(wèn)題、軸對(duì)稱問(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ì)稱問(wèn)題是1(12)；已知。NTEN2NTENS+NSHR已知。PROPS材料常數(shù)矩陣，如模量啊，粘度系數(shù)啊等等；作為已知量傳入，已知。DROT—對(duì)finitestrain問(wèn)題，應(yīng)變應(yīng)該排除旋轉(zhuǎn)部分，該矩陣提供了旋轉(zhuǎn)矩陣，詳見(jiàn)下面的解釋。已知。PNEWDT—可用來(lái)控制時(shí)間步的變化。如果設(shè)置為小于1的數(shù)，則程序放棄當(dāng)前計(jì)算，并用新的時(shí)間增量DTIMEXPNEWDT作為新的時(shí)間增量計(jì)算；這對(duì)時(shí)間相關(guān)的材料如聚合物等有用；如果設(shè)為大余的數(shù)，則下一個(gè)增量步加大)TIME為DTIMEXPNEWDTo可以更新。其他變量含義可參看手冊(cè)，暫時(shí)用不到。CINCLUDE'ABAPARAM.INC'定義了一些參數(shù)，變量什么的，不用管CCHARACTER*8CMNAMECDIMENSIONSTRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)矩陣的尺寸聲明CCLOCALARRAYSC C EELAS-ELASTICSTRAINSC EPLAS-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 C PROPS(1)- EC PROPS(2)- NUC PROPS(3..) -SYIELD AN HARDENING DATAC CALLSHARDSUB FOR CURVE OF YIELD STRESSVS.PLASTICSTRAINC CCELASTICPROPERTIESC獲取楊氏模量，泊松比，作為已知量由ROPS向量傳入EMOD=PROPS(1)EENU=PROPS(2)vEBULK3=EMOD/(ONE-TWO*ENU)3KEG2=EMOD/(ONE+ENU)2GEG=EG2/TWO G

EG3=THREE*EG3GELAM=(EBULK3-EG2)/THREE入DOK1=1,NTENSDOK2=1,NTENSDDSDDE(K1,K2)=ZEROENDDOENDDO彈性部分，Jacobian矩陣很容易計(jì)算「九「九+2G九九 九+2G九 九J=九九九+2GGGG注意，在ABAQUS中，剪切應(yīng)變采用工程剪切應(yīng)變的定Y二u+u，所以剪切iji,jj,i部分模量是G而不是2G!CCELASTICSTIFFNESSCDOK1=1,NDIDOK2=1,NDIDDSDDE(K2,K1)=ELAMENDDODDSDDE(K1,K1)=EG2+ELAMENDDODOK1=NDI+1,NTENSDDSDDE(K1,K1)=EGENDDOCCRECOVERELASTICANDPLASTICSTRAINSANDROTATEFORWARDCALSORECOVEREQUIVALENTPLASTICSTRAINC讀取彈性應(yīng)變分量，塑性應(yīng)變分量，并旋轉(zhuǎn)(調(diào)用了OTSIG)，分別保存在EELAS和EPLAS中；CALLROTSIG(STATEV(1),DROT,EELAS,2,NDI,NSHR)CALLROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)讀取等效塑性應(yīng)變EQPLAS=STATEV(1+2*NTENS)先假設(shè)沒(méi)有發(fā)生塑性流動(dòng)，按完全彈性變形計(jì)算試算應(yīng)力Ao=J.AsG=Q+Ann+1 nCCCALCULATEPREDICTORSTRESSANDELASTICSTRAINCDOK1=1,NTENSDOK2=1,NTENSSTRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)ENDDOEELAS(K1)=EELAS(K1)+DSTRAN(K1)ENDDOC計(jì)算Mises應(yīng)力CCALCULATEEQUIVALENTVONMISESSTRESSCSMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**21+(STRESS(3)-STRESS(1))**2DOK1=NDI+1,NTENSSMISES=SMISES+SIX*STRESS(K1)**2ENDDOSMISES=SQRT(SMISES/TWO)C 根據(jù)當(dāng)前等效塑性應(yīng)變，調(diào)用IARDSUB得到對(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.(ONE+TOLER)*SYIEL0)THENCC ACTIVELYYIELDINGC SEPARATETHEHYDROSTATICFROM THEDEVIATORICSTRESSC CALCULATETHEFLOWDIRECTIONCSHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREEDOK1=1,NDIFLOW(K1)=(STRESS(K1)-SHYDRO)/SMISESENDDODOK1=NDI+1,NTENSFLOW(K1)=STRESS(K1)/SMISESENDDOC根據(jù)J2理論并應(yīng)用Newton-Rampson方法求得等效塑性應(yīng)變?cè)隽緾SOLVEFOREQUIVALENTVONMISESSTRESSCANDEQUIVALENTPLASTICSTRAININCREMENTUSINGNEWTONITERATIONCSYIELD=SYIEL0DEQPL=ZERODOKEWTON=1,NEWTONRHS=SMISES-EG3*DEQPL-SYIELDDEQPL=DEQPL+RHS/(EG3+HARD)CALLHARDSUB(SYIELD,HARD,EQPLAS+DEQPL,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)變分量n+1CUPDATESTRESS,ELASTICANDPLASTICSTRAINSANDCEQUIVALENTPLASTICSTRAINCDOK1=1,NDISTRESS(K1)=FLOW(K1)*SYIELD+SHYDROEPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPLEELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPLENDDODOK1=NDI+1,NTENSSTRESS(K1)=FLOW(K1)*SYIELDEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPLEELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPLENDDOEQPLAS=EQPLAS+DEQPLCCCALCULATEPLASTICDISSIPATIONCSPD=DEQPL*(SYIEL0+SYIELD)/TWOCC計(jì)算塑性變形下的Jacobian矩陣FORMULATETHEJACOBIAN(MATERIALTANGENT)CFIRSTCALCULATEEFFECTIVEMODULICEFFG=EG*SYIELD/SMISESEFFG2=TWO*EFFGEFFG3=THREE/TWO*EFFG2EFFLAM=(EBULK3-EFFG2)/THREEEFFHRD=EG3*HARD/(EG3+HARD)-EFFG3c...if(props(7).lt..001)goto99c...DOK1=1,NDIDOK2=1,NDIDDSDDE(K2,K1)=EFFLAMENDDODDSDDE(K1,K1)=EFFG2+EFFLAMENDDODOK1=NDI+1,NTENSDDSDDE(K1,K1)=EFFGENDDODOK1=1,NTENSDOK2=1,NTENSDDSDDE(K2,K1)=DDSDDE(K2,K1)+EFFHRD*FLOW(K2)*FLOW(K1)ENDDOENDDOc...99continuec...ENDIFC將彈性應(yīng)變，塑性應(yīng)變分量保存到狀態(tài)變量中，并傳到下一個(gè)增量步CSTOREELASTICAND(EQUIVALENT)PLASTICSTRAINSCINSTATEVARIABLEARRAYCDOK1=1,NTENSSTATEV(K1)=EELAS(K1)STATEV(K1+NTENS)=EPLAS(K1)ENDDOSTATEV(1+2*NTENS)=EQPLASRETURNENDc...c...子程序，根據(jù)等效塑性應(yīng)變，利用插值的