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

1、讲师：孔祥宏版本： Abq 6.14难度：关键词：强度分析， UMATUMAT 子程序在复合材料强度分析中的应用本例使用UMAT 用户子程序进行复合材料单层板的应力分析和渐进损伤压缩强度分析，介绍UMA T 用户子程序编写方法及在Abaqus/CAE 中的设置。本章使用最大应变强度理论作为复合材料单层板的失效准则，相应的Fortran 程序简单易读，便于理解UAMT 子程序的工作原理。知识要点：强度分析UMAT 用户子程序最大应变理论刚度折减&.1 本章内容简介本章通过两个实例介绍UMAT 用户子程序在复合材料单层板的应力分析和强度分析中的应用。在第一个实例中，对一个简单的复合材料单层

2、板进行应力分析，UMAT 子程序主要计算应力，不进行强度分析，本例用于验证UMA T 子程序的计算精度。在第二个实例中，对复合材料单层板进行渐进损伤强度分析，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/MPa136G

3、13/GPa4.16Sxz/MPa136G23/GPa3.0SYZ/MPa95.6&.2 实例一：UMAT 用户子程序应力分析在使用 UMAT 用户子程序进行高级应用之前，应该先了解UMAT子程序，熟悉UMAT 子程序的工作原理，了解UMAT 中的参数、变量的含义。为了便于读者快速了解和使用UMA T，本例通过复合材料单层板的应力分析来介绍一个简单的UMAT 子程序。读者可将本例中的单层板替换为层压板，进行对比分析。&.2.1 问题描述复合材料单层板几何尺寸为15mm× 10mm× 0.15mm，纤维方向为 45° ，单层板的3D 实体模型如图&a

4、mp;-1 所示， X 轴方向为0°方向，左侧面施加X 轴向对称边界条件，下侧面施加Y 轴向对称边界条件，垂直于 Z轴且 Z=0 的平面施加Z轴向对称边界条件，右侧面施加100MPa的拉力。&-1 单层板边界条件及加载情况本例中单位系统为mm、 MPa。字母12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455&.2.2 UMAT 用户子程序UMAT 用户子程序UMA T-Stress.for的全部代码如下，C 及“ !

5、”之后为注释内容。SUBROUTINUEMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1 RPL,DDSDDT,DRPLDE,DRPLDT,2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)CINCLUDE'ABA_PARAM.INC'CCHARACT*E8R0 CMNA

6、MEDIMENSIO NSTRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),4 JSTEP(4)DIMENSIO NEG(6), XNU(3,3), STRAND(6), C(6,6), STRESS0(6)C*C EGE1,E2,E3,G12,G13,G23C

7、 XNUNU12 ， NU21,NU13， NU31， NU23， NU32C STRANDSTRAINT AT THE END OF THE INCREMENTC C6X6 STIFFNESS MATRIXC STRESS0STRESS AT THE BEGINNING OF THE INCREMENTC*C INITIALIZE XNU & C MATRIXXNU=0C=0C GET THE MATERIAL PROPERTIES-ENGINEERING CONSTANTSEG(1) = PROPS(1)EG(2) = PROPS(2)EG(3) = EG(2)XNU(1,2) =

8、 PROPS(3)XNU(2,1) = XNU(1,2)*EG(2)/EG(1)XNU(1,3) = XNU(1,2)XNU(3,1) = XNU(1,3)*EG(3)/EG(1)XNU(2,3) = PROPS(4)XNU(3,2) = XNU(2,3)*EG(3)/EG(2)EG(4) = PROPS(5)EG(5) = EG(4)EG(6) = PROPS(6)C*!E1,YOUNG'S MODULUS IN DIRECTION 1!E2,YOUNG'S MODULUS IN DIRECTION 2!E3,YOUNG'S MODULUS IN DIRECTION

