1、UMAT子程序在复合材料强度分析中的应用UMAT子程序在复合材料强度分析中的应用本例使用UMAT用户子程序进行复合材料单层板的应力分析和渐进损伤压缩强度分析，介绍UMAT用户子程序编写方法及在Abaqus/CAE中的设置。本章使用最大应变强度理论作为复合材料单层板的失效准则，相应的Fortran程序简单易读，便于理解UAMT子程序的工作原理。知识要点： 强度分析 UMAT用户子程序 最大应变理论 刚度折减讲师：孔祥宏版本：Abq 难度：关键词：强度分析，UMAT&.1 本章内容简介本章通过两个实例介绍UMAT用户子程序在复合材料单层板的应力分析和强度分析中的应用。在第一个实例中，对一个简单的复
2、合材料单层板进行应力分析，UMAT子程序主要计算应力，不进行强度分析，本例用于验证UMAT子程序的计算精度。在第二个实例中，对复合材料单层板进行渐进损伤强度分析，UMAT子程序用于应力计算、强度分析和刚度折减。本章所用复合材料为T700/BA9916，材料属性如表&-1所示。表&-1 T700/BA9916材料属性参数值强度值E1/GPa114XT/MPa2688E2/GPaXC/MPa1458E3/GPaYT/MPa12YC/MPa23613ZT/MPa23ZC/MPa175G12/GPaSXY/MPa136G13/GPaSXZ/MPa136G23/GPaSYZ/MPa&.2 实例一：UMA
3、T用户子程序应力分析在使用UMAT用户子程序进行高级应用之前，应该先了解UMAT子程序，熟悉UMAT子程序的工作原理，了解UMAT中的参数、变量的含义。为了便于读者快速了解和使用UMAT，本例通过复合材料单层板的应力分析来介绍一个简单的UMAT子程序。读者可将本例中的单层板替换为层压板，进行对比分析。&.问题描述复合材料单层板几何尺寸为15mm10mm，纤维方向为45，单层板的3D实体模型如图&-1所示，X轴方向为0方向，左侧面施加X轴向对称边界条件，下侧面施加Y轴向对称边界条件，垂直于Z轴且Z=0的平面施加Z轴向对称边界条件，右侧面施加100MPa的拉力。图&-1 单层板边界条件及加载情况本
4、例中单位系统为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,
5、KSPT,JSTEP,KINC)6C 7 INCLUDE 8C 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(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),14 4 JSTEP(4)15 16 DIMENSION E
6、G(6), XNU(3,3), STRAND(6), C(6,6), STRESS0(6)17C*18C EG.E1,E2,E3,G12,G13,G2319C XNU.NU12，NU21,NU13，NU31，NU23，NU3220C STRAND.STRAINT AT THE END OF THE INCREMENT21C C.6X6 STIFFNESS MATRIX22C STRESS0.STRESS AT THE BEGINNING OF THE INCREMENT23C*24C INITIALIZE XNU & C MATRIX25 XNU=026 C=027C GET THE MATE
7、RIAL 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 RATIO POI_1232 XNU(2,1) = XNU(1,2)*EG(2)/EG(1) !POISONS RATIO POI_2133 XNU(
8、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 MODULUS IN 12 PLANE38 EG(5) = EG(4) !G13,SHEAR MODULUS IN 13 PLANE39 EG(6) =
9、 PROPS(6) !G23,SHEAR MODULUS IN 23 PLANE40C*41C 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)44C 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
10、(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*58C CA
11、LCULATE STRAIN59 DO I = 1, 6 60 STRAND(I) = STRAN(I)+DSTRAN(I)61 ENDDO62C 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 ENDDO70C CALCULATE SSE71 DO I = 1, 672 SSE = SSE+*(STRESS0(I)+STRESS(I)*DSTRAN(I)73 ENDDO
12、74C*75C 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部分变量名及其含义STRESS增量步开始时的应力(S11, S22, .)，用增量步结束时的应力计算结果对其更新STATEV(NSTATV)状态变量（状态变量个数），如果在材料中定义了
13、状态变量，则在UMAT中需要对其更新STRAN增量步开始时的应变(E11, E22, .)DSTRAN当前增量步的应变增量(E11, E22, .)NDI, NSHR, NTENS应力、应变的个数，NDI为正应力或正应变的个数，NSHR为剪应力或剪应变的个数，NTENS=NDI+NSHRPROPS, NPROPS材料参数、材料参数的个数DDSDDE雅克比矩阵，SSE, SPD, SCD特定的弹性应变能、塑性耗散、蠕变耗散，只对能量输出有影响，对其他计算结果无影响，在UMAT中需要对其更新CELENT单元特征长度第15到83行为用户自己编写的固定格式的Fortran程序，用于计算刚度矩阵、应力、
14、应变能、雅克比矩阵。由于本例中没有使用状态变量，因此不需要更新STATEV，只需要更新STRESS、DDSDDE和SSE即可。第16行定义了5个数组，其中EG、STRAND、STRESS0为一维数组，XNU、C为二维数组。第18到22行为注释部分，EG存放材料的3个弹性模量和3个剪切模量；STRAND存放当前增量步结束时的应变(E11, E22, .)；STRESS0存放增量步开始时的应力(S11, S22, .)；XNU为33的二维矩阵，存放泊松比12、21、13、31、23、32；C为66的刚度矩阵。第25、26行初始化二维数组XNU、C，使其每个元素都为0。第28到39行，读取材料常数，
15、计算泊松比。泊松比的计算公式如下。 (&-1)第42到56行，计算刚度矩阵。刚度矩阵的计算公式如下。 (&-2)对于本例所用材料，由于E2=E3，G12=G13，12=13，所以21=31、23=32，可以将式(&-2)化简后在UMAT中计算刚度矩阵。本例为保证的可读性，刚度矩阵的计算没有化简，直接按照式(&-2)编写。第59到61行，计算当前增量步的应变，计算公式如式(&-3)，式中上标表示增量步序号。 (&-3)第63到69行，将当前分析步开始时的应力STRESS的值赋给STRESS0，然后计算当前增量步的应力并赋给STRESS，当前增量步的应力计算如式(&-4)，式中上标表示增量步序号，
16、应力i和应变i为列向量。 (&-4)第71到73行，计算应变能，如式(&-5)所示，式中上标表示增量步序号；下标表示应力、应变增量的分量序号，其中下标为1、2、3表示正应力、正应变增量，4、5、6表示剪应力、剪应变增量。 (&-5)第76到80行为更新雅克比矩阵。由于没有对刚度矩阵没变化，应力与应变的关系如式(&-4)所示，所以应力增量与应变增量的关系如式(&-6)所示，所以雅克比矩阵就等于刚度矩阵。 (&-6)第76到80行可简写为一行，如下所示：DDSDDE(1:6,1:6) = C(1:6,1:6)，或DDSDDE = CUMAT中的剪应变为工程剪应变。&.复合材料单层板应力分析1、创建
17、部件及划分网格 创建部件在Part模块，单击工具区的(Create Part)，在Create Part对话框中，Name后面输入Lam-C，Modeling Space选择3D，Type选择Deformable，在Base Feature区域选择Solid、Extrusion，Approximate size使用默认的200，单击Continue.进入绘图模式。单击工具区的(Create Lines: Rectangle (4 Lines)，在提示区输入第1个点的坐标(0,0)后按回车键，再输入第2个点的坐标(15,10)后按回车，再按Esc键或单击鼠标中键。单击提示区的Done或鼠标中键，
18、在Edit Base Extrusion对话框Depth后面输入，单击OK完成。 划分网格在环境栏Module后面选择Mesh，进入Mesh模块。环境栏中Object选择Part: Lam-C。单击工具区的(Seed Part)，在Global Seeds对话框中Approximate global size后面输入，单击OK。单击工具区的(Mesh Part)，单击提示区的Yes或鼠标中键，完成网格划分，如图&-2所示，板的厚度方向只划分为1层单元。图&-2 划分网格单击工具区的(Assign Element Type)，在Element Type对话框中，选择依次选择Standard，Li
19、near，3D Stress，在Hex标签页中勾选Reduced integration，在Element Controls区域Hourglass control选择Enhanced，即选择C3D8R单元，单击OK完成。单击工具区的(Assign Stack Direction)，在视图区选择部件平行于X-Y平面的面，单击鼠标中键或提示区的Yes完成。2、创建材料并给部件赋材料属性 创建材料在环境栏Module后面选择Property，进入Property模块。单击工具区的(Create Material)，在Edit Material对话框中，Name后面输入UMat-T700；单击Gene
20、ralUser Material，在User Material区域中Data区域的Mechanical Constants一栏依次输入114000， 8610，4160，3000，单击OK完成。单击工具区的(Create Material)，在Edit Material对话框中，Name后面输入Mat-T700；单击MechanicalElasticityElastic，在Elastic区域中Type选择Engineering Constants，在Data区域输入从左到右依次输入114000，8610，8610，4160，4160，3000，单击OK完成。材料UMat-T700用于UMAT用
21、户子程序，Mat-T700用于做对比分析。输入数据时，每输入完一行后按回车（Enter）键，光标会自动移到下一行。也可以通过右键快捷菜单添加或删除一行。本例UMAT子程序较简单，不需要使用状态变量，因此在材料UMat-T700中没有定义Depvar。 给部件赋材料属性单击工具区的(Create Composite Layup)，在打开的对话框中Name使用默认名称，Initial ply count后面输入1，Element Type选择Solid，单击Continue.；在Edit Composite Layup对话框中，Layup Orientation区域的Definition选择Coo
22、rdinate system，单击Definition下一行的(Create Datum CSYS)。在Create Datum CSYS对话框中使用默认名称Datum csys-1，类型选择Rectangular，单击Continue.，在提示区输入原点坐标(0,0,0)后按回车键，再输入(1,0,0)后按回车键，最后输入(0,1,0)后按回车键，单击Create Datum CSYS对话框的Cancel。在Edit Composite Layup对话框中单击(Select CSYS.)，在视图区选择刚创建的Datum csys-1，Stacking Direction选择Element d
23、irection 3，Rotation axis选择Axis 3。在Plies标签页中，双击Region，在视图区选择部件后单击鼠标中键；右单击Material，在快捷菜单中单击Edit Material.，在Select Material对话框中选择UMat-T700，单击OK；右单击Element Relative Thickness，在快捷菜单中单击Edit Thickness.，在Thickness对话框中Specify Value后面输入1后单击OK；在Rotation Angle一栏输入0；Integration Points使用默认的1，单击OK完成。在定义复合材料铺层时，视图区
24、部件上会显示铺层方向，在Edit Composite Layup对话框中的Display标签页中可以设置所需显示的方向，在视图区部件上白色箭头及字母S表示Stacking Direction。3、装配在环境栏Module后面选择Assembly，进入Assembly模块。单击工具区的(Create Instance)，在Create Instance对话框中选择Parts: Lam-C，单击OK完成。4、创建分析步、设置输出变量 创建分析步在环境栏Module后面选择Step，进入Step模块。单击工具区的(Create Step)，在Create Step对话框中，在Initial分析步之后
25、插入Static, General分析步，单击Continue.；在Edit Step对话框中使用默认设置，单击OK完成。 设置输出变量单击工具区的(Field Output Manager)，在Field Output Requests Manager对话框中，选中F-Output-1，单击Edit.；在Edit Field Output Request对话框中，设置如图&-3所示，输出整个模型最后一个增量步的S、E、U，单击OK完成。图&-3 场输出变量设置单击工具区的(History Output Manager)，在History Output Requests Manager对话框中
26、，选中H-Output-1，单击Edit.；在Edit History Output Request对话框中，设置输出整个模型的内能和应变能，即ALLIE和ALLSE，单击OK完成。5、创建边界条件及施加载荷 创建边界条件在环境栏Module后面选择Load，进入Load模块。单击工具区的（Create Boundary Condition），在Create Boundary Condition对话框中，Name后面输入BC-X，Step选择Initial，Category选择Mechanical，Types for Selected Step选择Symmetry.，单击Continue.；在
27、视图区选择装配实例Lam-C-1左侧端面，即垂直于X轴且X=0的侧面，单击鼠标中键或提示区的Done，在Edit Boundary Condition对话框中选择XSYMM，单击OK完成。类似操作，选择Lam-C-1的下侧端面，即垂直于Y轴且Y=0的侧面，创建边界条件BC-Y，边界类型为YSYMM；选择Lam-C-1垂直于Z轴且Z=0的侧面，创建边界条件BC-Z，边界类型为ZSYMM。边界条件及加载情况见图&-1。 施加载荷单击工具区的(Create Load)，在Create Load对话框中，Name使用默认的Load-1，Step选择Step-1，Category选择Mechanical
28、，Types for Selected Step选择Pressure，单击Continue.；在视图区选择Lam-C-1的右侧端面，即垂直于X轴且X=15的侧面，单击鼠标中键，在Edit Load对话框中Magnitude后面输入-100，单击OK完成。6、创建分析作业并提交分析 创建分析作业在环境栏Module后面选择Job，进入Job模块。单击工具区的（Job Manager），在Job Manager对话框中单击Create.；在Create Job对话框中，Name后面输入Job-Lam-Stress-Umat，Source选择Model-1，单击Continue.；在Edit Job
29、对话框的General标签页中，单击User subroutine file后面的(Select.)，在相应路径下找到并选择文件；单击Edit Job对话框中的OK完成。在Edit Job对话框的Parallelization标签页中可以设置多核并行计算。 提交分析在Job Manager对话框中，选中Job-Lam-Stress-Umat分析作业，单击Submit提交计算。当Job-Lam-Stress-Umat的状态（Status）由Running变为Completed时，计算完成，单击Results进入可视化后处理模块。 保存模型单击工具栏的File工具条中的（Save Model Da
30、tabase），在Save Model Database As对话框的File Name后面输入Laminate-Umat，单击OK完成。7、修改材料在环境栏Module后面选择Property，进入Property模块。单击工具区的(Composite Layup Manager)，在打开的对话框中选择CompositeLayup-1，单击Edit.；在Edit Composite Layup对话框的Plies标签页中，右单击Material，单击快捷菜单中的Edit Material.；在Select Material对话框中选择Mat-T700，单击OK；在Edit Composite Layup对话框中单击OK完成。8、再次创建分析作业并提交分析 创建分析作业在环境栏Module后面选择Job，进入Job模块。单击工具区的（Job Manager），在Job Mana
