【正文】 兩者對(duì)于塑性階段應(yīng)力更新的算法不同。常剛度法采用的是：而切線剛度法采用的是：把握好了兩者的本質(zhì)上的區(qū)別，對(duì)于兩者的算法設(shè)計(jì)和程序開發(fā)問題便迎刃而解5. UMAT程序設(shè)計(jì)和編碼 本章將嚴(yán)格按照前一章推導(dǎo)的公式展開程序設(shè)計(jì)和編碼，為了便于編程，本文將本構(gòu)關(guān)系做了抽象化處理，即將其描述成一個(gè)含參數(shù)的表達(dá)式，改變參數(shù)即可應(yīng)用于不同的模型，這樣做的好處是能保證程序的復(fù)用性，這也是本文反復(fù)強(qiáng)調(diào)的使用UMAT的原則。. 本構(gòu)關(guān)系描述 本文采用各向同性硬化彈塑性材料，材料參數(shù)如下： (51) 彈性部分：，彈性模量E=200000Mpa，泊松比Mu= 圖51 彈性部分本構(gòu)關(guān)系塑性部分：，為了研究方便，取A=700，B=，C=400圖52 塑性部分本構(gòu)關(guān)系將2個(gè)曲線統(tǒng)一到同一個(gè)坐標(biāo)系（為方便顯示，x軸標(biāo)注時(shí)擴(kuò)大了1000倍）圖53 本構(gòu)關(guān)系由此可以求出兩條線的交點(diǎn)即初始屈服點(diǎn)的應(yīng)力Yield0=400Mpa（注意兩條曲線相差了一個(gè)屈服應(yīng)變，因?yàn)閮烧咂鋵?shí)不是一個(gè)坐標(biāo)系）綜上定義的材料常數(shù)見表51：表51彈性模量E200000泊松比Mu屈服應(yīng)力Yield0400A700BC400注意：上面A,B,C的取值只是為了便于理解和分析本文的材料模型，為了保證程序的通用性，本文的參數(shù)在程序中的A,B,C一律用變量表示。. 常剛度法程序設(shè)計(jì)算法設(shè)計(jì)1：定義程序需要用到的常數(shù)和變量2：讀取ABAQUS定義的材料常數(shù)和狀態(tài)變量(這里只定義了一個(gè)狀態(tài)變量)，材料常數(shù)為，彈性模量E，泊松比Mu，屈服應(yīng)力Yield0，參數(shù)A,B,C,并且計(jì)算出剪切模量G,狀態(tài)變量為等效塑性應(yīng)變EQPLAS3：讀取應(yīng)力分量，計(jì)算平均應(yīng)力，應(yīng)力偏量以及Mises等效應(yīng)力平均應(yīng)力： 應(yīng)力偏量：Mises等效應(yīng)力：4：根據(jù)3計(jì)算的Mises等效應(yīng)力和2讀取的屈服應(yīng)力Yield0比較，如果Mises等效應(yīng)力小于屈服應(yīng)力，表明此時(shí)材料未屈服，那么轉(zhuǎn)到5，否則轉(zhuǎn)到65：雅可比矩陣,初始化為0，計(jì)算彈性矩陣，按照彈性理論更新應(yīng)力6：雅可比矩陣,初始化為01） ：計(jì)算切線模量H39。,注意到當(dāng)?shù)刃苄詰?yīng)時(shí)對(duì)應(yīng)于本構(gòu)關(guān)系的屈服點(diǎn)，此時(shí)的H不能通過(guò)上式計(jì)算，可以取此時(shí)的H為彈性模量2） ：計(jì)算等效塑性應(yīng)變?cè)隽坎⒏翬QPLAS=DEQPLAS+EQPLAS更新狀態(tài)變量3）計(jì)算流動(dòng)方向 4） 計(jì)算塑性應(yīng)變?cè)隽?）更新應(yīng)力算法流程圖圖54 常剛度法算法流程圖. 常剛度法程序編碼根據(jù)算法流程，用FORTRAN77 固定格式編制了常剛度法的計(jì)算程序如下： SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT, 1 DRPLDE,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) INCLUDE 39。39。 CHARACTER*8 CMNAME DIMENSION 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),DPLAS(6) C 常剛度法 CCCCCCCCCCCCCCCCCCCCCCCCCC1：定義變量c 定義常數(shù) PARAMETER(ZERO=,ONE=,TWO=,THREE=)C 定義材料常數(shù) DOUBLE PRECISION E,Mu,Yield0,A,B,CC 定義中間變量 DOUBLE PRECISION H,W,CEGMA_EQ,CEGMA_X,CEGMA_Y,CEGMA_Z,TAO_XY, 1TAO_YZ,TAO_ZX,CEGMA_CP,CEGMA_SX,CEGMA_SY,CEGMA_SZ,TAO_SXY,TAO_SYZ 2,TAO_SZXC 定義狀態(tài)變量 DOUBLE PRECISION EQPLASC 定義更新變量 DOUBLE PRECISION DEQPLASCCCCCCCCCCCCCCCCCCCCCCCCCC2：讀取變量C 彈性模量 E=PROPS(1)C 泊松比 Mu=PROPS(2)C 剪切模量 G=E/TWO/(ONE+Mu)C 屈服應(yīng)力 Yield0=PROPS(3) C 參數(shù) A=PROPS(4) B=PROPS(5) C=PROPS(6)C 讀取狀態(tài)參數(shù) EQPLAS=STATEV(1) Yield=A*(EQPLAS+Yield0/E)**B+CCCCCCCCCCCCCCCCCCCCCCCCCCC3：讀取應(yīng)力分量，計(jì)算平均應(yīng)力，應(yīng)力偏量以及Mises等效應(yīng)力C 6個(gè)應(yīng)力分量 CEGMA_X= STRESS(1) CEGMA_Y= STRESS(2) CEGMA_Z= STRESS(3) TAO_XY= STRESS(4) TAO_ZX= STRESS(5) TAO_YZ= STRESS(6)C 平均應(yīng)力 CEGMA_CP=(CEGMA_X+CEGMA_Y+CEGMA_Z)/THREEC 6個(gè)應(yīng)力偏量 CEGMA_SX=CEGMA_XCEGMA_CP CEGMA_SY=CEGMA_YCEGMA_CP CEGMA_SZ=CEGMA_ZCEGMA_CP TAO_SXY=TAO_XY TAO_SZX=TAO_ZX TAO_SYZ=TAO_YZc Mises等效應(yīng)力 CEGMA_EQ=SQRT(THREE/TWO*(CEGMA_SX**2+CEGMA_SY**2+CEGMA_SZ**2+ 1TWO*(TAO_SXY**2+TAO_SYZ**2+TAO_SZX**2))) C 雅可比矩陣,初始化默認(rèn)為 DO k1=1,NTENS DO k2=1,NTENS DDSDDE(K1,K2)=ZERO END DO END DOc 計(jì)算彈性矩陣 DO K1=1,NDI DO K2=1,NDI DDSDDE(K1,k2)=E*(ONEMu)/(ONE+Mu)/(ONETWO*Mu)*(Mu/(ONEMu)) END DO DDSDDE(K1,K1)=E*(ONEMu)/(ONE+Mu)/(ONETWO*Mu) END DO DO K1=NDI+1,NTENS DDSDDE(K1,K1)=E*(ONEMu)/(ONE+Mu)/(ONETWO*Mu)*(ONETWO*Mu)/TWO/ 1(ONEMu) END DO CCCCCCCCCCCCCCCCCCCCCCCCCC4：根據(jù)計(jì)算的Mises等效應(yīng)力和讀取的屈服應(yīng)力Yield0比較，CCCCCCCCCCCCCCCCCCCCCCCCCC如果Mises等效應(yīng)力小于屈服應(yīng)力，表明此時(shí)材料未屈服，那么轉(zhuǎn)到，否則轉(zhuǎn)到 IF() THENCCCCCCCCCCCCCCCCCCCCCCCCCC5：沒有屈服，按照彈性理論計(jì)算C 按照彈性理論更新應(yīng)力 DO K1=1,NTENS DO K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) END DO END DO ELSECCCCCCCCCCCCCCCCCCCCCCCCCC6：屈服發(fā)生，按照塑性理論計(jì)算C 切線模量H，根據(jù)本構(gòu)關(guān)系求導(dǎo) IF() THEN H=E ELSE H=A*B*EQPLAS**(BONE) END IFC 求W W=*G/TWO/CEGMA_EQ**2/(H+THREE*G)C 等效塑性應(yīng)變?cè)隽?DEQPLAS=(*G*CEGMA_EQ/(TWO*CEGMA_EQ**2*H+*G*(CEGMA_SX 1**2+CEGMA_SY**2+CEGMA_SZ**2+TWO*TAO_XY**2+TWO*TAO_ZX**2+TWO*TAO_YZ 2**2)))*(CEGMA_SX*DSTRAN(1)+CEGMA_SY*DSTRAN(2)+CEGMA_SZ* 3DSTRAN(3)+TAO_XY*DSTRAN(4)+TAO_ZX*DSTRAN(5)+TAO_YZ*DSTRAN(6))C write(10,*) DEQPLAS,DEQPLAS,EQPLAS,HC 更新狀態(tài)變量 STATEV(1)=EQPLAS+DEQPLASC 計(jì)算塑性應(yīng)變?cè)隽? DPLAS(1)= DEQPLAS*THREE*CEGMA_SX/TWO/CEGMA_EQ DPLAS(2)= DEQPLAS*THREE*CEGMA_SY/TWO/CEGMA_EQ DPLAS(3)= DEQPLAS*THREE*CEGMA_SZ/TWO/CEGMA_EQ DPLAS(4)= DEQPLAS*THREE*TAO_XY/CEGMA_EQ DPLAS(5)= DEQPLAS*THREE*TAO_ZX/CEGMA_EQ DPLAS(6)= DEQPLAS*THREE*TAO_ZY/CEGMA_EQ C 按照塑性理論更新應(yīng)力 DO K1=1,NTENS DO K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*(DSTRAN(K1)DPLAS(K1)) END DO END DO END IF CCCCCCCCCCCCCCCCCCCCCCCCCC計(jì)算完成 RETURN END. 切線剛度法程序設(shè)計(jì)算法設(shè)計(jì)1：定義程序需要用到的常數(shù)和變量2：讀取ABAQUS定義的材料常數(shù)和狀態(tài)變量(這里只定義了一個(gè)狀態(tài)變量)，材料常數(shù)為，彈性模量E，泊松比Mu，屈服應(yīng)力Yield0，參數(shù)A,B,C,并且計(jì)算出剪切模量G,狀態(tài)變量為等效塑性應(yīng)變EQPLAS3：讀取應(yīng)力分量，計(jì)算平均應(yīng)力，應(yīng)力偏量以及Mises等效應(yīng)力平均應(yīng)力： 應(yīng)力偏量：Mises等效應(yīng)力：4：根據(jù)3計(jì)算的Mises等效應(yīng)力和2讀取的屈服應(yīng)力Yield0比較，如果Mises等效應(yīng)力小于屈服應(yīng)力，表明此時(shí)材料未屈服，那么轉(zhuǎn)到5，否則轉(zhuǎn)到65：雅可比矩陣,初始化為0，然后計(jì)算彈性矩陣，按照彈性理論更新應(yīng)力6：雅可比矩陣,初始化為01） 計(jì)算切線模量H注意到當(dāng)?shù)刃苄詰?yīng)時(shí)對(duì)應(yīng)于本構(gòu)關(guān)系的屈服點(diǎn)，此時(shí)的H不能通過(guò)上式計(jì)算，可以取此時(shí)的H為彈性模量2） 計(jì)算w根據(jù)前一章推導(dǎo)的公式3） 計(jì)算等效塑性應(yīng)變?cè)隽緿EQPLAS，并更新狀態(tài)變量根據(jù)前一章推導(dǎo)的公式：帶入后可得等效塑性應(yīng)變?cè)隽緿EQPLAS然后EQPLAS=DEQPLAS+EQPLAS更新狀態(tài)變量4） 計(jì)算雅可比矩陣首先初始化默認(rèn)為0，然后用下式計(jì)算雅可比矩陣注意：abaqus的剪應(yīng)力方向跟彈性力學(xué)規(guī)定的方向不一致，所以上式的最后兩行的應(yīng)交換5） 更新應(yīng)力算法流程圖 圖55 切線剛度法算法流程圖. 切線剛度法程序編碼根據(jù)算法流程，用FORTRAN77 固定格式編制了切線剛度法的計(jì)算程序如下： SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT, 1 DRPLDE,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) INCLUDE 39。39。 CHARACTER*8 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS), 1 DDSDDT(N