9、3!POISON'S RATIO POI_12!POISON'S RATIO POI_21!POISON'S RATIO POI_13!POISON'S RATIO POI_31!POISON'S RATIO POI_23!POISON'S RATIO POI_32!G12,SHEAR MODULUS IN 12 PLANE!G13,SHEAR MODULUS IN 13 PLANE!G23,SHEAR MODULUS IN 23 PLANEC FILL THE 6X6 STIFFNESS MATRIX C(6,6)RNU = 1/(1-XNU(

10、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)C STIFFNESS MATRIX C(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

11、(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)*RNU56C(3,2) = C(2,3)57 C*58 C CALCULATE STRAIN59 DOI = 1, 660 STRAND(I) = STRAN(I)+DSTRAN(I)61 ENDDO62 C CALCULATE STRESS63 DOI = 1, 664 STRESS0(I) = STRESS(I)65 STRESS(I)

12、= 066 DOJ = 1, 667 STRESS(I) = STRESS(I)+C(I,J)*STRAND(J)68 ENDDO69 ENDDO70 C CALCULATE SSE71 DOI = 1, 672 SSE = SSE+0.5*(STRESS0(I)+STRESS(I)*DSTRAN(I)73 ENDDO74 C*75 C UPDATE DDSDDE76 DOI = 1, 677 DOJ = 1, 678 DDSDDE(I,J) = C(I,J)79 ENDDO80 ENDDO81 RETURN82 END第 1 到 14行及第 81、 82行为 UMAT 子程序固定格式，其中，

13、第1 到 5 行括号内的变量为UMAT 子程序中可以使用的变量，第10 到14行定义各变量数组的维数和长度。部分主要变量的含义如表&-2 所表 &-2 UMA T 部分变量名及其含义STRESS增量步开始时的应力(S11, S22, .)， 用增量步结束时的应力计算结果对其更新STATEV(NSTATV)状态变量(状态变量个数)，如果在材料中定义了状态变量，则在UMA T 中需要对其更新STRAN增量步开始时的应变(E11, E22, .)DSTRAN当前增量步的应变增量( E11, E22, .)NDI, NSHR, NTENS应力、应变的个数，NDI 为正应力或正应变的个数

14、，NSHR 为剪应力或剪应变的个数，NTENS=NDI+NSHRPROPS, NPROPS材料参数、材料参数的个数DDSDDE雅克比矩阵， / SSE, SPD, SCD特定的弹性应变能、塑性耗散、蠕变耗散，只对能量输出有影响，对其他计算结果无影响，在 UMAT 中需要对其更新单元特征长度CELENT第 15 到 83行为用户自己编写的固定格式的Fortran 程序， 用于计算刚度矩阵、应力、应变能、雅克比矩阵。由于本例中没有使用状态变量，因此不需要更新STATEV，只需要更新STRESS、 DDSDDE 和SSE 即可。第 16 行定义了5 个数组， 其中EG、 STRAND 、 STRES

15、S0为一维数组， 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 行，读取材料常数，计算泊松比。泊松比的计算公式如下。第 42 到 56 行，计算刚度矩阵。刚度矩阵

16、的计算公式如下。C11E1(12332 )C22E2(11331 )C33E3(11221)C44G12C55G13C66G23E2 ( 121332 )E3 ( 131223)E3 ( 232113)(&-2)C12E1(213123)C13E1(312132)C23E2 ( 321231)C21C12C31C13C23112211331123322132132对于本例所用材料，由于E2=E3， G12=G13， 12= 13，所以 21= 31、 23= 32， 可以将式(&-2) 化简后在UMAT 中计算刚度矩阵。本例为保证UMAT-Stress.for 的可读性，刚度矩

17、阵的计算没有化简，直接按照式(&-2) 编写。第 59 到 61 行，计算当前增量步的应变，计算公式如式(&-3) ，式中上标表示增量步序号。i1i i 1(&-3)第 63 到 69 行，将当前分析步开始时的应力STRESS 的值赋给STRESS0，然后计算当前增量步的应力并赋给STRESS，当前增量步的应力计算如式(&-4) ，式中上标表示增量步序号，应力i 和应变i为列向量。UMAT 中 的 剪应变为工程剪 应变。 i C i(&-4)第 71 到 73 行，计算应变能，如式(&-5) 所示，式中上标表示增量步序号；下标表示应力、应变增量的

18、分量序号，其中下标为1 、 2、 3表示正应力、正应变增量，4、 5、 6 表示剪应力、剪应变增量。6Ei 1Ei 0.5 ( jiji 1)ij 1(&-5)j1第 76 到 80 行为更新雅克比矩阵。由于没有对刚度矩阵没变化，应力与应变的关系如式(&-4) 所示， 所以应力增量与应变增量的关系如式 (&-6) 所示，所以雅克比矩阵就等于刚度矩阵。i C i i(&-6)iC第 76 到 80 行可简写为一行，如下所示：DDSDDE(1:6,1:6) = C(1:6,1:6) ，或 DDSDDE = C&.2.3 复合材料单层板应力分析1 、创建部件及

19、划分网格创建部件在 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 键或单击鼠标中

20、键。单击提示区的Done 或鼠标中键，在 Edit Base Extrusion 对话框 Depth 后面输入0.15，单击OK 完成。划分网格在环境栏Module 后面选择Mesh， 进入Mesh模块。环境栏中Object选择 Part: Lam-C 。单击工具区的(Seed Part)， 在 Global Seeds 对话框中Approximateglobal size 后面输入0.5，单击OK 。单击工具区的(Mesh Part)，单击提示区的Yes或鼠标中键，完成网格划分，如图&-2 所示，板的厚度方向只划分为1 层单元。输入数据时，每输入完一行后按回车(Enter)键， 光标

21、会自动移到下一行。也可以通过右键快捷菜单添加或删除一行。本例UMAT子程序较简单，不需要使用状态变量，因此在材料UMat-T700 中图 &-2 划分网格单击工具区的(Assign Element Type) ， 在 Element Type 对话框中，选择依次选择Standard， Linear， 3D Stress，在Hex 标签页中勾选Reduced integration ，在 Element Controls 区域 Hourglass control 选择 Enhanced，即选择C3D8R 单元，单击OK 完成。单击工具区的(Assign Stack Direction)

22、，在视图区选择部件平行于 X-Y 平面的面，单击鼠标中键或提示区的Yes 完成。2、创建材料并给部件赋材料属性创建材料在环境栏Module 后面选择Property，进入Property 模块。单击工具区的(Create Material) ， 在 Edit Material 对话框中，Name后面输入UMat-T700 ；单击General User Material ，在 User Material区域中 Data 区域的 Mechanical Constants 一栏依次输入114000， 8610，0.3， 0.45， 4160， 3000，单击OK 完成。单击工具区的(Create

23、Material) ， 在 Edit Material 对话框中，Name后面输入Mat-T700 ；单击 Mechanical Elasticity Elastic，在Elastic区域中Type 选择Engineering Constants，在 Data 区域输入从左到右依次输入114000，8610， 8610， 0.3， 0.3，0.45， 4160，4160，3000，单击 OK 完成。材料 UMat-T700 用于 UMAT 用户子程序，Mat-T700 用于做对比没有定义Depvar。分析。给部件赋材料属性单击工具区的(Create Composite Layup) ，在打开的

24、对话框中Name 使用默认名称，Initial ply count 后面输入1， Element Type 选择Solid，单击Continue. ；在Edit Composite Layup 对话框中，Layup在定义复合材料铺层时，视图区部件上会显示铺层方向，在 EditComposite Layup 对话框中的Display 标签页中可以设置所需显示的方向，在视图区部件上白色箭头及字母S 表示StackingDirection 。Orientation 区域的 Definition 选择 Coordinate system，单击Definition下一行的(Create Datum CS

25、YS)。在 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 选择 Elementdirection 3 ， Rotation axis 选择 Axis 3

26、。在 Plies 标签页中，双击Region，在视图区选择部件后单击鼠标中键；右单击Material ，在快捷菜单中单击Edit Material. ，在 SelectMaterial 对话框中选择UMat-T700 ， 单击OK； 右单击 Element RelativeThickness，在快捷菜单中单击Edit Thickness.，在Thickness 对话框中 Specify Value 后面输入1 后单击 OK； 在 Rotation Angle 一栏输入0；Integration Points 使用默认的1，单击 OK 完成。3 、装配在环境栏Module 后面选择Assembl

27、y，进入 Assembly 模块。单击工具区的(Create Instance)， 在 Create Instance对话框中选择Parts: Lam-C，单击 OK 完成。4、创建分析步、设置输出变量创建分析步在环境栏Module 后面选择Step，进入Step 模块。单击工具区的(Create Step)，在Create Step 对话框中，在Initial分析步之后插入Static, General 分析步，单击Continue. ；在Edit Step对话框中使用默认设置，单击OK 完成。设置输出变量单击工具区的(Field Output Manager) ，在 Field Outpu

28、t RequestsManager 对话框中，选中 F-Output-1 ， 单击 Edit. ； 在 Edit Field OutputRequest对话框中，设置如图&-3 所示， 输出整个模型最后一个增量步的S、 E 、 U ，单击 OK 完成。图 &-3 场输出变量设置单击工具区的(History Output Manager) ，在 History OutputRequests Manager 对话框中，选中H-Output-1 ，单击Edit. ；在EditHistory Output Request 对话框中，设置输出整个模型的内能和应变能，即 ALLIE 和 A

29、LLSE ，单击 OK 完成。5、创建边界条件及施加载荷创建边界条件在环境栏Module 后面选择Load，进入 Load 模块。单击工具区的( Create Boundary Condition ) ， 在 Create BoundaryCondition 对话框中，Name 后面输入BC-X ， Step 选择 Initial ， Category边界条件及 加载情况见图 &-1 。选择 Mechanical ， Types for Selected Step 选择 Symmetry. ，单击 Continue.；在视图区选择装配实例Lam-C-1 左侧端面，即垂直于X轴且 X=0

30、 的侧面，单击鼠标中键或提示区的Done，在Edit BoundaryCondition 对话框中选择XSYMM ，单击 OK 完成。类似操作，选择Lam-C-1 的下侧端面，即垂直于Y 轴且 Y=0 的侧面，创建边界条件BC-Y ，边界类型为YSYMM ；选择 Lam-C-1 垂直于 Z 轴且 Z=0 的侧面，创建边界条件BC-Z，边界类型为ZSYMM 。施加载荷单击工具区的(Create Load)，在 Create Load 对话框中，Name使用默认的Load-1 ， Step 选择Step-1， Category 选择 Mechanical ， Typesfor Selected S

31、tep 选择Pressure， 单击Continue. ； 在视图区选择Lam-C-1的右侧端面，即垂直于X 轴且 X=15 的侧面，单击鼠标中键，在EditLoad 对话框中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

32、 对话框的General标签页中，单击 User subroutine file 后面的(Select.)， 在相应路径下找到并选择UMA T-Stress.for 文件； 单击 Edit Job 对话框中的OK 完成。提交分析在 Job Manager 对话框中，选中Job-Lam-Stress-Umat 分析作业，单击 Submit 提交计算。当 Job-Lam-Stress-Umat 的状态 ( Status) 由 Running 变为 Completed时，计算完成，单击Results进入可视化后处理模块。保存模型单击工具栏的File 工具条中的( Save Model Databas

33、e) ， 在 SaveModel 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 ，

34、 单击 OK； 在 Edit CompositeLayup 对话框中单击OK 完成。在 Edit Job对 话框的 Parallelization 标 签页中可以设置 多核并行计算。8 、再次创建分析作业并提交分析创建分析作业在环境栏Module 后面选择Job，进入Job 模块。单击工具区的( Job Manager)，在Job Manager 对话框中单击Create.；在Create Job 对话框中，Name 后面输入Job-Lam-Stress，Source 选择 Model-1 ，单击Continue. ；在Edit Job 对话框使用默认设置，单击OK 完成。提交分析在 Job

35、Manager 对话框中，选中Job-Lam-Stress 分析作业，单击Submit 提交计算。在 Edit Job对 话框的 Parallelization 标 签页中可以设置 多核并行计算。当 Job-Lam-Stress 的状态 ( Status) 由 Running 变为 Completed 时，Results 进入可视化后处理模块。保存模型单击工具栏的File 工具条中的(Save Model Database)保存模型。9、可视化后处理显示云图在视图区显示Job-Lam-Stress-Umat.odb 。长按工具区的(Plot Contours on Deformed Shape

36、) ， 显示隐藏工具(Plot Contours on Undeformed Shape) ，在 Field Output 工具条S11。单击菜单栏Result Section Points. ，打开SectionPoints 对话框，Selection method 选择 Plies，在 Plies 区域选择PLY-1 ，Apply ，在视图区显示该铺层的S11 应力云图。相同操作，可以显示各铺层的各应力分离的云图。Job-Lam-Stress-Umat.odb 和 Job-Lam-Stress.odb的各应力分量的云 &-4 所示。(a) S11， Job-Lam-Stress-U

37、mat(b) S11 ， Job-Lam-Stress(d) S22， Job-Lam-Stress(c) S22， Job-Lam-Stress-Umat读者在阅 读 Abaqus 帮助 文 件 Abaqus Example Problems Guide 中 1.4.6 Failure of blunt notched fiber metal laminates 一例的 exa_fml_ortho_ damage_umat.f 文件时注意 SSE 的计算公(f) S12 ， Job-Lam-Stress(e) S12， Job-Lam-Stress-Umat&-4 各应力分量对比单击

38、菜单栏Result History Output.，可以直 接 打 开 History Output 对话框。输出应变能在视图区显示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 对话框中选择X

39、YData-1 ， 单击Edit.； 在 Edit XYData 对话框中可以看到应变能数据，分析结束时整个模型的应变能为9.671mJ。同样操作，观察 Job-Lam-Stress.odb 中整个模型的应变能数据，同样为9.673mJ。通过对比，验证了 UMAT-Stress.for中第 72行 SSE计算的准确性。&.2.4 Inp 文件解释本例中 Job-Lam-Stress-Umat.inp 节选如下。*Heading* * Job name: Job-Lam-Stress-Umat Model name: Model-1* * Generated by: Abaqus/CAE

40、 6.14-1*Preprint, echo=NO, model=NO, history=NO, contact=NO* * * PARTS* * 部件 Lam-C-1 的节点、单元数据* Part, name=Lam-C*Node1,15.,10.,0.15 节点编号及坐标* 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* 定义复合材料铺层*Solid Section, el

41、set=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* * 单元控制，使用沙漏控制

42、* 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=6 114000.,8610., 0.3,

43、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: Press

44、ure*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 Ste

45、p&.2.5 应用 UMAT 子程序应力分析小结本例所介绍的应力分析是UMAT 用户子程序最简单的应用，读者在了解 UMAT 子程序的基础上，可以结合相应的强度理论、刚度折减方法对 UMAT 子程序中的刚度矩阵、雅克比矩阵进行计算和更新，从而达到刚度折减和渐进损伤强度分析的目的。&.3 实例二：UMAT 用户子程序渐进损伤强度分析本例在实例一使用的UMAT 子程序的基础上，通过增加失效判定， 对 UMAT 子程序中的刚度矩阵进行折减，达到渐进损伤强度分析的目的。为了使本例的UMAT 子程序简单易读，本例对复合材料单层板进行纤维方向的拉伸强度分析，仅考虑纤维方向拉伸破坏。本例使用

46、本章实例一的模型，稍做修改，用于本例的强度分析。6.3.1 问题描述本例中单位系统为mm、 MPa。复合材料单层板的3D 实体模型的几何尺寸及边界条件同本章的实例一的模型。将右侧面100MPa 的拉力替换为X 轴正方向0.4mm 的位移载荷。6.3.2 UMAT 用户子程序本例使用的UMAT 用户子程序UMAT-Strength.for 的全部代码如下，字母C 及“!”之后为注释内容。1 SUBROUTINUEMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,2 1RPL,DDSDDT,DRPLDE,DRPLDT,3 2STRAN,DSTRAN,TIME,DTIME,T

47、EMP,DTEMP,PREDEF,DPRED,CMNAME,4 3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,5 4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)6C7INCLUDE'ABA_PARAM.INC'8C9CHARACTE*8R0 CMNAME10DIMENSIONSTRESS(NTENS),STATEV(NSTATV),111DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),122STRAN(NTE

48、NS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),133PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),14 4JSTEP(4)1516DIMENSIONSTRESS0(6),STRAND(6),C(6,6),CD(6,6),DCDE(6,6)17 PARAMETE(RALPHA = 1000.0, LAMBDA = 0.0, DMAX = 0.9999, DRND = 3)18 C*19 C STRESS0STRESS AT THE BEGINNING OF THE INCREMENT

