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

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

《UMAT子程序在复合材料强度分析中的应用.doc》由会员分享，可在线阅读，更多相关《UMAT子程序在复合材料强度分析中的应用.doc（27页珍藏版）》请在冰豆网上搜索。

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

本例使用UMAT用户子程序进行复合材料单层板的应力分析和渐进损伤压缩强度分析，介绍UMAT用户子程序编写方法及在Abaqus/CAE中的设置。

本章使用最大应变强度理论作为复合材料单层板的失效准则，相应的Fortran程序简单易读，便于理解UAMT子程序的工作原理。

知识要点：

Ø强度分析

ØUMAT用户子程序

Ø最大应变理论

Ø刚度折减

讲师：

孔祥宏

版本：

Abq

难度：

关键词：

强度分析，UMAT

&.1本章内容简介

本章通过两个实例介绍UMAT用户子程序在复合材料单层板的应力分析和强度分析中的应用。

在第一个实例中，对一个简单的复合材料单层板进行应力分析，UMAT子程序主要计算应力，不进行强度分析，本例用于验证UMAT子程序的计算精度。

在第二个实例中，对复合材料单层板进行渐进损伤强度分析，UMAT子程序用于应力计算、强度分析和刚度折减。

本章所用复合材料为T700/BA9916，材料属性如表&-1所示。

表&-1T700/BA9916材料属性

参数

值

强度

值

E1/GPa

114

XT/MPa

2688

E2/GPa

XC/MPa

1458

E3/GPa

YT/MPa

μ12

YC/MPa

236

μ13

ZT/MPa

μ23

ZC/MPa

175

G12/GPa

SXY/MPa

136

G13/GPa

SXZ/MPa

136

G23/GPa

SYZ/MPa

&.2实例一：

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及“!

”之后为注释内容。

1SUBROUTINEUMAT（STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

21RPL,DDSDDT,DRPLDE,DRPLDT,

32STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,

43NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,

54CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC）

6C

7INCLUDE''

8C

9CHARACTER*80CMNAME

10DIMENSIONSTRESS（NTENS）,STATEV（NSTATV）,

111DDSDDE（NTENS,NTENS）,DDSDDT（NTENS）,DRPLDE（NTENS）,

122STRAN（NTENS）,DSTRAN（NTENS）,TIME

（2）,PREDEF

（1）,DPRED

（1）,

133PROPS（NPROPS）,COORDS（3）,DROT（3,3）,DFGRD0（3,3）,DFGRD1（3,3）,

144JSTEP（4）

15

16DIMENSIONEG（6）,XNU（3,3）,STRAND（6）,C（6,6）,STRESS0（6）

17C****************************

18CEG.....E1,E2,E3,G12,G13,G23

19CXNU.....NU12，NU21,NU13，NU31，NU23，NU32

20CSTRAND.....STRAINTATTHEENDOFTHEINCREMENT

21CC.....6X6STIFFNESSMATRIX

22CSTRESS0.....STRESSATTHEBEGINNINGOFTHEINCREMENT

23C****************************

24CINITIALIZEXNU&CMATRIX

25XNU=0

26C=0

27CGETTHEMATERIALPROPERTIES---ENGINEERINGCONSTANTS

28EG

（1）=PROPS

（1）!

E1,YOUNG'SMODULUSINDIRECTION1

29EG

（2）=PROPS

（2）!

E2,YOUNG'SMODULUSINDIRECTION2

30EG（3）=EG

（2）!

E3,YOUNG'SMODULUSINDIRECTION3

31XNU（1,2）=PROPS（3）!

POISON'SRATIOPOI_12

32XNU（2,1）=XNU（1,2）*EG

（2）/EG

（1）!

POISON'SRATIOPOI_21

33XNU（1,3）=XNU（1,2）!

POISON'SRATIOPOI_13

34XNU（3,1）=XNU（1,3）*EG（3）/EG

（1）!

POISON'SRATIOPOI_31

35XNU（2,3）=PROPS（4）!

POISON'SRATIOPOI_23

36XNU（3,2）=XNU（2,3）*EG（3）/EG

（2）!

POISON'SRATIOPOI_32

37EG（4）=PROPS（5）!

G12,SHEARMODULUSIN12PLANE

38EG（5）=EG（4）!

