一起學習戶定義材料User

一起學習UMATZHANGchunyu1、 什么時候用用戶定義材料(User-definedmaterial,UMAT)?很簡單，當ABAQUS沒有提供我們需要的材料模型時。所以，在決定自己定義一種新的材料模型之前，最好對ABAQUS已經(jīng)提供的模型心中有數(shù)，并且盡量使用現(xiàn)有的模型，因為這些模型已經(jīng)經(jīng)過詳細的驗證，并被廣泛接受。2、 好學嗎？需要哪些基礎知識？先看一下ABAQUS手冊(ABAQUSAnalysisUser'sManual)里的一段話:Warning: Theuseofthisoptiongenerallyrequiresconsiderableexpertise.Theuseriscautionedthattheimplementationofanyrealisticconstitutivemodelrequiresextensivedevelopmentandtesting.Initialtestingonasingleelementmodelwithprescribedtractionloadingisstronglyrecommended.但這并不意味著非力學專業(yè)，或者力學基礎知識不很豐富者就只能望洋興嘆，因為我們的任務不是開發(fā)一套完整的有限元軟件，而只是提供一個描述材料力學性能的本構方程(Constitutiveequation)而已。當然，最基本的一些概念和知識還是要具備的，比如應力(stress)，應變(strain)及其分量；volumetriepart和deviatoricpart；模量(modulus)、泊松比(Poisson'sratio)、拉美常數(shù)(Lameconstant)；矩陣的加減乘除甚至求逆；還有一些高等數(shù)學知識如積分、微分等。3、 UMAT的基本任務？我們知道，有限元計算(增量方法)的基本問題是：已知第n步的結果(應力，應變等)g,s；然后給出一個應變增量ds,n n n+1計算新的應力。UMAT要完成這一計算，并要計算Jacobian矩陣DDSDDE(I,J)= 。 是應力增量矩陣(張量或許更合適)， 一J是應變增量矩陣。DDSDDE(I,J)定義了第J個應變分量的微小變化對第I個應力分量帶來的變化。該矩陣只影響收斂速度，不影響計算結果的準確性(當然，不收斂自然得不到結果)。4、怎樣建立自己的材料模型?本構方程就是描述材料應力應變(增量)關系的數(shù)學公式，不是憑空想象出來的，而是根據(jù)實驗結果作出的合理歸納。比如對彈性材料，實驗發(fā)現(xiàn)應力和應變同步線性增長，所以用一個簡單的數(shù)學公式描述。為了解釋彈塑性材料的實驗現(xiàn)象，又提出了一些彈塑性模型，并用數(shù)學公式表示出來。對各向同性材料(Isotropicmaterial),經(jīng)常采用的辦法是先研究材料單向應力-應變規(guī)律(如單向拉伸、壓縮試驗)，并用一數(shù)學公式加以描述，然后把講該規(guī)律推廣到各應力分量。這叫做“泛化“(generalization)。5、一個完整的例子及解釋下面這個UMAT取自ABAQUS手冊，是一個用于大變形下的彈塑性材料模型。希望我的注釋能幫助初學者理解。需要了解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-應力矩陣，在增量步的開始，保為并作為已知量傳入UMAT;在增量步的結束應該保存更新的應o；STRAW-當前應變e，已知。DSTRAN—應變增量，已知。STATEV-狀態(tài)變量矩陣，用來保存用戶自己定義的一些變量，如累計塑性應變，粘彈性應變等等。增量步開始時作為已知量傳入，增量步結束應該更新；DDSDD曲△刃&士。需要更新DTIME—時間增量dt。已知。NDI—正應力、應變個數(shù)對三維問題、軸對稱問題自然是3(11,22,33)，平面問題是2(11,22)；已知。NSHR—剪應力、應變個數(shù)，三維問題時3(12,13,23)，軸對稱問題是1(12)；已知。NTEN2NTENS+NSHR已知。PROPS材料常數(shù)矩陣，如模量啊，粘度系數(shù)啊等等；作為已知量傳入，已知。DROT—對finitestrain問題，應變應該排除旋轉部分，該矩陣提供了旋轉矩陣，詳見下面的解釋。已知。PNEWDT—可用來控制時間步的變化。如果設置為小于1的數(shù)，則程序放棄當前計算，并用新的時間增量DTIMEXPNEWDT作為新的時間增量計算；這對時間相關的材料如聚合物等有用；如果設為大余的數(shù)，則下一個增量步加大)TIME為DTIMEXPNEWDTo可以更新。其他變量含義可參看手冊，暫時用不到。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局部變量，用來暫時保存彈性應變、塑性應變分量以及流動方向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矩陣很容易計算「九「九+2G九九 九+2G九 九J=九九九+2GGGG注意，在ABAQUS中，剪切應變采用工程剪切應變的定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讀取彈性應變分量，塑性應變分量，并旋轉(調(diào)用了OTSIG)，分別保存在EELAS和EPLAS中；CALLROTSIG(STATEV(1),DROT,EELAS,2,NDI,NSHR)CALLROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)讀取等效塑性應變EQPLAS=STATEV(1+2*NTENS)先假設沒有發(fā)生塑性流動，按完全彈性變形計算試算應力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計算Mises應力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ù)當前等效塑性應變，調(diào)用IARDSUB得到對應的屈服應力CGETYIELDSTRESSFROMTHESPECIFIEDHARDENINGCURVECNVALUE=NPROPS/2-1CALLHARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)CCDETERMINEIFACTIVELYYIELDINGC 如果Mises應力大余屈服應力，屈服發(fā)生，計算流動方向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理論并應用Newton-Rampson方法求得等效塑性應變增量CSOLVEFOREQUIVALENTVONMISESSTRESSCANDEQUIVALENTPLASTICSTRAININCREMENTUSINGNEWTONITERATIONCSYIELD=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更新應力了 ，應變分量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計算塑性變形下的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將彈性應變，塑性應變分量保存到狀態(tài)變量中，并傳到下一個增量步CSTOREELASTICAND(EQUIVALENT)PLASTICSTRAINSCINSTATEVARIABLEARRAYCDOK1=1,NTENSSTATEV(K1)=EELAS(K1)STATEV(K1+NTENS)=EPLAS(K1)ENDDOSTATEV(1+2*NTENS)=EQPLASRETURNENDc...c...子程序，根據(jù)等效塑性應變，利用插值的