49、20 C STRANDSTRAIN AT THE END OF THE INCREMENT21 C C6X6 STIFFNESS MATRIX22 C CD6X6 DAMAGED STIFFNESS MATRIX23 C DCDED CD/D E24252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384C STATEV(1)DAMAGE VARIABLE DC*C GET THE MATERIAL PROPER

50、TIESE1 = PROPS(1)E2 = PROPS(2)XNU12 = PROPS(3)XNU21 = XNU12*E2/E1XNU23 = PROPS(4)G12 = PROPS(5)G23 = PROPS(6)STH = PROPS(7)C*!E1,YOUNG'S MODULUS IN DIRECTION 1!E2 = E3,YOUNG'S MODULUS IN DIRECTION 2 & 3!POISON'S RATIO POI_12,XNU13 = XNU12!POISON'S RATIO POI_21,XNU31 = XNU21!POISO

51、N'S RATIO POI_23,XNU32 = XNU23!G12 = G13,SHEAR MODULUS IN 12 & 13 PLANE!G23,SHEAR MODULUS IN 23 PLANE!FAILURE STRESS IN 1 DIRECTION IN TENSIONC STIFFNESS MATRIX C(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) =

52、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)C CALCULATE THE STRAIN AT THE END OF THE INCREMENTDOI = 1, 6STRAND(I) = STRAN(I) + DSTRAN(I)ENDDOC*C CALCULATE THE FAILURE COEFFICIENTSTRANF = STH/E1IF (STRAND(1) > 0) THENF = STRAND(1)/STRANFELSEF = 0ENDIFC CALCULATE D,DAMAGE VARIABLED = 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