G13,SHEARMODULUSIN13PLANE

39EG（6）=PROPS（6）!

G23,SHEARMODULUSIN23PLANE

40C****************************

41CFILLTHE6X6STIFFNESSMATRIXC（6,6）

42RNU=1/（1-XNU（1,2）*XNU（2,1）-XNU（1,3）*XNU（3,1）-

431XNU（3,2）*XNU（2,3）-2*XNU（1,3）*XNU（2,1）*XNU（3,2））

44CSTIFFNESSMATRIXC（6,6）

45C（1,1）=EG

（1）*（1-XNU（2,3）*XNU（3,2））*RNU

46C（2,2）=EG

（2）*（1-XNU（1,3）*XNU（3,1））*RNU

47C（3,3）=EG（3）*（1-XNU（1,2）*XNU（2,1））*RNU

48C（4,4）=EG（4）

49C（5,5）=EG（5）

50C（6,6）=EG（6）

51C（1,2）=EG

（1）*（XNU（2,1）+XNU（3,1）*XNU（2,3））*RNU

52C（2,1）=C（1,2）

53C（1,3）=EG

（1）*（XNU（3,1）+XNU（2,1）*XNU（3,2））*RNU

54C（3,1）=C（1,3）

55C（2,3）=EG

（2）*（XNU（3,2）+XNU（1,2）*XNU（3,1））*RNU

56C（3,2）=C（2,3）

57C****************************

58CCALCULATESTRAIN

59DOI=1,6

60STRAND（I）=STRAN（I）+DSTRAN（I）

61ENDDO

62CCALCULATESTRESS

63DOI=1,6

64STRESS0（I）=STRESS（I）

65STRESS（I）=0

66DOJ=1,6

67STRESS（I）=STRESS（I）+C（I,J）*STRAND（J）

68ENDDO

69ENDDO

70CCALCULATESSE

71DOI=1,6

72SSE=SSE+*（STRESS0（I）+STRESS（I））*DSTRAN（I）

73ENDDO

74C****************************

75CUPDATEDDSDDE

76DOI=1,6

77DOJ=1,6

78DDSDDE（I,J）=C（I,J）

79ENDDO

80ENDDO

81RETURN

82END

第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+NSHR

PROPS,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=C

UMAT中的剪应变为工程剪应变。

&.复合材料单层板应力分析

1、创建部件及划分网格

l创建部件

在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完成。

l划分网格

在环境栏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完成。

2、创建材料并给部件赋材料属性

l创建材料

在环境栏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。

l给部件赋材料属性

单击工具区的（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。

3、装配

在环境栏Module后面选择Assembly，进入Assembly模块。

单击工具区的（CreateInstance），在CreateInstance对话框中选择Parts:

Lam-C，单击OK完成。

4、创建分析步、设置输出变量

l创建分析步

在环境栏Module后面选择Step，进入Step模块。

单击工具区的（CreateStep），在CreateStep对话框中，在Initial分析步之后插入Static,General分析步，单击Continue...；在EditStep对话框中使用默认设置，单击OK完成。

l设置输出变量

单击工具区的（FieldOutputManager），在FieldOutputRequestsManager对话框中，选中F-Output-1，单击Edit...；在EditFieldOutputRequest对话框中，设置如图&-3所示，输出整个模型最后一个增量步的S、E、U，单击OK完成。

图&-3场输出变量设置

单击工具区的（HistoryOutputManager），在HistoryOutputRequestsManager对话框中，选中H-Output-1，单击Edit...；在EditHistoryOutputRequest对话框中，设置输出整个模型的内能和应变能，即ALLIE和ALLSE，单击OK完成。

5、创建边界条件及施加载荷

l创建边界条件

在环境栏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。

l施加载荷

单击工具区的（CreateLoad），在CreateLoad对话框中，Name使用默认的Load-1，Step选择Step-1，Category选择Mechanical，TypesforSelectedStep选择Pressure，单击Continue...；在视图区选择Lam-C-1的右侧端面，即垂直于X轴且X=15的侧面，单击鼠标中键，在EditLoad对话框中Magnitude后面输入-100，单击OK完成。

6、创建分析作业并提交分析

l创建分析作业

在环境栏Module后面选择Job，进入Job模块。

单击工具区的（JobManager），在Job