【正文】 TENS),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 屈服應力 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：讀取應力分量，計算平均應力，應力偏量以及Mises等效應力C 6個應力分量 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 平均應力 CEGMA_CP=(CEGMA_X+CEGMA_Y+CEGMA_Z)/THREEC 6個應力偏量 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等效應力 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等效應力和2讀取的屈服應力Yield0比較，CCCCCCCCCCCCCCCCCCCCCCCCCC如果Mises等效應力小于屈服應力，表明此時材料未屈服，那么轉到5，否則轉到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 按照彈性理論更新應力 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ù)本構關系求導 IF() THEN H=E ELSE H=A*B*EQPLAS**(BONE) END IFC 求W W=*G/TWO/CEGMA_EQ**2/(H+THREE*G)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+ 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 按照塑性理論更新應力 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)試 程序設計是否合理，程序編碼是否正確并不容易從程序的代碼中簡單的看出，而是通過在調(diào)試中發(fā)現(xiàn)的，可以說，程序的調(diào)試同程序編碼同等重要，通常我們在做UMAT程序的時候，是寫好程序直接到ABAQUS提交，這樣的方法本身并沒有什么錯誤，但是從程序調(diào)試的角度看，這是一種不太好的方法，因為這樣做出錯的可能性極大，不管是語法錯誤還是算法錯誤，一旦出錯，我們很難發(fā)現(xiàn)到底錯誤在哪里，到底為什么會出錯，所以在這種情況下，我們迫切需要采用一種好的調(diào)試手段，以能快速發(fā)現(xiàn)程序中的BUG，修正設計，下面介紹一下本文采用的調(diào)試方法：調(diào)試方法主要分為2個方面（1） 程序語法檢查（2） 程序算法檢查 對于語法檢查主要是通過編譯器檢查錯誤的方式進行的，但是事實上，如果直接將UMAT子程序拿來給FORTRAN編譯器編譯，不管是否有語法錯誤，都是不能編譯通過的，這是因為每個UMAT程序都有一句INCLUDE 39。39。，但是事實上，翻遍ABAQUS的目錄，我們都不可能找到這個文件，，然后復制到UMAT子程序的同一個目錄下，就可以編譯通過了，通過這種方法，凡是語法錯誤，我們可以在編譯階段發(fā)現(xiàn)和修改，避免在ABAQUS里反復提交的繁瑣手段。 對于程序沒有語法錯誤卻有算法錯誤，這種情況也是很常見的，對程序的算法做檢查，本文采用的主要方式是通過輸出變量到文本的方法，具體就是在程序內(nèi)用OPEN語句打開一個文本文件，然后將程序中關心的中間變量用WRITE語句輸?shù)皆撐谋?，提交給ABAQUS以后，ABAQUS執(zhí)行到該處就會將其輸出到文本，這樣，只需要查看最后的文本文件，就可以發(fā)現(xiàn)算法中是否純在錯誤，修正算法。 通過以上的方法，程序錯誤就可以在調(diào)試中快速解決。6. 程序驗證 任何有關材料模型的開發(fā)都需要經(jīng)過實踐的檢驗和驗證，尤其對于工程問題，通過驗證可以發(fā)現(xiàn)程序中的BUG，并做相應的修改，本文的驗證主要采用綜合比較法，整個驗證過程如下：（1）確定一個具體的模型和它的本構關系。（2）由于ABAQUS自帶的彈塑性材料本構關系是經(jīng)過大量工程實踐檢驗并被認為是一種正確的材料模型，所以對這個模型采用ABAQUS自帶的彈塑性材料本構關系來計算，得出其Mises等效應力值的范圍以及應力應變關系，該結果作為一個真實正確的解。（3）采用常剛度算法的UMAT子程序計算，得出其Mises等效應力值的范圍以及應力應變關系，然后同ABAQUS計算的相應結果做比較，看兩者是否一致，驗證常剛度法的正確性。（4）采用切線剛度算法的UMAT子程序計算，得出其Mises等效應力值的范圍以及應力應變關系；然后同ABAQUS計算的相應結果做比較，驗證切線剛度法的正確性。（5）提高非線性的程度，分別用常剛度法和切線剛度法計算，比較兩者的效率，得出相應的結論。. 問題描述 用ABAQUS建立一個左端固結的懸臂梁尺寸為20x20x80，簡單拉伸，右端作用一個均布力420MPa，材料采用前面所述的內(nèi)容，彈性模量200000MPa，如圖61所示 圖 61 驗證模型圖. 本構關系 采用Mises材料模型，屈服應力400MPa，塑性階段的應力應變關系為前文所述：，其中A=700 B= C=400. ABAQUS自帶材料模型計算 塑性階段應力應變關系用WPS表格軟件描點得到，然后采用表格形式輸入到A