【正文】 程圖. 切線剛度法程序編碼根據(jù)算法流程，用FORTRAN77 固定格式編制了切線剛度法的計算程序如下： 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。 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) 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)力分量，計算平均應(yīng)力，應(yīng)力偏量以及Mises等效應(yīng)力C 6個應(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個應(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)))CCCCCCCCCCCCCCCCCCCCCCCCCC4：根據(jù)3計算的Mises等效應(yīng)力和2讀取的屈服應(yīng)力Yield0比較，CCCCCCCCCCCCCCCCCCCCCCCCCC如果Mises等效應(yīng)力小于屈服應(yīng)力，表明此時材料未屈服，那么轉(zhuǎn)到5，否則轉(zhuǎn)到6 IF() THENCCCCCCCCCCCCCCCCCCCCCCCCCC5：沒有屈服，按照彈性理論計算C 雅可比矩陣,初始化默認為0 DO k1=1,NTENS DO k2=1,NTENS DDSDDE(K1,K2)=ZERO END DO END DOc 計算彈性矩陣 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 DOC 按照彈性理論更新應(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ā)生，按照塑性理論計算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)變增量 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+ 2TWO*TAO_YZ**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 更新狀態(tài)變量 STATEV(1)=EQPLAS+DEQPLASC 雅可比矩陣,初始化默認為0 DO k1=1,NTENS DO k2=1,NTENS DDSDDE(K1,K2)=ZERO END DO END DOC 雅可比矩陣更新 DDSDDE(1,1)=E/(ONE+Mu)*((ONEMu)/(ONETWO*Mu)W*CEGMA_SX**2) DDSDDE(2,1)=E/(ONE+Mu)*((Mu)/(ONETWO*Mu)W*CEGMA_SX*CEGMA_Y) DDSDDE(3,1)=E/(ONE+Mu)*((Mu)/(ONETWO*Mu)W*CEGMA_SX*CEGMA_Z) DDSDDE(4,1)=E/(ONE+Mu)*(W*CEGMA_SZ*TAO_SXY) DDSDDE(5,1)=E/(ONE+Mu)*(W*CEGMA_SZ*TAO_SZX) DDSDDE(6,1)=E/(ONE+Mu)*(W*CEGMA_SZ*TAO_SYZ) DDSDDE(2,2)=E/(ONE+Mu)*((ONEMu)/(ONETWO*Mu)W*CEGMA_Y**2) DDSDDE(3,2)=E/(ONE+Mu)*((Mu)/(ONETWO*Mu)W*CEGMA_SY*CEGMA_SZ) DDSDDE(4,2)=E/(ONE+Mu)*(W*CEGMA_SY*TAO_SXY) DDSDDE(5,2)=E/(ONE+Mu)*(W*CEGMA_SY*TAO_SZX) DDSDDE(6,2)=E/(ONE+Mu)*(W*CEGMA_SY*TAO_SYZ) DDSDDE(3,3)=E/(ONE+Mu)*((ONEMu)/(ONETWO*Mu)W*CEGMA_SZ**2) DDSDDE(4,3)=E/(ONE+Mu)*(W*CEGMA_SZ*TAO_SXY) DDSDDE(5,3)=E/(ONE+Mu)*(W*CEGMA_SZ*TAO_SZX) DDSDDE(6,3)=E/(ONE+Mu)*(W*CEGMA_SZ*TAO_SYZ) DDSDDE(4,4)=E/(ONE+Mu)*(ONE/TWOW*TAO_SXY**2) DDSDDE(5,4)=E/(ONE+Mu)*(W*TAO_SXY*TAO_SZX) DDSDDE(6,4)=E/(ONE+Mu)*(W*TAO_SXY*TAO_SYZ) DDSDDE(5,5)=E/(ONE+Mu)*(ONE/TWOW*TAO_SZX**2) DDSDDE(6,5)=E/(ONE+Mu)*(W*TAO_SYZ*TAO_SZX) DDSDDE(6,6)=E/(ONE+Mu)*(ONE/TWOW*TAO_SYZ**2) DDSDDE(1,2)=DDSDDE(2,1) DDSDDE(1,3)=DDSDDE(3,1) DDSDDE(1,4)=DDSDDE(4,1) DDSDDE(1,5)=DDSDDE(5,1) DDSDDE(1,6)=DDSDDE(6,1) DDSDDE(2,3)=DDSDDE(3,2) DDSDDE(2,4)=DDSDDE(4,2) DDSDDE(2,5)=DDSDDE(5,2) DDSDDE(2,6)=DDSDDE(6,2) DDSDDE(3,4)=DDSDDE(4,3) DDSDDE(3,5)=DDSDDE(5,3) DDSDDE(3,6)=DDSDDE(6,3) DDSDDE(4,5)=DDSDDE(5,4) DDSDDE(4,6)=DDSDDE(6,4) DDSDDE(5,6)=DDSDDE(6,5)C 按照塑性理論更新應(yīng)力 DO K1=1,NTENS DO K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) END DO END DO END IF CCCCCCCCCCCCCCCCCCCCCCCCCC計算完成 RETURN END. 程序的調(diào)試 程序設(shè)計是否合理，程序編碼是否正確并不容易從程序的代碼中簡單的看出，而是通過在調(diào)試中發(fā)現(xiàn)的，可以說，程序的調(diào)試同程序編碼同等重要，通常我們在做UMAT程序的時候，是寫好程序直接到ABAQUS提交，這樣的方法本身并沒有什么錯誤，但是從程序調(diào)試的角度看，這是一種不太好的方法，因為這樣做出錯的可能性極大，不管是語法錯誤還是算法錯誤，一旦出錯，我們很難發(fā)現(xiàn)到底錯誤在哪里，到底為什么會出錯，所以在這種情況下，我們迫切需要采用一種好的調(diào)試手段，以能快速發(fā)現(xiàn)程序中的BUG，修正設(shè)計，下面介紹一下本文采用的調(diào)試方法：調(diào)試方法主要分為2個方面（1） 程序語法檢查（2） 程序算法檢查 對于語法檢查主要是通過編譯器檢查錯誤的方式進行的，但是事實上，如果直接將UMAT子程序拿來給FORTRAN編譯器編譯，不管是否有語法錯誤，都是不能編譯通過的，這是因為每個UMAT程序都有一句INCLUDE 39。但是事實上，翻遍ABAQUS的目錄，我們都不可能找到這個文件，然后復(fù)制到UMAT子程序的同一個目錄下，就可以編譯通過了，通過這種方法，凡是語法錯誤，我們可以在編譯階段發(fā)現(xiàn)和修改，避免在ABAQUS里反復(fù)提交的繁瑣手段。 通過以上的方法，程序錯誤就可以在調(diào)試中快速解決。（2）由于ABAQUS自帶的彈塑性材料本構(gòu)關(guān)系是經(jīng)過大量工程實踐檢驗并被認為是一種正確的材料模型，所以對這個模型采用ABAQUS自帶的彈塑性材料本構(gòu)關(guān)系來計算，得出其Mises等效應(yīng)力值的范圍以及應(yīng)力應(yīng)變關(guān)系，該結(jié)果作為一個真實正確的解。（4）采用切線剛度算法的UMAT子程序計算，得出其Mises等效應(yīng)力值的范圍以及應(yīng)力應(yīng)變關(guān)系；然后同ABAQUS計算的相應(yīng)結(jié)果做比較，驗證切線剛度法的正確性。. 問題描述 用ABAQUS建立一個左端固結(jié)的懸臂梁尺寸為20x20x80，簡單拉伸，右端作用一個均布力420MPa，材料采用前面所述的內(nèi)容，彈性模量200000MPa，如圖61所示 圖 61 驗證模型圖. 本構(gòu)關(guān)系 采用Mises材料模型，屈服應(yīng)力400MPa，塑性階段的應(yīng)力應(yīng)變關(guān)系為前文所述：，其中A=700 B= C=400. ABAQUS自帶材料模型計算 塑性階段應(yīng)力應(yīng)變關(guān)系用WPS表格軟件描點得到，然后采用表格形式輸入到AB