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

UMAT子程序在复合材料强度分析中的应用UMAT子程序在复合材料强度分析中的应用UMAT子程序在复合材料强度分析中的应用UMAT子程序在复合材料强度分析中的应用编制仅供参考审核批准生效日期地址：电话：传真:邮编:UMAT子程序在复合材料强度分析中的应用本例使用UMAT用户子程序进行复合材料单层板的应力分析和渐进损伤压缩强度分析，介绍UMAT用户子程序编写方法及在Abaqus/CAE中的设置。本章使用最大应变强度理论作为复合材料单层板的失效准则，相应的Fortran程序简单易读，便于理解UAMT子程序的工作原理。知识要点：强度分析UMAT用户子程序最大应变理论刚度折减讲师：孔祥宏版本：Abq难度：关键词：强度分析，UMAT

本章内容简介本章通过两个实例介绍UMAT用户子程序在复合材料单层板的应力分析和强度分析中的应用。在第一个实例中，对一个简单的复合材料单层板进行应力分析，UMAT子程序主要计算应力，不进行强度分析，本例用于验证UMAT子程序的计算精度。在第二个实例中，对复合材料单层板进行渐进损伤强度分析，UMAT子程序用于应力计算、强度分析和刚度折减。本章所用复合材料为T700/BA9916，材料属性如表&-1所示。表&-1T700/BA9916材料属性参数值强度值E1/GPa114XT/MPa2688E2/GPaXC/MPa1458E3/GPaYT/MPaμ12YC/MPa236μ13ZT/MPaμ23ZC/MPa175G12/GPaSXY/MPa136G13/GPaSXZ/MPa136G23/GPaSYZ/MPa实例一：UMAT用户子程序应力分析在使用UMAT用户子程序进行高级应用之前，应该先了解UMAT子程序，熟悉UMAT子程序的工作原理，了解UMAT中的参数、变量的含义。为了便于读者快速了解和使用UMAT，本例通过复合材料单层板的应力分析来介绍一个简单的UMAT子程序。读者可将本例中的单层板替换为层压板，进行对比分析。&.问题描述复合材料单层板几何尺寸为15mm×10mm×，纤维方向为45°，单层板的3D实体模型如图&-1所示，X轴方向为0°方向，左侧面施加X轴向对称边界条件，下侧面施加Y轴向对称边界条件，垂直于Z轴且Z=0的平面施加Z轴向对称边界条件，右侧面施加100MPa的拉力。图&-1单层板边界条件及加载情况本例中单位系统为mm、MPa。&.UMAT用户子程序本例使用的UMAT用户子程序的全部代码如下，字母C及“!”之后为注释内容。SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1RPL,DDSDDT,DRPLDE,DRPLDT,2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)CINCLUDE''CCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),4JSTEP(4)DIMENSIONEG(6),XNU(3,3),STRAND(6),C(6,6),STRESS0(6)C****************************CEG.....E1,E2,E3,G12,G13,G23CXNU.....NU12，NU21,NU13，NU31，NU23，NU32CSTRAND.....STRAINTATTHEENDOFTHEINCREMENTCC.....6X6STIFFNESSMATRIXCSTRESS0.....STRESSATTHEBEGINNINGOFTHEINCREMENTC****************************CINITIALIZEXNU&CMATRIXXNU=0C=0CGETTHEMATERIALPROPERTIES---ENGINEERINGCONSTANTSEG(1)=PROPS(1)!E1,YOUNG'SMODULUSINDIRECTION1EG(2)=PROPS(2)!E2,YOUNG'SMODULUSINDIRECTION2EG(3)=EG(2)!E3,YOUNG'SMODULUSINDIRECTION3XNU(1,2)=PROPS(3)!POISON'SRATIOPOI_12XNU(2,1)=XNU(1,2)*EG(2)/EG(1)!POISON'SRATIOPOI_21XNU(1,3)=XNU(1,2)!POISON'SRATIOPOI_13XNU(3,1)=XNU(1,3)*EG(3)/EG(1)!POISON'SRATIOPOI_31XNU(2,3)=PROPS(4)!POISON'SRATIOPOI_23XNU(3,2)=XNU(2,3)*EG(3)/EG(2)!POISON'SRATIOPOI_32EG(4)=PROPS(5)!G12,SHEARMODULUSIN12PLANEEG(5)=EG(4)!G13,SHEARMODULUSIN13PLANEEG(6)=PROPS(6)!G23,SHEARMODULUSIN23PLANEC****************************CFILLTHE6X6STIFFNESSMATRIXC(6,6)RNU=1/(1-XNU(1,2)*XNU(2,1)-XNU(1,3)*XNU(3,1)-1XNU(3,2)*XNU(2,3)-2*XNU(1,3)*XNU(2,1)*XNU(3,2))CSTIFFNESSMATRIXC(6,6)C(1,1)=EG(1)*(1-XNU(2,3)*XNU(3,2))*RNUC(2,2)=EG(2)*(1-XNU(1,3)*XNU(3,1))*RNUC(3,3)=EG(3)*(1-XNU(1,2)*XNU(2,1))*RNUC(4,4)=EG(4)C(5,5)=EG(5)C(6,6)=EG(6)C(1,2)=EG(1)*(XNU(2,1)+XNU(3,1)*XNU(2,3))*RNUC(2,1)=C(1,2)C(1,3)=EG(1)*(XNU(3,1)+XNU(2,1)*XNU(3,2))*RNUC(3,1)=C(1,3)C(2,3)=EG(2)*(XNU(3,2)+XNU(1,2)*XNU(3,1))*RNUC(3,2)=C(2,3)C****************************CCALCULATESTRAINDOI=1,6STRAND(I)=STRAN(I)+DSTRAN(I)ENDDOCCALCULATESTRESSDOI=1,6STRESS0(I)=STRESS(I)STRESS(I)=0DOJ=1,6STRESS(I)=STRESS(I)+C(I,J)*STRAND(J)ENDDOENDDOCCALCULATESSEDOI=1,6SSE=SSE+*(STRESS0(I)+STRESS(I))*DSTRAN(I)ENDDOC****************************CUPDATEDDSDDEDOI=1,6DOJ=1,6DDSDDE(I,J)=C(I,J)ENDDOENDDORETURNEND第1到14行及第81、82行为UMAT子程序固定格式，其中，第1到5行括号内的变量为UMAT子程序中可以使用的变量，第10到14行定义各变量数组的维数和长度。部分主要变量的含义如表&-2所示。表&-2UMAT部分变量名及其含义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为3×3的二维矩阵，存放泊松比ν12、ν21、ν13、ν31、ν23、ν32；C为6×6的刚度矩阵。第25、26行初始化二维数组XNU、C，使其每个元素都为0。第28到39行，读取材料常数，计算泊松比。泊松比的计算公式如下。(&-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)，式中上标表示增量步序号，应力σ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中的剪应变为工程剪应变。&.复合材料单层板应力分析创建部件及划分网格创建部件在Part模块，单击工具区的(CreatePart)，在CreatePart对话框中，Name后面输入Lam-C，ModelingSpace选择3D，Type选择Deformable，在BaseFeature区域选择Solid、Extrusion，Approximatesize使用默认的200，单击Continue...进入绘图模式。单击工具区的(CreateLines:Rectangle(4Lines))，在提示区输入第1个点的坐标(0,0)后按回车键，再输入第2个点的坐标(15,10)后按回车，再按Esc键或单击鼠标中键。单击提示区的Done或鼠标中键，在EditBaseExtrusion对话框Depth后面输入，单击OK完成。划分网格在环境栏Module后面选择Mesh，进入Mesh模块。环境栏中Object选择Part:Lam-C。单击工具区的(SeedPart)，在GlobalSeeds对话框中Approximateglobalsize后面输入，单击OK。单击工具区的(MeshPart)，单击提示区的Yes或鼠标中键，完成网格划分，如图&-2所示，板的厚度方向只划分为1层单元。图&-2划分网格单击工具区的(AssignElementType)，在ElementType对话框中，选择依次选择Standard，Linear，3DStress，在Hex标签页中勾选Reducedintegration，在ElementControls区域Hourglasscontrol选择Enhanced，即选择C3D8R单元，单击OK完成。单击工具区的(AssignStackDirection)，在视图区选择部件平行于X-Y平面的面，单击鼠标中键或提示区的Yes完成。创建材料并给部件赋材料属性创建材料在环境栏Module后面选择Property，进入Property模块。单击工具区的(CreateMaterial)，在EditMaterial对话框中，Name后面输入UMat-T700；单击General→UserMaterial，在UserMaterial区域中Data区域的MechanicalConstants一栏依次输入114000，8610，，，4160，3000，单击OK完成。单击工具区的(CreateMaterial)，在EditMaterial对话框中，Name后面输入Mat-T700；单击Mechanical→Elasticity→Elastic，在Elastic区域中Type选择EngineeringConstants，在Data区域输入从左到右依次输入114000，8610，8610，，，，4160，4160，3000，单击OK完成。材料UMat-T700用于UMAT用户子程序，Mat-T700用于做对比分析。输入数据时，每输入完一行后按回车（Enter）键，光标会自动移到下一行。也可以通过右键快捷菜单添加或删除一行。本例UMAT子程序较简单，不需要使用状态变量，因此在材料UMat-T700中没有定义Depvar。给部件赋材料属性单击工具区的(CreateCompositeLayup)，在打开的对话框中Name使用默认名称，Initialplycount后面输入1，ElementType选择Solid，单击Continue...；在EditCompositeLayup对话框中，LayupOrientation区域的Definition选择Coordinatesystem，单击Definition下一行的(CreateDatumCSYS)。在CreateDatumCSYS对话框中使用默认名称Datumcsys-1，类型选择Rectangular，单击Continue...，在提示区输入原点坐标(0,0,0)后按回车键，再输入(1,0,0)后按回车键，最后输入(0,1,0)后按回车键，单击CreateDatumCSYS对话框的Cancel。在EditCompositeLayup对话框中单击(SelectCSYS...)，在视图区选择刚创建的Datumcsys-1，StackingDirection选择Elementdirection3，Rotationaxis选择Axis3。在Plies标签页中，双击Region，在视图区选择部件后单击鼠标中键；右单击Material，在快捷菜单中单击EditMaterial...，在SelectMaterial对话框中选择UMat-T700，单击OK；右单击ElementRelativeThickness，在快捷菜单中单击EditThickness...，在Thickness对话框中SpecifyValue后面输入1后单击OK；在RotationAngle一栏输入0；IntegrationPoints使用默认的1，单击OK完成。在定义复合材料铺层时，视图区部件上会显示铺层方向，在EditCompositeLayup对话框中的Display标签页中可以设置所需显示的方向，在视图区部件上白色箭头及字母S表示StackingDirection。装配在环境栏Module后面选择Assembly，进入Assembly模块。单击工具区的(CreateInstance)，在CreateInstance对话框中选择Parts:Lam-C，单击OK完成。创建分析步、设置输出变量创建分析步在环境栏Module后面选择Step，进入Step模块。单击工具区的(CreateStep)，在CreateStep对话框中，在Initial分析步之后插入Static,General分析步，单击Continue...；在EditStep对话框中使用默认设置，单击OK完成。设置输出变量单击工具区的(FieldOutputManager)，在FieldOutputRequestsManager对话框中，选中F-Output-1，单击Edit...；在EditFieldOutputRequest对话框中，设置如图&-3所示，输出整个模型最后一个增量步的S、E、U，单击OK完成。图&-3场输出变量设置单击工具区的(HistoryOutputManager)，在HistoryOutputRequestsManager对话框中，选中H-Output-1，单击Edit...；在EditHistoryOutputRequest对话框中，设置输出整个模型的内能和应变能，即ALLIE和ALLSE，单击OK完成。创建边界条件及施加载荷创建边界条件在环境栏Module后面选择Load，进入Load模块。单击工具区的（CreateBoundaryCondition），在CreateBoundaryCondition对话框中，Name后面输入BC-X，Step选择Initial，Category选择Mechanical，TypesforSelectedStep选择Symmetry...，单击Continue...；在视图区选择装配实例Lam-C-1左侧端面，即垂直于X轴且X=0的侧面，单击鼠标中键或提示区的Done，在EditBoundaryCondition对话框中选择XSYMM，单击OK完成。类似操作，选择Lam-C-1的下侧端面，即垂直于Y轴且Y=0的侧面，创建边界条件BC-Y，边界类型为YSYMM；选择Lam-C-1垂直于Z轴且Z=0的侧面，创建边界条件BC-Z，边界类型为ZSYMM。边界条件及加载情况见图&-1。施加载荷单击工具区的(CreateLoad)，在CreateLoad对话框中，Name使用默认的Load-1，Step选择Step-1，Category选择Mechanical，TypesforSelectedStep选择Pressure，单击Continue...；在视图区选择Lam-C-1的右侧端面，即垂直于X轴且X=15的侧面，单击鼠标中键，在EditLoad对话框中Magnitude后面输入-100，单击OK完成。创建分析作业并提交分析创建分析作业在环境栏Module后面选择Job，进入Job模块。单击工具区的（JobManager），在JobManager对话框中单击Create...；在CreateJob对话框中，Name后面输入Job-Lam-Stress-Umat，Source选择Model-1，单击Continue...；在EditJob对话框的General标签页中，单击Usersubroutinefile后面的(Select...)，在相应路径下找到并选择文件；单击EditJob对话框中的OK完成。在EditJob对话框的Parallelization标签页中可以设置多核并行计算。提交分析在JobManager对话框中，选中Job-Lam-Stress-Umat分析作业，单击Submit提交计算。当Job-Lam-Stress-Umat的状态（Status）由Running变为Completed时，计算完成，单击Results进入可视化后处理模块。保存模型单击工具栏的File工具条中的（SaveModelDatabase），在SaveModelDatabaseAs对话框的FileName后面输入Laminate-Umat，单击OK完成。修改材料在环境栏Module后面选择Property，进入Property模块。单击工具区的(CompositeLayupManager)，在打开的对话框中选择CompositeLayup-1，单击Edit...；在EditCompositeLayup对话框的Plies标签页中，右单击Material，单击快捷菜单中的EditMaterial...；在SelectMaterial对话框中选择Mat-T700，单击OK；在EditCompositeLayup对话框中单击OK完成。再次创建分析作业并提交分析创建分析作业在环境栏Module后面选择Job，进入Job模块。单击工具区的（JobManager），在JobManager对话框中单击Create...；在CreateJob对话框中，Name后面输入Job-Lam-Stress，Source选择Model-1，单击Continue...；在EditJob对话框使用默认设置，单击OK完成。在EditJob对话框的Parallelization标签页中可以设置多核并行计算。提交分析在JobManager对话框中，选中Job-Lam-Stress分析作业，单击Submit提交计算。当Job-Lam-Stress的状态（Status）由Running变为Completed时，计算完成，单击Results进入可视化后处理模块。保存模型单击工具栏的File工具条中的(SaveModelDatabase)保存模型。可视化后处理显示云图在视图区显示。长按工具区的(PlotContoursonDeformedShape)，显示隐藏工具后单击(PlotContoursonUndeformedShape)，在FieldOutput工具条中设置输出S11。单击菜单栏Result→SectionPoints...，打开SectionPoints对话框，Selectionmethod选择Plies，在Plies区域选择PLY-1，单击Apply，在视图区显示该铺层的S11应力云图。相同操作，可以显示各铺层的各应力分离的云图。和的各应力分量的云图如图&-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帮助文件AbaqusExampleProblemsGuide中Failureofbluntnotchedfibermetallaminates一例的文件时注意SSE的计算公式。输出应变能在视图区显示。单击工具区的(XYDataManager)，在打开的对话框中单击Create...；在CreateXYData对话框中选择ODBhistoryoutput，单击Continue...；在HistoryOutput对话框中选择ALLSE，单击SaveAs...；在SaveXYDataAs对话框中Name使用默认的XYData-1，单击OK；在XYDataManager对话框中选择XYData-1，单击Edit...；在EditXYData对话框中可以看到应变能数据，分析结束时整个模型的应变能为。同样操作，观察中整个模型的应变能数据，同样为。通过对比，验证了中第72行SSE计算的准确性。单击菜单栏Result→HistoryOutput...，可以直接打开HistoryOutput对话框。&.Inp文件解释本例中节选如下。*Heading**Jobname:Job-Lam-Stress-UmatModelname:Model-1**Generatedby:Abaqus/CAE*Preprint,echo=NO,model=NO,history=NO,contact=NO****PARTS**部件Lam-C-1的节点、单元数据*Part,name=Lam-C*Node1,15.,10.,……节点编号及坐标*Element,type=C3D8R1,64,65,23,22,43,44,2,1……单元编号及节点编号*Orientation,name=Ori-11.,0.,0.,0.,1.,0.3,0.**Section:CompositeLayup-1-1**定义复合材料铺层*SolidSection,elset=CompositeLayup-1-1,composite,orientation=Ori-1,controls=EC-1,stackdirection=3,layup=CompositeLayup-11,1,UMat-T700,0.,Ply-1*EndPart****ASSEMBLY*Assembly,name=Assembly**使用部件Lam-C创建装配实例Lam-C-1*Instance,name=Lam-C-1,part=Lam-C*EndInstance*EndAssembly**ELEMENTCONTROLS**单元控制，使用沙漏控制*SectionControls,name=EC-1,hourglass=ENHANCED1.,1.,1.****MATERIALS**使用工程常数创建材料Mat-T700*Material,name=Mat-T700*Elastic,type=ENGINEERINGCONSTANTS114000.,8610.,8610.,,,,4160.,4160.3000.,**使用用户材料创建材料UMat-T700*Material,name=UMat-T700*UserMaterial,constants=6114000.,8610.,,,4160.,3000.****BOUNDARYCONDITIONS**创建对称边界条件BC-X，BC-Y和BC-Z略**Name:BC-XType: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-1Type:Pressure*Dsload_PickedSurf6,P,-100.****OUTPUTREQUESTS**设置输出变量*Restart,write,frequency=0****FIELDOUTPUT:F-Output-1**场输出变量U、E、S*Output,field,frequency=99999*NodeOutputU,*ElementOutput,directions=YESE,S****HISTORYOUTPUT:H-Output-1**设置历史输出变量*Output,history*EnergyOutputALLIE,ALLSE*EndStep&.应用UMAT子程序应力分析小结本例所介绍的应力分析是UMAT用户子程序最简单的应用，读者在了解UMAT子程序的基础上，可以结合相应的强度理论、刚度折减方法对UMAT子程序中的刚度矩阵、雅克比矩阵进行计算和更新，从而达到刚度折减和渐进损伤强度分析的目的。实例二：UMAT用户子程序渐进损伤强度分析本例在实例一使用的UMAT子程序的基础上，通过增加失效判定，对UMAT子程序中的刚度矩阵进行折减，达到渐进损伤强度分析的目的。为了使本例的UMAT子程序简单易读，本例对复合材料单层板进行纤维方向的拉伸强度分析，仅考虑纤维方向拉伸破坏。本例使用本章实例一的模型，稍做修改，用于本例的强度分析。&.问题描述复合材料单层板的3D实体模型的几何尺寸及边界条件同本章的实例一的模型。将右侧面100MPa的拉力替换为X轴正方向的位移载荷。本例中单位系统为mm、MPa。&.UMAT用户子程序本例使用的UMAT用户子程序的全部代码如下，字母C及“!”之后为注释内容。SUBROUTINEUMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1RPL,DDSDDT,DRPLDE,DRPLDT,2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)CINCLUDE''CCHARACTER*80CMNAMEDIMENSIONSTRESS(NTENS),STATEV(NSTATV),1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),4JSTEP(4)DIMENSIONSTRESS0(6),STRAND(6),C(6,6),CD(6,6),DCDE(6,6)PARAMETER(ALPHA=,LAMBDA=,DMAX=,DRND=3)C****************************CSTRESS0.....STRESSATTHEBEGINNINGOFTHEINCREMENTCSTRAND.....STRAINATTHEENDOFTHEINCREMENTCC.....6X6STIFFNESSMATRIXCCD.....6X6DAMAGEDSTIFFNESSMATRIXCDCDE.....DCD/DECSTATEV(1).....DAMAGEVARIABLEDC****************************CGETTHEMATERIALPROPERTIESE1=PROPS(1)!E1,YOUNG'SMODULUSINDIRECTION1E2=PROPS(2)!E2=E3,YOUNG'SMODULUSINDIRECTION2&3XNU12=PROPS(3)!POISON'SRATIOPOI_12,XNU13=XNU12XNU21=XNU12*E2/E1!POISON'SRATIOPOI_21,XNU31=XNU21XNU23=PROPS(4)!POISON'SRATIOPOI_23,XNU32=XNU23G12=PROPS(5)!G12=G13,SHEARMODULUSIN12&13PLANEG23=PROPS(6)!G23,SHEARMODULUSIN23PLANESTH=PROPS(7)!FAILURESTRESSIN1DIRECTIONINTENSIONC****************************CSTIFFNESSMATRIXC(6,6)RNU=1/(1-2*XNU12*XNU21-XNU23**2-2*XNU12*XNU21*XNU23)C=0C(1,1)=E1*(1-XNU23**2)*RNUC(2,2)=E2*(1-XNU12*XNU21)*RNUC(3,3)=C(2,2)C(4,4)=G12C(5,5)=G12C(6,6)=G23C(1,2)=E1*(XNU21+XNU21*XNU23)*RNUC(2,1)=C(1,2)C(1,3)=C(1,2)C(3,1)=C(1,2)C(2,3)=E2*(XNU23+XNU12*XNU21)*RNUC(3,2)=C(2,3)CCALCULATETHESTRAINATTHEENDOFTHEINCREMENTDOI=1,6STRAND(I)=STRAN(I)+DSTRAN(I)ENDDOC****************************CCALCULATETHEFAILURECOEFFICIENTSTRANF=STH/E1IF(STRAND(1)>0)THENF=STRAND(1)/STRANFELSEF=0ENDIFCCALCULATED,DAMAGEVARIABLED=STATEV(1)DDDE=0IF(F>1)THENDV=1-EXP(ALPHA*(1-F))IF(DV>D)THEND=D*LAMBDA/(LAMBDA+1)+DV/(LAMBDA+1)DDDE=ALPHA*(1-DV)/STRANF/(1+LAMBDA)D=ANINT(D*10**DRND)/10**DRNDDDDE=ANINT(DDDE*10**DRND)/10**DRNDENDIFIF(D>DMAX)THEND=DMAXENDIFENDIFSTATEV(1)=D!UPDATEDCDAMAGEDSTIFFNESSMATRIXCD(6,6)CD=CCD(1,1)=(1-D)*C(1,1)CD(1,2)=(1-D)*C(1,2)CD(1,3)=CD(1,2)CD(2,1)=CD(1,2)CD(3,1)=CD(1,2)CD(4,4)=(1-D)*C(4,4)CD(5,5)=(1-D)*C(5,5)C****************************CCALCULATESTRESSDOI=1,6STRESS0(I)=STRESS(I)STRESS(I)=DOJ=1,6STRESS(I)=STRESS(I)+CD(I,J)*STRAND(J)ENDDOENDDOCCALCULATESSE,ELASTICSTRAINENERGYDOI=1,NTENSSSE=SSE+*(STRESS0(I)+STRESS(I))*DSTRAN(I)ENDDOC****************************CUPDATETHEJACOBIANDCDE=0DCDE(1,1)=-DDDE*CD(1,1)DCDE(1,2)=-DDDE*CD(1,2)DCDE(1,3)=DCDE(1,2)DCDE(2,1)=DCDE(1,2)DCDE(3,1)=DCDE(1,2)DCDE(4,4)=-DDDE*CD(4,4)DCDE(5,5)=-DDDE*CD(5,5)DDSDDE=CDDOI=1,6DOJ=1,6ATEMP=DCDE(I,J)*STRAND(J)ENDDODDSDDE(I,1)=DDSDDE(I,1)+ATEMPENDDORETURNEND第1到15行及第118、119行为UMAT子程序固定格式，第16到117行为本程序的主体部分。第16行定义了5个数组，其含义见第18到25行之间的注释部分。由于本例仅考虑纤维方向拉伸破坏，且使用最大应变强度理论，所以DCDE为折减后的刚度矩阵Cd对纤维方向正应变ε1的偏导矩阵，即矩阵Cd中所以元素分别对ε1求偏导。本例只使用一个状态变量STATEV(1)记录单元纤维方向拉伸破坏的破坏变量d。第17行定义了4个常数参数，其中ALPHA(即α)、LAMBDA(即λ)用于计算破坏变量d，DMAX(即dmax)为破坏变量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行，计算破坏后的刚度矩阵Cd，由于近考虑纤维方向拉伸破坏，因此刚度折减如式(&-10)所示。(&-10)第90到100行，计算应力STRESS和应变能SSE。第103到110行，计算DCDE，即，对于本例，实际上只计算。第111到117行，计算雅克比矩阵。雅克比矩阵的基本计算公式如下：(&-11)式(&-11)中J为NTENS阶方阵，其中Jij可表示为：(&-12)在式(&-11)中，(&-13)式(&-13)中，(i=1,2,...)为一个列向量，该向量各元素分别加到雅克比矩阵J的第i列的相应元素上。&.复合材料单层板强度分析修改材料及复合材料铺层修改材料打开本章实例一创建的。在环境栏Module后面选择Property，进入Property模块。单击工具区的(MaterialManager)，在MaterialManager对话框中选择UMat-T700，单击Edit...；在EditMaterial对话框中单击General→Depvar，在Numberofsolution-dependentstatevariable后面输入1；单击MaterialBehaviors区域的UserMaterial，勾选Useunsymmetricmaterialstiffnessmatirx，在Data区域MechanicalConstants一栏的最后添加2688，即材料纤维方向的拉伸强度2688MPa；单击OK完成。修改复合材料铺层单击工具区的(CompositeLayupManager)，在打开的对话框中选择CompositeLayup-1，单击Edit...；在Plies标签页中，将Ply-1的Material改为UMat-T700，RotationAngle改为0，单击OK完成。创建参考点和约束创建参考点在环境栏Module后面选择Interaction，进入Interaction模块。单击工具区的(CreateReferencePoint)，在提示区输入(15,5,，按回车键完成，单击单击工具区的退出创建参考点。创建集合单击菜单栏Tools→Set→Manager...，在SetManager对话框中单击Create...；在CreateSet对话框中Name后面输入Set-RP，Type选择Geometry，单击Continue...；在视图区选择RP-1，单击鼠标中键完成。再单击SetManager对话框中单击Create...；在CreateSet对话框中Name后面输入Set-Nds，Type选择Node，单击Continue...；在视图区选择部件右侧面（垂直于X轴且X=15的侧面）上的所有节点，单击鼠标中键完成。选择一个面上的所有节点时，可以在提示区将选择方法设为byangle。创建约束单击工具区的(CreateConstraint)，在打开的对话框中Name使用默认名称，Type选择Equation，单击Continue...；在EditConstraint对话框设置如图&-5所示，单击OK完成。图&-5设置约束修改分析步、输出变量修改分析步在环境栏Module后面选择Step，进入Step模块。单击工具区的(StepManager)，在StepManager对话框中选择Step-1，单击Edit...；在EditStep对话框的Basic标签页中，Nlgeom选择On；在Incrementation标签页设置如图&-6所示，单击OK完成。图&-6增量步设置修改输出变量单击工具区的(FieldOutputManager)，在FieldOutputRequestsManager对话框中，选中F-Output-1，单击Edit...；在EditFieldOutputRequest对话框中，设置如图&-7所示，单击OK完成。图&-7场输出变量设置单击工具区的(HistoryOutputManager)，在HistoryOutputRequestsManager对话框中，选中H-Output-1，单击Edit...；在EditHistoryOutputRequest对话框中设置如图&-8，单击OK完成。图&-8历史输出变量设置修改载荷删除旧载荷在环境栏Module后面选择Load，进入Load模块。单击工具区的（LoadManager），在LoadManager对话框中选择Load-1，单击Delete...，单击提示框中的Yes，关闭LoadManager对话框。可以再模型树Loads中选择Load-1，按Delete键删除。施加位移载荷单击工具区的(CreateBoundaryCondition)，在打开的对话框中，Name后面输入U1，Step选择Step-1，Category选