1、UMAT子程序在復合材料強度分析中的應用本例使用UMAT用戶子程序進行復合材料單層板的應力分析和漸進損傷壓縮強度分析，介紹UMAT用戶子程序編寫方法及在Abaqus/CAE中的設置。本章使用最大應變強度理論作為復合材料單層板的失效準則，相應的Fortran程序簡單易讀，便于理解UAMT子程序的工作原理。知識要點： 強度分析 UMAT用戶子程序 最大應變理論 剛度折減講師：孔祥宏版本：Abq 6.14難度：關鍵詞：強度分析，UMAT&.1 本章內容簡介本章通過兩個實例介紹UMAT用戶子程序在復合材料單層板的應力分析和強度分析中的應用。在第一個實例中，對一個簡單的復合材料單層板進行應力分析，UMA

2、T子程序主要計算應力，不進行強度分析，本例用于驗證UMAT子程序的計算精度。在第二個實例中，對復合材料單層板進行漸進損傷強度分析，UMAT子程序用于應力計算、強度分析和剛度折減。本章所用復合材料為T700/BA9916，材料屬性如表&-1所示。表&-1 T700/BA9916材料屬性參數值強度值E1/GPa114XT/MPa2688E2/GPa8.61XC/MPa1458E3/GPa8.61YT/MPa69.5120.3YC/MPa236130.3ZT/MPa55.5230.45ZC/MPa175G12/GPa4.16SXY/MPa136G13/GPa4.16SXZ/MPa136G23/GPa

3、3.0SYZ/MPa95.6&.2 實例一：UMAT用戶子程序應力分析在使用UMAT用戶子程序進行高級應用之前，應該先了解UMAT子程序，熟悉UMAT子程序的工作原理，了解UMAT中的參數、變量的含義。為了便于讀者快速了解和使用UMAT，本例通過復合材料單層板的應力分析來介紹一個簡單的UMAT子程序。讀者可將本例中的單層板替換為層壓板，進行對比分析。&.2.1問題描述復合材料單層板幾何尺寸為15mm10mm0.15mm，纖維方向為45，單層板的3D實體模型如圖&-1所示，X軸方向為0方向，左側面施加X軸向對稱邊界條件，下側面施加Y軸向對稱邊界條件，垂直于Z軸且Z=0的平面施加Z軸向對稱邊界條件

4、，右側面施加100MPa的拉力。圖&-1 單層板邊界條件及加載情況本例中單位系統為mm、MPa。&.2.2 UMAT用戶子程序本例使用的UMAT用戶子程序UMAT-Stress.for的全部代碼如下，字母C及“!”之后為注釋內容。1 SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,2 1 RPL,DDSDDT,DRPLDE,DRPLDT,3 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,4 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,

5、DROT,PNEWDT,5 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)6 C 7 INCLUDE ABA_PARAM.INC8 C 9 CHARACTER*80 CMNAME10 DIMENSION STRESS(NTENS),STATEV(NSTATV),11 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),12 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),13 3 PROPS(NPROPS),COORDS(

6、3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),14 4 JSTEP(4)15 16 DIMENSION EG(6), XNU(3,3), STRAND(6), C(6,6), STRESS0(6)17 C*18 C EG.E1,E2,E3,G12,G13,G2319 C XNU.NU12，NU21,NU13，NU31，NU23，NU3220 C STRAND.STRAINT AT THE END OF THE INCREMENT21 C C.6X6 STIFFNESS MATRIX22 C STRESS0.STRESS AT THE BEGINNING OF THE

7、 INCREMENT23 C*24 C INITIALIZE XNU & C MATRIX25 XNU=026 C=027 C GET THE MATERIAL PROPERTIES-ENGINEERING CONSTANTS28 EG(1) = PROPS(1) !E1,YOUNGS MODULUS IN DIRECTION 129 EG(2) = PROPS(2) !E2,YOUNGS MODULUS IN DIRECTION 230 EG(3) = EG(2) !E3,YOUNGS MODULUS IN DIRECTION 331 XNU(1,2) = PROPS(3) !POISONS

