UMAT子程序在复合材料强度分析中的应用进损伤压缩强度分析，介绍UMAT用户子程序编写方法及在Abaqus/CAE中的设置。
本章使用最大应变强度理论作为复合材料单层板的失效准则，相应的Fortran程序简单易读，便于理解UAMT知识要点：强度分析UMAT用户子程序最大应变理论刚度折减&.1 本章内容简介本章通过两个实例介绍UMAT用户子程序在复合材料单层板的应力分析和强度分析中的应用。
在第一个实例中，对一个简单的复合材料单层板进行应力分析，UMAT子程序主要计算应力，不进行强度分析，本例用于验证UMAT子程序的计算精度。
在第二个实例中，对复合材料单层板进行渐进损伤强度分析，UMAT子程序用于应力计算、强度分析和刚度折减。
本章所用复合材料为T700/BA9916，材料属性如表&-1所示。
表&-1 T700/BA9916材料属性参数值强度值E1/GPa114X T/MPa2688E2/GPa X C/MPa1458E3/GPa Y T/MPaμ12Y C/MPa236μ13Z T/MPaμ23Z C/MPa175G12/GPa S XY/MPa136G13/GPa S XZ/MPa136G23/GPa S YZ/MPa&.2 实例一：UMAT 用户子程序应力分析&.问题描述复合材料单层板几何尺寸为15mm ×10mm ×，纤维方向为45°，单层板的3D 实体模型如图&-1所示，X 轴方向为0°方向，左侧面施加X 轴向对称边界条件，下侧面施加Y 轴向对称边界条件，垂直于Z 轴且Z=0的平面施加Z 轴向对称边界条件，右侧面施加100MPa 的拉力。
图&-1 单层板边界条件及加载情况本例中单位系统为mm 、MPa 。
&. UMAT 用户子程序本例使用的UMAT 用户子程序的全部代码如下，字母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,DROT,PNEWDT,5 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)6 C7 INCLUDE '' 8 C9 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(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3), 14 4 JSTEP(4) 1516 DIMENSION EG(6), XNU(3,3), STRAND(6), C(6,6), STRESS0(6) 17C****************************在使用UMAT 用户子程序进行高级应用之前，应该先了解UMAT 子程序，熟悉UMAT 子程序的工作原理，了解UMAT 中的参数、变量的含义。
为了便于读者快速了解和使用UMAT ，本例通过复合材料单层板的应力分析来介绍一个简单的UMAT 子程序。
读者可将本例中的单层板替换为层压板，进行对比分析。
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 INCREMENT23C****************************24 C INITIALIZE XNU & C MATRIX25 XNU=026 C=027 C GET THE MATERIAL PROPERTIES---ENGINEERING CONSTANTS28 EG(1) = PROPS(1) !E1,YOUNG'S MODULUS IN DIRECTION 129 EG(2) = PROPS(2) !E2,YOUNG'S MODULUS IN DIRECTION 230 EG(3) = EG(2) !E3,YOUNG'S MODULUS IN DIRECTION 331 XNU(1,2) = PROPS(3) !POISON'S RATIO POI_1232 XNU(2,1) = XNU(1,2)*EG(2)/EG(1) !POISON'S RATIO POI_2133 XNU(1,3) = XNU(1,2) !POISON'S RATIO POI_1334 XNU(3,1) = XNU(1,3)*EG(3)/EG(1) !POISON'S RATIO POI_3135 XNU(2,3) = PROPS(4) !POISON'S RATIO POI_2336 XNU(3,2) = XNU(2,3)*EG(3)/EG(2) !POISON'S RATIO POI_3237 EG(4) = PROPS(5) !G12,SHEAR MODULUS IN 12 PLANE38 EG(5) = EG(4) !G13,SHEAR MODULUS IN 13 PLANE39 EG(6) = PROPS(6) !G23,SHEAR MODULUS IN 23 PLANE40C****************************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) = 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,3)55 C(2,3) = EG(2)*(XNU(3,2)+XNU(1,2)*XNU(3,1))*RNU56 C(3,2) = C(2,3)57C****************************58 C CALCULATE STRAIN59DO I = 1, 660 STRAND(I) = STRAN(I)+DSTRAN(I)61ENDDO62 C CALCULATE STRESS63DO I = 1, 664 STRESS0(I) = STRESS(I)65 STRESS(I) = 066DO J = 1, 667 STRESS(I) = STRESS(I)+C(I,J)*STRAND(J)68ENDDO69ENDDO70 C CALCULATE SSE71DO I = 1, 672 SSE = SSE+*(STRESS0(I)+STRESS(I))*DSTRAN(I)73ENDDO74C****************************75 C UPDATE DDSDDE76DO I = 1, 677DO J = 1, 678 DDSDDE(I,J) = C(I,J)79 ENDDO 80 ENDDO 81 RETURN 82 END第1到14行及第81、82行为UMAT 子程序固定格式，其中，第1到5行括号内的变量为UMAT 子程序中可以使用的变量，第10到14行定义各变量数组的维数和长度。
部分主要变量的含义如表&-2所示。
表&-2 UMAT 部分变量名及其含义第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, ...)；XNU为3×3的二维矩阵，存放泊松比ν12、ν21、ν13、ν31、ν23、ν32；C 为6×6的刚度矩阵。
第25、26行初始化二维数组XNU 、C ，使其每个元素都为0。
第28到39行，读取材料常数，计算泊松比。
泊松比的计算公式如下。
ji ji ij jji iij E E ννE νE ν⋅==， (&-1)第42到56行，计算刚度矩阵。
刚度矩阵的计算公式如下。
322113322331132112233213311221132123331123222323121333221311133213122233121112236613551244211233331132223223111211)()()()()()()1()1()1(νννννννννγγνννγνννγνννγνννγνννγνννγννγννγνν⋅⋅⋅-⋅-⋅-⋅-=⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎧===⋅⋅-⋅=⋅⋅-⋅=⋅⋅-⋅=⋅⋅-⋅=⋅⋅-⋅=⋅⋅-⋅====⋅⋅-⋅=⋅⋅-⋅=⋅⋅-⋅=CC C C C C E E C E E C E E C G C G C G C E C E C E C (&-2) 对于本例所用材料，由于E 2=E 3，G 12=G 13，ν12=ν13，所以ν21=ν31、ν23=ν32，可以将式(&-2)化简后在UMAT 中计算刚度矩阵。