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/GPa114XT/MPa2688E2/GPa8.61XC/MPa1458E3/GPa8.61YT/MPa69.5120.3YC/MPa236130.3ZT/MPa55.5230.45ZC/MPa175G12/GPa4.16SXY/MPa136G13/GPa4.16SXZ/MPa136G23/GPa3.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轴向对称边界条件，右侧面施加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 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(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 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 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 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) = 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)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为一维数组，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行，读取材料常数，计算泊松比。泊松比的计算公式如下。 (&-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、Extrusion，Approximate size使用默认的200，单击Continue.进入绘图模式。单击工具区的(Create 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)，单击提示区的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、 创建材料并给部件赋材料属性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 Composite Layup)，在打开的对话框中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)后按回车键，最后输入(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 Relative 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)，在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-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；选择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的右侧端面，即垂直于X轴且X=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变为Completed时，计算完成，单击Results进入可视化后处理模块。l 保存模型单击工具栏的File工具条中的（Save Model Database），在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、 再次创建分析作业并提交分析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变为Completed时，计算完成，单击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 Example Problems Guide中1.4.6 Failure of blunt notched fiber metal laminates一例的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 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: CompositeLayup-1-1* 定义复合材料铺层*Solid 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* 单元控制，使用沙漏控制*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-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_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的拉力替换为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,NTENS,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),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, DMAX = 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.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) !POISONS 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 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) = 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