UMAT子程序在复合材料强度分析中的应用

-UMAT子程序在复合材料强度分析中的应用本例使用UMAT用户子程序进展复合材料单层板的应力分析和渐进损伤压缩强度分析，介绍UMAT用户子程序编写方法及在Abaqus/CAE中的设置。本章使用最大应变强度理论作为复合材料单层板的失效准则，相应的Fortran程序简单易读，便于理解UAMT子程序的工作原理。知识要点： 强度分析 UMAT用户子程序 最大应变理论 刚度折减讲师：孔祥宏版本：Abq 6.14难度：关键词：强度分析，UMAT&.1 本章容简介本章通过两个实例介绍UMAT用户子程序在复合材料单层板的应力分析和强度分析中的应用。在第一个实例中，对一个简单的复合材料单层板进展应力分析，UMAT子程序主要计算应力，不进展强度分析，本例用于验证UMAT子程序的计算精度。在第二个实例中，对复合材料单层板进展渐进损伤强度分析，UMAT子程序用于应力计算、强度分析和刚度折减。本章所用复合材料为T700/BA9916，材料属性如表&-1所示。表&-1 T700/BA9916材料属性参数值强度值E1/GPa114*T/MPa2688E2/GPa8.61*C/MPa1458E3/GPa8.61YT/MPa69.5120.3YC/MPa236130.3ZT/MPa55.5230.45ZC/MPa175G12/GPa4.16S*Y/MPa136G13/GPa4.16S*Z/MPa136G23/GPa3.0SYZ/MPa95.6&.2 实例一：UMAT用户子程序应力分析在使用UMAT用户子程序进展高级应用之前，应该先了解UMAT子程序，熟悉UMAT子程序的工作原理，了解UMAT中的参数、变量的含义。为了便于读者快速了解和使用UMAT，本例通过复合材料单层板的应力分析来介绍一个简单的UMAT子程序。读者可将本例中的单层板替换为层压板，进展比照分析。&.2.1问题描述复合材料单层板几何尺寸为15mm10mm0.15mm，纤维方向为45，单层板的3D实体模型如图&-1所示，*轴方向为0方向，左侧面施加*轴向对称边界条件，下侧面施加Y轴向对称边界条件，垂直于Z轴且Z=0的平面施加Z轴向对称边界条件，右侧面施加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,DROT,PNEWDT,5 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)6 C 7 INCLUDEABA_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(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),14 4 JSTEP(4)1516 DIMENSION EG(6), *NU(3,3), STRAND(6), C(6,6), STRESS0(6)17 C*18 C EG.E1,E2,E3,G12,G13,G2319 C *NU.NU12，NU21,NU13，NU31，NU23，NU3220 C STRAND.STRAINT AT THE END OF THE INCREMENT21 C C.6*6 STIFFNESS MATRI*22 C STRESS0.STRESS AT THE BEGINNING OF THE INCREMENT23 C*24 C INITIALIZE *NU & C MATRI*25 *NU=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 *NU(1,2) = PROPS(3) !POISONS RATIO POI_1232 *NU(2,1) = *NU(1,2)*EG(2)/EG(1) !POISONS RATIO POI_2133 *NU(1,3) = *NU(1,2) !POISONS RATIO POI_1334 *NU(3,1) = *NU(1,3)*EG(3)/EG(1) !POISONS RATIO POI_3135 *NU(2,3) = PROPS(4) !POISONS RATIO POI_2336 *NU(3,2) = *NU(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) = PROPS(6) !G23,SHEAR MODULUS IN 23 PLANE40 C*41 C FILL THE 6*6 STIFFNESS MATRI* C(6,6)42 RNU = 1/(1-*NU(1,2)*NU(2,1)-*NU(1,3)*NU(3,1)-43 1 *NU(3,2)*NU(2,3)-2*NU(1,3)*NU(2,1)*NU(3,2)44 C STIFFNESS MATRI* C(6,6)45 C(1,1) = EG(1)*(1-*NU(2,3)*NU(3,2)*RNU46 C(2,2) = EG(2)*(1-*NU(1,3)*NU(3,1)*RNU47 C(3,3) = EG(3)*(1-*NU(1,2)*NU(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)*(*NU(2,1)+*NU(3,1)*NU(2,3)*RNU52 C(2,1) = C(1,2)53 C(1,3) = EG(1)*(*NU(3,1)+*NU(2,1)*NU(3,2)*RNU54 C(3,1) = C(1,3)55 C(2,3) = EG(2)*(*NU(3,2)+*NU(1,2)*NU(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 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局部变量名及其含义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特定的弹性应变能、塑性耗散、蠕变耗散，只对能量输出有影响，对其他计算结果无影响，在UMAT中需要对其更新CELENT单元特征长度第15到83行为用户自己编写的固定格式的Fortran程序，用于计算刚度矩阵、应力、应变能、雅克比矩阵。由于本例中没有使用状态变量，因此不需要更新STATEV，只需要更新STRESS、DDSDDE和SSE即可。第16行定义了5个数组，其中EG、STRAND、STRESS0为一维数组，*NU、C为二维数组。第18到22行为注释局部，EG存放材料的3个弹性模量和3个剪切模量；STRAND存放当前增量步完毕时的应变(E11, E22, .)；STRESS0存放增量步开场时的应力(S11, S22, .)；*NU为33的二维矩阵，存放泊松比12、21、13、31、23、32；C为66的刚度矩阵。第25、26行初始化二维数组*NU、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)，式中上标表示增量步序号。 (&-3)第63到69行，将当前分析步开场时的应力STRESS的值赋给STRESS0，然后计算当前增量步的应力并赋给STRESS，当前增量步的应力计算如式(&-4)，式中上标表示增量步序号，应力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中的剪应变为工程剪应变。&.2.3复合材料单层板应力分析1、 创立部件及划分网格l 创立部件在Part模块，单击工具区的(Create Part)，在Create Part对话框中，Name后面输入Lam-C，Modeling Space选择3D，Type选择Deformable，在Base Feature区域选择Solid、E*trusion，Appro*imate size使用默认的200，单击Continue.进入绘图模式。单击工具区的(Create Lines: Rectangle (4 Lines)，在提示区输入第1个点的坐标(0,0)后按回车键，再输入第2个点的坐标(15,10)后按回车，再按Esc键或单击鼠标中键。单击提示区的Done或鼠标中键，在Edit Base E*trusion对话框Depth后面输入0.15，单击OK完成。l 划分网格在环境栏Module后面选择Mesh，进入Mesh模块。环境栏中Object选择Part: Lam-C。单击工具区的(Seed Part)，在Global Seeds对话框中Appro*imate global size后面输入0.5，单击OK。单击工具区的(Mesh Part)，单击提示区的Yes或鼠标中键，完成网格划分，如图&-2所示，板的厚度方向只划分为1层单元。图&-2 划分网格单击工具区的(Assign Element Type)，在Element Type对话框中，选择依次选择Standard，Linear，3D Stress，在He*标签页中勾选Reduced integration，在Element Controls区域Hourglass control选择Enhanced，即选择C3D8R单元，单击OK完成。单击工具区的(Assign Stack Direction)，在视图区选择部件平行于*-Y平面的面，单击鼠标中键或提示区的Yes完成。2、 创立材料并给部件赋材料属性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；单击MechanicalElasticityElastic，在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 posite Layup)，在翻开的对话框中Name使用默认名称，Initial ply count后面输入1，Element Type选择Solid，单击Continue.；在Edit posite Layup对话框中，Layup Orientation区域的Definition选择Coordinate 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 posite Layup对话框中单击(Select CSYS.)，在视图区选择刚创立的Datum csys-1，Stacking Direction选择Element direction 3，Rotation a*is选择A*is 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完成。在定义复合材料铺层时，视图区部件上会显示铺层方向，在Edit posite Layup对话框中的Display标签页中可以设置所需显示的方向，在视图区部件上白色箭头及字母S表示Stacking Direction。3、 装配在环境栏Module后面选择Assembly，进入Assembly模块。单击工具区的(Create Instance)，在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.；在Edit 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 Boundary Condition，在Create Boundary Condition对话框中，Name后面输入BC-*，Step选择Initial，Category选择Mechanical，Types for Selected Step选择Symmetry.，单击Continue.；在视图区选择装配实例Lam-C-1左侧端面，即垂直于*轴且*=0的侧面，单击鼠标中键或提示区的Done，在Edit Boundary Condition对话框中选择*SYMM，单击OK完成。类似操作，选择Lam-C-1的下侧端面，即垂直于Y轴且Y=0的侧面，创立边界条件BC-Y，边界类型为YSYMM；选择Lam-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的右侧端面，即垂直于*轴且*=15的侧面，单击鼠标中键，在Edit Load对话框中Magnitude后面输入-100，单击OK完成。6、 创立分析作业并提交分析l 创立分析作业在环境栏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 提交分析在Job Manager对话框中，选中Job-Lam-Stress-Umat分析作业，单击Submit提交计算。当Job-Lam-Stress-Umat的状态Status由Running变为pleted时，计算完成，单击Results进入可视化后处理模块。l 保存模型单击工具栏的File工具条中的Save Model Database，在Save Model Database As对话框的File Name后面输入Laminate-Umat，单击OK完成。7、 修改材料在环境栏Module后面选择Property，进入Property模块。单击工具区的(posite Layup Manager)，在翻开的对话框中选择positeLayup-1，单击Edit.；在Edit posite Layup对话框的Plies标签页中，右单击Material，单击快捷菜单中的Edit Material.；在Select Material对话框中选择Mat-T700，单击OK；在Edit posite Layup对话框中单击OK完成。8、 再次创立分析作业并提交分析l 创立分析作业在环境栏Module后面选择Job，进入Job模块。单击工具区的Job Manager，在Job Manager对话框中单击Create.；在Create Job对话框中，Name后面输入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变为pleted时，计算完成，单击Results进入可视化后处理模块。l 保存模型单击工具栏的File工具条中的(Save Model Database)保存模型。9、 可视化后处理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-Umat.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 E*ample Problems Guide中1.4.6 Failure of blunt notched fiber metal laminates一例的e*a_fml_ortho_damage_umat.f文件时注意SSE的计算公式。l 输出应变能在视图区显示Job-Lam-Stress-Umat.odb。单击工具区的(*Y Data Manager)，在翻开的对话框中单击Create.；在Create *Y Data对话框中选择ODB history output，单击Continue.；在History Output对话框中选择ALLSE，单击Save As.；在Save *Y Data As对话框中Name使用默认的*YData-1，单击OK；在*Y Data Manager对话框中选择*YData-1，单击Edit.；在Edit *Y 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: Abaqus/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: positeLayup-1-1* 定义复合材料铺层*Solid Section, elset=positeLayup-1-1, posite, orientation=Ori-1, controls=EC-1, stack direction=3, layup=positeLayup-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* 单元控制，使用沙漏控制*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.3, 0.45,4160.,3000.* * BOUNDARY CONDITIONS* 创立对称边界条件BC-*，BC-Y和BC-Z略* Name: BC-* Type: Symmetry/Antisymmetry/Encastre*Boundary_PickedSet4, *SYMM* -* 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_PickedSurf6, 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子程序应力分析小结本例所介绍的应力分析是UMAT用户子程序最简单的应用，读者在了解UMAT子程序的根底上，可以结合相应的强度理论、刚度折减方法对UMAT子程序中的刚度矩阵、雅克比矩阵进展计算和更新，从而到达刚度折减和渐进损伤强度分析的目的。&.3 实例二：UMAT用户子程序渐进损伤强度分析本例在实例一使用的UMAT子程序的根底上，通过增加失效判定，对UMAT子程序中的刚度矩阵进展折减，到达渐进损伤强度分析的目的。为了使本例的UMAT子程序简单易读，本例对复合材料单层板进展纤维方向的拉伸强度分析，仅考虑纤维方向拉伸破坏。本例使用本章实例一的模型，稍做修改，用于本例的强度分析。&.3.1问题描述复合材料单层板的3D实体模型的几何尺寸及边界条件同本章的实例一的模型。将右侧面100MPa的拉力替换为*轴正方向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,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,5 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)6 C 7 INCLUDEABA_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(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, DMA* = 0.9999, DRND = 3)18 C*19 C STRESS0.STRESS AT THE BEGINNING OF THE INCREMENT20 C STRAND.STRAIN AT THE END OF THE INCREMENT21 C C.6*6 STIFFNESS MATRI*22 C CD.6*6 DAMAGED STIFFNESS MATRI*23 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 *NU12 = PROPS(3) !POISONS RATIO POI_12,*NU13 = *NU1230 *NU21 = *NU12*E2/E1 !POISONS RATIO POI_21,*NU31 = *NU21 31 *NU23 = PROPS(4) !POISONS RATIO POI_23,*NU32 = *NU2332 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 TENSION35 C*36 C STIFFNESS MATRI* C(6,6)37 RNU = 1/(1-2*NU12*NU21-*NU23*2-2*NU12*NU21*NU23)38 C = 039 C(1,1) = E1*(1-*NU23*2)*RNU40 C(2,2) = E2*(1-*NU12*NU21)*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*(*NU21+*NU21*NU23)*RNU46 C(2,1) = C(1,2)47 C(1,3) = C(1,2)48 C(3,1) = C(1,2)49 C(2,3) = E2*(*NU23+*NU12*NU21)*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 VARIABLE64 D = STATEV(1)65 DDDE = 066 IF (F1) THEN67 DV = 1 - E*P(ALPHA*(1-F)68 IF (DV D) THEN69 D = D*LAMBDA/(LAMBDA+1)+DV/(LAMBDA+1)70 DDDE = ALPHA*(1-DV)/STRANF/(1+LAMBDA)71 D = ANINT(D*10*DRND)/10*DRND72 DDDE = ANINT(DDDE*10*DRND)/10*DRND73 ENDIF74 IF (DDMA*) THEN75 D=DMA*76 ENDIF77 ENDIF78 STATEV(1) = D !UPDATE D79 C DAMAGED STIFFNESS MATRI* CD(6,6)80 CD = C 81 CD(1,1) = (1-D)*C(1,1)82 CD(1,2) = (1-D)*C(1,2)83 CD(1,3) = CD(1,2)84 CD(2,1) = CD(1,2)85 CD(3,1) = CD(1,2)86 CD(4,4) = (1-D)*C(4,4)87 CD(5,5) = (1-D)*C(5,5) 88 C* 89 C CALCULATE STRESS 90 DO I = 1, 691 STRESS0(I) = STRESS(I)92 STRESS(I) = 0.093 DO J = 1, 694 STRESS(I) = STRESS(I)+CD(I,J)*STRAND(J)95 ENDDO96 ENDDO97 C CALCULATE SSE,ELASTIC STRAIN ENERGY98 DO I = 1, NTENS99 SSE = SSE+0.5*(STRESS0(I)+STRESS(I)*DSTRAN(I)100 ENDDO101 C* 102 C UPDATE THE JACOBIAN103 DCDE=0 104 DCDE(1,1) = -DDDE*CD(1,1)105 DCDE(1,2) = -DDDE*CD(1,2)106 DCDE(1,3) = DCDE(1,2)107 DCDE(2,1) = DCDE(1,2)108 DCDE(3,1) = DCDE(1,2)109 DCDE(4,4) = -DDDE*CD(4,4)110 DCDE(5,5) = -DDDE*CD(5,5)111 DDSDDE = CD112 DO I = 1,6113 DO J = 1,6114 ATEMP = DCDE(I,J)*STRAND(J)115 ENDDO116 DDSDDE(I,1) = DDSDDE(I,1)+ATEMP117 ENDDO118 RETURN119 END第1到15行及第118、119行为UMAT子程序固定格式，第16到117行为本程序的主体局部。第16行定义了5个数组，其含义见第18到25行之间的注释局部。由于本例仅考虑纤维方向拉伸破坏，且使用最大应变强度理论，所以DCDE为折减后的刚度矩阵Cd对纤维方向正应变1的偏导矩阵，即矩阵Cd中所以元素分别对1求偏导。本例只使用一个状态变量STATEV(1)记录单元纤维方向拉伸破坏的破坏变量d。第17行定义了4个常数参数，其中ALPHA(即)、LAMBDA(即)用于计算破坏变量d，DMA*(即dma*)为破坏变量d的最大值，DRND(即dr)为破坏变量d保存小数的位数。第27到34行，读取材料参数。其中，第34行STH为材料纤维方向的拉伸强度。第37到50行，计算材料无损伤时的刚度矩阵C。第52到54行，计算当前增量步完毕时的应变 = (1, 2, 3, 12, 13, 23)，其中剪应变为工程剪应变。第57到62行，使用最大应变强度理论计算纤维方向拉伸失效系数f。纤维方向的拉伸失效应变0及失效系数f的计算公式如下： (&-7)式(&-7)中为材料纤维方向拉伸强度，E1为材料纤维方向弹性模量，1为当前分析步完毕时纤维方向的拉伸正应变。第64到78行，计算破坏变量D即d和破坏系数对应变的偏导数导数DDDE即。破坏变量计算公式如下： (&-8)式(&-8)中dv为当前增量步完毕时的破坏变量，di为前一增量步完毕时或当前增量步开场时的破坏变量，di+1为当前增量步完毕时的破坏变量。由于仅考虑纤维方向拉伸破坏，式(&-7)中失效系数f仅为1的函数，破坏变量d仅为f的函数，因此DDDE可以写为： (&-9)式(&-9)中的d对应式(&-8)中的di+1。如果d对为的函数，则为一个向量，即。第71、72行，按照参数DRND的值对D和DDDE保存指定的小数位数。由于使用状态变量STATEV(1)保存破坏变量，因此要用式(&-8)中的di+1对STATEV(1)进展更新。第80到87行，