8、 RATIO POI_1232 XNU(2,1) = XNU(1,2)*EG(2)/EG(1) !POISONS RATIO POI_2133 XNU(1,3) = XNU(1,2) !POISONS RATIO POI_1334 XNU(3,1) = XNU(1,3)*EG(3)/EG(1) !POISONS RATIO POI_3135 XNU(2,3) = PROPS(4) !POISONS RATIO POI_2336 XNU(3,2) = XNU(2,3)*EG(3)/EG(2) !POISONS RATIO POI_3237 EG(4) = PROPS(5) !G12,SHEAR

9、MODULUS IN 12 PLANE38 EG(5) = EG(4) !G13,SHEAR MODULUS IN 13 PLANE39 EG(6) = PROPS(6) !G23,SHEAR MODULUS IN 23 PLANE40 C*41 C FILL THE 6X6 STIFFNESS MATRIX C(6,6)42 RNU = 1/(1-XNU(1,2)*XNU(2,1)-XNU(1,3)*XNU(3,1)-43 1 XNU(3,2)*XNU(2,3)-2*XNU(1,3)*XNU(2,1)*XNU(3,2)44 C STIFFNESS MATRIX C(6,6)45 C(1,1)

10、 = EG(1)*(1-XNU(2,3)*XNU(3,2)*RNU46 C(2,2) = EG(2)*(1-XNU(1,3)*XNU(3,1)*RNU47 C(3,3) = EG(3)*(1-XNU(1,2)*XNU(2,1)*RNU48 C(4,4) = EG(4)49 C(5,5) = EG(5)50 C(6,6) = EG(6)51 C(1,2) = EG(1)*(XNU(2,1)+XNU(3,1)*XNU(2,3)*RNU52 C(2,1) = C(1,2)53 C(1,3) = EG(1)*(XNU(3,1)+XNU(2,1)*XNU(3,2)*RNU54 C(3,1) = C(1,

11、3)55 C(2,3) = EG(2)*(XNU(3,2)+XNU(1,2)*XNU(3,1)*RNU56 C(3,2) = C(2,3)57 C*58 C CALCULATE STRAIN59 DO I = 1, 6 60 STRAND(I) = STRAN(I)+DSTRAN(I)61 ENDDO62 C CALCULATE STRESS63 DO I = 1, 664 STRESS0(I) = STRESS(I)65 STRESS(I) = 066 DO J = 1, 667 STRESS(I) = STRESS(I)+C(I,J)*STRAND(J)68 ENDDO69 ENDDO70

12、 C CALCULATE SSE71 DO I = 1, 672 SSE = SSE+0.5*(STRESS0(I)+STRESS(I)*DSTRAN(I)73 ENDDO74 C*75 C UPDATE DDSDDE76 DO I = 1, 677 DO J = 1, 678 DDSDDE(I,J) = C(I,J)79 ENDDO80 ENDDO81 RETURN82 END第1到14行及第81、82行為UMAT子程序固定格式，其中，第1到5行括號內的變量為UMAT子程序中可以使用的變量，第10到14行定義各變量數組的維數和長度。部分主要變量的含義如表&-2所示。表&-2 UMAT部分變量

13、名及其含義STRESS增量步開始時的應力(S11, S22, .)，用增量步結束時的應力計算結果對其更新STATEV(NSTATV)狀態變量（狀態變量個數），如果在材料中定義了狀態變量，則在UMAT中需要對其更新STRAN增量步開始時的應變(E11, E22, .)DSTRAN當前增量步的應變增量(E11, E22, .)NDI, NSHR, NTENS應力、應變的個數，NDI為正應力或正應變的個數，NSHR為剪應力或剪應變的個數，NTENS=NDI+NSHRPROPS, NPROPS材料參數、材料參數的個數DDSDDE雅克比矩陣，SSE, SPD, SCD特定的彈性應變能、塑性耗散、蠕變耗散

14、，只對能量輸出有影響，對其他計算結果無影響，在UMAT中需要對其更新CELENT單元特征長度第15到83行為用戶自己編寫的固定格式的Fortran程序，用于計算剛度矩陣、應力、應變能、雅克比矩陣。由于本例中沒有使用狀態變量，因此不需要更新STATEV，只需要更新STRESS、DDSDDE和SSE即可。第16行定義了5個數組，其中EG、STRAND、STRESS0為一維數組，XNU、C為二維數組。第18到22行為注釋部分，EG存放材料的3個彈性模量和3個剪切模量；STRAND存放當前增量步結束時的應變(E11, E22, .)；STRESS0存放增量步開始時的應力(S11, S22, .)；XN

15、U為33的二維矩陣，存放泊松比12、21、13、31、23、32；C為66的剛度矩陣。第25、26行初始化二維數組XNU、C，使其每個元素都為0。第28到39行，讀取材料常數，計算泊松比。泊松比的計算公式如下。 (&-1)第42到56行，計算剛度矩陣。剛度矩陣的計算公式如下。 (&-2)對于本例所用材料，由于E2=E3，G12=G13，12=13，所以21=31、23=32，可以將式(&-2)化簡后在UMAT中計算剛度矩陣。本例為保證UMAT-Stress.for的可讀性，剛度矩陣的計算沒有化簡，直接按照式(&-2)編寫。第59到61行，計算當前增量步的應變，計算公式如式(&-3)，式中上標表

16、示增量步序號。 (&-3)第63到69行，將當前分析步開始時的應力STRESS的值賦給STRESS0，然后計算當前增量步的應力并賦給STRESS，當前增量步的應力計算如式(&-4)，式中上標表示增量步序號，應力i和應變i為列向量。 (&-4)第71到73行，計算應變能，如式(&-5)所示，式中上標表示增量步序號；下標表示應力、應變增量的分量序號，其中下標為1、2、3表示正應力、正應變增量，4、5、6表示剪應力、剪應變增量。 (&-5)第76到80行為更新雅克比矩陣。由于沒有對剛度矩陣沒變化，應力與應變的關系如式(&-4)所示，所以應力增量與應變增量的關系如式(&-6)所示，所以雅克比矩陣就等于

17、剛度矩陣。 (&-6)第76到80行可簡寫為一行，如下所示：DDSDDE(1:6,1:6) = C(1:6,1:6)，或DDSDDE = CUMAT中的剪應變為工程剪應變。&.2.3復合材料單層板應力分析1、 創建部件及劃分網格l 創建部件在Part模塊，單擊工具區的(Create Part)，在Create Part對話框中，Name后面輸入Lam-C，Modeling Space選擇3D，Type選擇Deformable，在Base Feature區域選擇Solid、Extrusion，Approximate size使用默認的200，單擊Continue.進入繪圖模式。單擊工具區的(Cr

18、eate Lines: Rectangle (4 Lines)，在提示區輸入第1個點的坐標(0,0)后按回車鍵，再輸入第2個點的坐標(15,10)后按回車，再按Esc鍵或單擊鼠標中鍵。單擊提示區的Done或鼠標中鍵，在Edit Base Extrusion對話框Depth后面輸入0.15，單擊OK完成。l 劃分網格在環境欄Module后面選擇Mesh，進入Mesh模塊。環境欄中Object選擇Part: Lam-C。單擊工具區的(Seed Part)，在Global Seeds對話框中Approximate global size后面輸入0.5，單擊OK。單擊工具區的(Mesh Part)，單

19、擊提示區的Yes或鼠標中鍵，完成網格劃分，如圖&-2所示，板的厚度方向只劃分為1層單元。圖&-2 劃分網格單擊工具區的(Assign Element Type)，在Element Type對話框中，選擇依次選擇Standard，Linear，3D Stress，在Hex標簽頁中勾選Reduced integration，在Element Controls區域Hourglass control選擇Enhanced，即選擇C3D8R單元，單擊OK完成。單擊工具區的(Assign Stack Direction)，在視圖區選擇部件平行于X-Y平面的面，單擊鼠標中鍵或提示區的Yes完成。2、 創建材料

20、并給部件賦材料屬性l 創建材料在環境欄Module后面選擇Property，進入Property模塊。單擊工具區的(Create Material)，在Edit Material對話框中，Name后面輸入UMat-T700；單擊GeneralUser Material，在User Material區域中Data區域的Mechanical Constants一欄依次輸入114000， 8610，0.3，0.45，4160，3000，單擊OK完成。單擊工具區的(Create Material)，在Edit Material對話框中，Name后面輸入Mat-T700；單擊MechanicalElas

21、ticityElastic，在Elastic區域中Type選擇Engineering Constants，在Data區域輸入從左到右依次輸入114000，8610，8610，0.3，0.3，0.45，4160，4160，3000，單擊OK完成。材料UMat-T700用于UMAT用戶子程序，Mat-T700用于做對比分析。輸入數據時，每輸入完一行后按回車（Enter）鍵，光標會自動移到下一行。也可以通過右鍵快捷菜單添加或刪除一行。本例UMAT子程序較簡單，不需要使用狀態變量，因此在材料UMat-T700中沒有定義Depvar。l 給部件賦材料屬性單擊工具區的(Create Composite L

22、ayup)，在打開的對話框中Name使用默認名稱，Initial ply count后面輸入1，Element Type選擇Solid，單擊Continue.；在Edit Composite Layup對話框中，Layup Orientation區域的Definition選擇Coordinate system，單擊Definition下一行的(Create Datum CSYS)。在Create Datum CSYS對話框中使用默認名稱Datum csys-1，類型選擇Rectangular，單擊Continue.，在提示區輸入原點坐標(0,0,0)后按回車鍵，再輸入(1,0,0)后按回車鍵，

23、最后輸入(0,1,0)后按回車鍵，單擊Create Datum CSYS對話框的Cancel。在Edit Composite Layup對話框中單擊(Select CSYS.)，在視圖區選擇剛創建的Datum csys-1，Stacking Direction選擇Element direction 3，Rotation axis選擇Axis 3。在Plies標簽頁中，雙擊Region，在視圖區選擇部件后單擊鼠標中鍵；右單擊Material，在快捷菜單中單擊Edit Material.，在Select Material對話框中選擇UMat-T700，單擊OK；右單擊Element Relativ

24、e Thickness，在快捷菜單中單擊Edit Thickness.，在Thickness對話框中Specify Value后面輸入1后單擊OK；在Rotation Angle一欄輸入0；Integration Points使用默認的1，單擊OK完成。在定義復合材料鋪層時，視圖區部件上會顯示鋪層方向，在Edit Composite Layup對話框中的Display標簽頁中可以設置所需顯示的方向，在視圖區部件上白色箭頭及字母S表示Stacking Direction。3、 裝配在環境欄Module后面選擇Assembly，進入Assembly模塊。單擊工具區的(Create Instance

25、)，在Create Instance對話框中選擇Parts: Lam-C，單擊OK完成。4、 創建分析步、設置輸出變量l 創建分析步在環境欄Module后面選擇Step，進入Step模塊。單擊工具區的(Create Step)，在Create Step對話框中，在Initial分析步之后插入Static, General分析步，單擊Continue.；在Edit Step對話框中使用默認設置，單擊OK完成。l 設置輸出變量單擊工具區的(Field Output Manager)，在Field Output Requests Manager對話框中，選中F-Output-1，單擊Edit.；在E

26、dit Field Output Request對話框中，設置如圖&-3所示，輸出整個模型最后一個增量步的S、E、U，單擊OK完成。圖&-3 場輸出變量設置單擊工具區的(History Output Manager)，在History Output Requests Manager對話框中，選中H-Output-1，單擊Edit.；在Edit History Output Request對話框中，設置輸出整個模型的內能和應變能，即ALLIE和ALLSE，單擊OK完成。5、 創建邊界條件及施加載荷l 創建邊界條件在環境欄Module后面選擇Load，進入Load模塊。單擊工具區的（Create

27、Boundary Condition），在Create Boundary Condition對話框中，Name后面輸入BC-X，Step選擇Initial，Category選擇Mechanical，Types for Selected Step選擇Symmetry.，單擊Continue.；在視圖區選擇裝配實例Lam-C-1左側端面，即垂直于X軸且X=0的側面，單擊鼠標中鍵或提示區的Done，在Edit Boundary Condition對話框中選擇XSYMM，單擊OK完成。類似操作，選擇Lam-C-1的下側端面，即垂直于Y軸且Y=0的側面，創建邊界條件BC-Y，邊界類型為YSYMM；選擇L

28、am-C-1垂直于Z軸且Z=0的側面，創建邊界條件BC-Z，邊界類型為ZSYMM。邊界條件及加載情況見圖&-1。l 施加載荷單擊工具區的(Create Load)，在Create Load對話框中，Name使用默認的Load-1，Step選擇Step-1，Category選擇Mechanical，Types for Selected Step選擇Pressure，單擊Continue.；在視圖區選擇Lam-C-1的右側端面，即垂直于X軸且X=15的側面，單擊鼠標中鍵，在Edit Load對話框中Magnitude后面輸入-100，單擊OK完成。6、 創建分析作業并提交分析l 創建分析作業在環境

29、欄Module后面選擇Job，進入Job模塊。單擊工具區的（Job Manager），在Job Manager對話框中單擊Create.；在Create Job對話框中，Name后面輸入Job-Lam-Stress-Umat，Source選擇Model-1，單擊Continue.；在Edit Job對話框的General標簽頁中，單擊User subroutine file后面的(Select.)，在相應路徑下找到并選擇UMAT-Stress.for文件；單擊Edit Job對話框中的OK完成。在Edit Job對話框的Parallelization標簽頁中可以設置多核并行計算。l 提交分析在

30、Job Manager對話框中，選中Job-Lam-Stress-Umat分析作業，單擊Submit提交計算。當Job-Lam-Stress-Umat的狀態（Status）由Running變為Completed時，計算完成，單擊Results進入可視化后處理模塊。l 保存模型單擊工具欄的File工具條中的（Save Model Database），在Save Model Database As對話框的File Name后面輸入Laminate-Umat，單擊OK完成。7、 修改材料在環境欄Module后面選擇Property，進入Property模塊。單擊工具區的(Composite Layu

31、p Manager)，在打開的對話框中選擇CompositeLayup-1，單擊Edit.；在Edit Composite Layup對話框的Plies標簽頁中，右單擊Material，單擊快捷菜單中的Edit Material.；在Select Material對話框中選擇Mat-T700，單擊OK；在Edit Composite Layup對話框中單擊OK完成。8、 再次創建分析作業并提交分析l 創建分析作業在環境欄Module后面選擇Job，進入Job模塊。單擊工具區的（Job Manager），在Job Manager對話框中單擊Create.；在Create Job對話框中，Name

32、后面輸入Job-Lam-Stress，Source選擇Model-1，單擊Continue.；在Edit Job對話框使用默認設置，單擊OK完成。在Edit Job對話框的Parallelization標簽頁中可以設置多核并行計算。l 提交分析在Job Manager對話框中，選中Job-Lam-Stress分析作業，單擊Submit提交計算。當Job-Lam-Stress的狀態（Status）由Running變為Completed時，計算完成，單擊Results進入可視化后處理模塊。l 保存模型單擊工具欄的File工具條中的(Save Model Database)保存模型。9、 可視化后處

33、理l 顯示云圖在視圖區顯示Job-Lam-Stress-Umat.odb。長按工具區的(Plot Contours on Deformed Shape)，顯示隱藏工具后單擊(Plot Contours on Undeformed Shape)，在Field Output工具條中設置輸出S11。單擊菜單欄ResultSection Points.，打開Section Points對話框，Selection method選擇Plies，在Plies區域選擇PLY-1，單擊Apply，在視圖區顯示該鋪層的S11應力云圖。相同操作，可以顯示各鋪層的各應力分離的云圖。Job-Lam-Stress-Uma

34、t.odb和Job-Lam-Stress.odb的各應力分量的云圖如圖&-4所示。(a) S11，Job-Lam-Stress-Umat(b) S11，Job-Lam-Stress(c) S22，Job-Lam-Stress-Umat(d) S22，Job-Lam-Stress(e) S12，Job-Lam-Stress-Umat(f) S12，Job-Lam-Stress圖&-4 各應力分量對比讀者在閱讀Abaqus幫助文件Abaqus Example Problems Guide中1.4.6 Failure of blunt notched fiber metal laminates一例的

35、exa_fml_ortho_damage_umat.f文件時注意SSE的計算公式。l 輸出應變能在視圖區顯示Job-Lam-Stress-Umat.odb。單擊工具區的(XY Data Manager)，在打開的對話框中單擊Create.；在Create XY Data對話框中選擇ODB history output，單擊Continue.；在History Output對話框中選擇ALLSE，單擊Save As.；在Save XY Data As對話框中Name使用默認的XYData-1，單擊OK；在XY Data Manager對話框中選擇XYData-1，單擊Edit.；在Edit XY

36、 Data對話框中可以看到應變能數據，分析結束時整個模型的應變能為9.671mJ。同樣操作，觀察Job-Lam-Stress.odb中整個模型的應變能數據，同樣為9.673mJ。通過對比，驗證了UMAT-Stress.for中第72行SSE計算的準確性。單擊菜單欄ResultHistory Output.，可以直接打開History Output對話框。&.2.4 Inp文件解釋本例中Job-Lam-Stress-Umat.inp節選如下。*Heading* Job name: Job-Lam-Stress-Umat Model name: Model-1* Generated by: Aba

37、qus/CAE 6.14-1*Preprint, echo=NO, model=NO, history=NO, contact=NO* PARTS* 部件Lam-C-1的節點、單元數據*Part, name=Lam-C*Node 1, 15., 10., 0.15 節點編號及坐標*Element, type=C3D8R 1, 64, 65, 23, 22, 43, 44, 2, 1 單元編號及節點編號*Orientation, name=Ori-1 1., 0., 0., 0., 1., 0.3, 0.* Section: CompositeLayup-1-1* 定義復合材料鋪層*Solid

38、Section, elset=CompositeLayup-1-1, composite, orientation=Ori-1, controls=EC-1, stack direction=3, layup=CompositeLayup-11, 1, UMat-T700, 0., Ply-1*End Part* ASSEMBLY*Assembly, name=Assembly* 使用部件Lam-C創建裝配實例Lam-C-1*Instance, name=Lam-C-1, part=Lam-C*End Instance*End Assembly* ELEMENT CONTROLS* 單元控制，

39、使用沙漏控制*Section Controls, name=EC-1, hourglass=ENHANCED1., 1., 1.* * MATERIALS* 使用工程常數創建材料Mat-T700*Material, name=Mat-T700*Elastic, type=ENGINEERING CONSTANTS114000.,8610.,8610., 0.3, 0.3, 0.45,4160.,4160.3000.,* 使用用戶材料創建材料UMat-T700*Material, name=UMat-T700*User Material, constants=6114000.,8610., 0.

40、3, 0.45,4160.,3000.* * BOUNDARY CONDITIONS* 創建對稱邊界條件BC-X，BC-Y和BC-Z略* Name: BC-X Type: Symmetry/Antisymmetry/Encastre*Boundary_PickedSet4, XSYMM* -* STEP: Step-1* 創建分析步Step-1*Step, name=Step-1, nlgeom=NO*Static1., 1., 1e-05, 1.* * LOADS* 在Step-1創建Pressure類型的載荷* Name: Load-1 Type: Pressure*Dsload_Pick

41、edSurf6, P, -100.* * OUTPUT REQUESTS* 設置輸出變量*Restart, write, frequency=0* * FIELD OUTPUT: F-Output-1* 場輸出變量U、E、S*Output, field, frequency=99999*Node OutputU, *Element Output, directions=YESE, S* * HISTORY OUTPUT: H-Output-1* 設置歷史輸出變量*Output, history*Energy OutputALLIE, ALLSE*End Step&.2.5 應用UMAT子程序應

42、力分析小結本例所介紹的應力分析是UMAT用戶子程序最簡單的應用，讀者在了解UMAT子程序的基礎上，可以結合相應的強度理論、剛度折減方法對UMAT子程序中的剛度矩陣、雅克比矩陣進行計算和更新，從而達到剛度折減和漸進損傷強度分析的目的。&.3 實例二：UMAT用戶子程序漸進損傷強度分析本例在實例一使用的UMAT子程序的基礎上，通過增加失效判定，對UMAT子程序中的剛度矩陣進行折減，達到漸進損傷強度分析的目的。為了使本例的UMAT子程序簡單易讀，本例對復合材料單層板進行纖維方向的拉伸強度分析，僅考慮纖維方向拉伸破壞。本例使用本章實例一的模型，稍做修改，用于本例的強度分析。&.3.1問題描述復合材料單

43、層板的3D實體模型的幾何尺寸及邊界條件同本章的實例一的模型。將右側面100MPa的拉力替換為X軸正方向0.4mm的位移載荷。本例中單位系統為mm、MPa。&.3.2 UMAT用戶子程序本例使用的UMAT用戶子程序UMAT-Strength.for的全部代碼如下，字母C及“!”之后為注釋內容。1 SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,2 1 RPL,DDSDDT,DRPLDE,DRPLDT,3 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,4 3 NDI,NSHR,NT

44、ENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,5 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)6 C 7 INCLUDE ABA_PARAM.INC8 C 9 CHARACTER*80 CMNAME10 DIMENSION STRESS(NTENS),STATEV(NSTATV),11 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),12 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRE

45、D(1),13 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),14 4 JSTEP(4)1516 DIMENSION STRESS0(6),STRAND(6),C(6,6),CD(6,6),DCDE(6,6)17 PARAMETER (ALPHA = 1000.0, LAMBDA = 0.0, DMAX = 0.9999, DRND = 3)18 C*19 C STRESS0.STRESS AT THE BEGINNING OF THE INCREMENT20 C STRAND.STRAIN AT THE END OF

46、THE INCREMENT21 C C.6X6 STIFFNESS MATRIX22 C CD.6X6 DAMAGED STIFFNESS MATRIX23 C DCDE.D CD/D E 24 C STATEV(1).DAMAGE VARIABLE D25 C*26 C GET THE MATERIAL PROPERTIES27 E1 = PROPS(1) !E1,YOUNGS MODULUS IN DIRECTION 128 E2 = PROPS(2) !E2 = E3,YOUNGS MODULUS IN DIRECTION 2 & 329 XNU12 = PROPS(3) !POISON

47、S RATIO POI_12,XNU13 = XNU1230 XNU21 = XNU12*E2/E1 !POISONS RATIO POI_21,XNU31 = XNU21 31 XNU23 = PROPS(4) !POISONS RATIO POI_23,XNU32 = XNU2332 G12 = PROPS(5) !G12 = G13,SHEAR MODULUS IN 12 & 13 PLANE33 G23 = PROPS(6) !G23,SHEAR MODULUS IN 23 PLANE34 STH = PROPS(7) !FAILURE STRESS IN 1 DIRECTION IN

48、 TENSION35 C*36 C STIFFNESS MATRIX C(6,6)37 RNU = 1/(1-2*XNU12*XNU21-XNU23*2-2*XNU12*XNU21*XNU23)38 C = 039 C(1,1) = E1*(1-XNU23*2)*RNU40 C(2,2) = E2*(1-XNU12*XNU21)*RNU41 C(3,3) = C(2,2)42 C(4,4) = G1243 C(5,5) = G12 44 C(6,6) = G2345 C(1,2) = E1*(XNU21+XNU21*XNU23)*RNU46 C(2,1) = C(1,2)47 C(1,3) =

49、 C(1,2)48 C(3,1) = C(1,2)49 C(2,3) = E2*(XNU23+XNU12*XNU21)*RNU50 C(3,2) = C(2,3)51 C CALCULATE THE STRAIN AT THE END OF THE INCREMENT52 DO I = 1, 653 STRAND(I) = STRAN(I) + DSTRAN(I)54 ENDDO55 C*56 C CALCULATE THE FAILURE COEFFICIENT 57 STRANF = STH/E158 IF (STRAND(1) 0) THEN59 F = STRAND(1)/STRANF60 ELSE61 F = 062 ENDIF63 C CALCULATE D,DAMAGE V