cohesive界面单元VUMAT子程序

1 基于cohesive单元的裂缝扩展模拟

1.1 cohesive单元简介

       cohesive单元（粘性单元）定义： 通过预置裂纹边（面）的方式来模拟二维（三维）裂纹，即在预判的裂纹区域加入一层厚度为0的cohesive单元。

        cohesive单元简界：如图1所示，cohesive单元由顶面，中面和底面组成。cohesive模型通过损伤起始准则以及演化方法来判断其损伤情况。损伤发生在cohesive单元上，当满足损伤判据时cohesive单元开始进入损伤，根据预先定义的损伤演化模型，当单元损伤累积到完全失效时cohesive单元的中层;一分为二，使物体产生几何上的不连续，从而形成开裂。

图1  cohesive单元示意图

1.2 断裂力学原理

      cohesive单元的损伤时由牵引-分离(traction-separation)定理描述的，在宏观上牵引力与位移的关系如图2所示。

ABAQUSprocess=/format,webp/resize,w_760" https: www.jishulink.com major abaqus" style="text-decoration: none; cursor: text; color: #333;" data-initial-src="https://img.jishulink.com/msimage/202303/46bee398bf9e0e51edbeb00a84b6b667.png">ABAQUS/format,webp/resize,w_760" data-initial-src="https://img.jishulink.com/msimage/202303/46bee398bf9e0e51edbeb00a84b6b667.png">

图2   牵引力-牵引位移示意图

      在损伤前，应力与应变满足如下关系:ABAQUSlass="p1">

        公式中，上标0单元分层损伤起始位移，f代表最终开裂位移，max代表加载过程中最大位移，并且满足损伤不可逆性。

       损伤起始准则用于判断材料是否开始损伤。常见的损伤判据如下表。

       当材料开始满足损伤起始准则之后，表现出整体刚度退化的现象。并且损伤演化定义了刚度退化的规律和裂缝的扩展，同时损伤演化分为位移和能量两种模式：位移即为图2的横轴，能量为图中曲线围成的面积。

2 ABAQUS中关于CohesivABAQUS和设置

2.1 牵引-分离准则

       通常假定在初始阶段粘结滑移单元的牵引- 分离曲线是线形变化的，一旦达到破坏条件，材料就会按照用户自定义的损伤演化准则发生破坏，图3就是典型的牵引-分离曲线。

图3 混合模式下cohesive单元牵引力-牵引位移示意图

2.2 弹性模量

      以三维cohesive单元为例，其牵引应力张量为t

式中:

tn--法向牵引应力(局部坐标3轴方向)

Ts,tt--切线牵引应力(局部坐标1, 2轴方向)

      ABAQUS中Cohesive单元的局部坐标方向如图2所示。

      以δn、δs、δt分别表示对应方向的位移分量，以T0表示Cohesive单元的初始厚度，则其名义应变分量可表示为:

      其材料本构关系可以用下式表达:

      在ABAQUS中指定材料参数时，如选择Traction，如图3所示，则刚度矩阵E的非对角线元素全部为0，即:

      此时仅需填写材料三个方向的弹性模量Enn、ESS、Ett.

      当在指定材料弹性类型时选择Coupled Traction，如图4所示，则刚度矩阵的的非对角元素不为0，此时由于刚度矩阵的对称性，需要填入除对角元素外的三个非对角元素Ens、Ent和Est，此外在此输入的弹性模量值受指定的材料初始厚度的影响，详见2材料厚度的选择。

   

图3: Traction选项                  

图4: Coupled Traction 选项

2.3 损伤准则

2.4 损伤演化

此部分主要为损伤演化因子D。

2.4.1 材料厚度的选择

假定由T-S曲线所得的材料弹性阶段的曲线斜率为Kp ，软件中指定的材料厚度为heff，以法向为例，在其法向位移为δn的情况下其法向应力为：

由材料力学可知：

结合上述公式可得：

      所以在ABAQUS中填写材料初始厚度时，如图5所示,有“Use analysis default" 、“Use nodal coordinates”和“Specify” 三种选择。当选择“Use analysis default”时，软件默认的材料初始厚度为1.0。此时E= K，软件输出结果中，应变=位移；当用户按实际材料厚度hgeo建立模型，并选择“Use nodal coordinates”时，此时材料弹性模量填写时，E= Khgeo；当选择“Specify" ,并由用户自己指定材料厚度时，假定指定的材料厚度为hsp，此时材料弹性模量填写时，E=Khsp。三种方式对应的输出变量中仅Cohesive单元的位移不同，其他变量(如，应力、应变)均无差别。

图5:材料初始厚度的选择

3 ABAQUS中Cohesive单元建模方法讲解

       应用 cohesive 单元模拟复合材料失效，包括两种模型：一种是基于 traction-separation描述；另一种是基于连续体描述。其中基于 traction-separation 描述的方法应用更加广泛。而在基于 traction-separation 描述的方法中，最常用的本构模型为图2所示的双线性本构模型。它给出了材料达到强度极限前的线弹性段和材料达到强度极限后的刚度线性降低软化阶段。注意图中纵坐标为应力，而横坐标为位移，因此线弹性段的斜率代表的实际是cohesive 单元的刚度。曲线下的面积即为材料断裂时的能量释放率 。因此在定义 cohesive 的力学性能时，实际就是要确定上述本构模型的具体形状：包括刚度、极限强度、以及临界断裂能量释放率，或者最终失效时单元的位移。常用的定义方法是给定上述参数中的前三项，也就确定了 cohesive 的本构模型。Cohesive 单元可理解为一种准二维单元，可以将它看作被一个厚度隔开的两个面，这两个面分别和其他实体单元连接。Cohesive 单元只考虑面外的力，包括法向的正应力以及 XZ ，YZ 两个方向的剪应力。

cohesive单元的刚度

基于 traction-separation 模型的界面单元的刚度可以通过一个简单杆的变形公式来理解:

其中 L为杆长，E为弹性刚度， A为初始截面积，P为载荷。公式（1）又可以写成:

其中S=P/A为名义应力，K=E/L为材料的刚度。

为了更好的理解K，我们把K=E/L写成:

      这里我们用L'来代替1,其中L可以理解为建模厚度,即建模时cohesive interface的几何厚度; L'为实际厚度，即cohesive interface的真实厚度，这个厚度在cohesive section中定义。E/L可以理解为几何刚度，即模型中

cohesive interface所具有的刚度;为cohesive interface 的真实刚度。当L'为1时，计算界面刚度就采用几何刚度E/L，当L'为0.001时，计算时界面刚度变为1000E/L。举个小例子，如果界面的实际厚度为0.01，而在建模时就是按照这个厚度建立的，在定义material-section 时又specify这层的厚度为0.01，实际上就等于把界面刚度提高了2个数量级，模拟结果当然是不对的，这时定义section时应采用默认厚度1。ABAQUS在cohesive建模中使用了很“人性化”的设计，实际问题中界面可能很薄，有的只有0.001mm,甚至更小。有些问题cohesive单元的interface 还可能是0厚度(比如crack问题)，而相对来说整体模型也许很大，如果不引入这两个厚度，我们就要在很大的模型中去创建这个很小的界面这是一个很麻烦的事情。引入这两个厚度,在建模时我们就可以用有限的厚度来代替这个很小的界面厚度，只要在section中定义这个L'就好了。(注:以上大部分内容来自仿真论坛:再议cohesive应用中对于一-些参数的理解)

4 一个解释

       另外有个我的经验公式:大体上energy > 0.5*( damage initiation) ^2/ (stiffness)这个公式不难理解，就是锐角三角形的总面积大于一条侧边下的面积，将traction-separation law画成图线你就一目了然了。不过根据不同的法则，会稍微有些区别的。”----以上的话引自dava的个人空间，这里我想解释下这个不等式，有些新手可能一下还看不明白。damage initiation 为开始破坏时的应力，即三角形的高； stiffness 为刚度，也就是斜率，即tan q ；所以侧边三角形的底边为damage initiation/stiffness ，0.5* (damage initiation) ^2/(stiffness)即为侧边下的三角形面积。实际上能量还要大于这个侧边下三角形的面积很多,因为斜率一般都很大。

5 关于材料参数

        定义cohesive的材料时，要填入材料的参数，这些材料参数是材料固有的特性，与几何没有关系，所以放心大胆的填入吧。材料参数是由试验得到的，如果不能做实验(多数情况如此)，就去查国际上相关的文献吧，数据甚至比你自己做试验都要详细，在填入数据时要注意单位的统一。再说句，断裂能为单位面积上的能量，如你的单位选取N (力的单位)长和M(度单位)，那么能量的单位为N/M。

下面举例来说明cohesive单元刚度的设置过程，以ABAQUS6.9为例:

       进入property界面，点击Material→Creat，在弹出的Edit Material对话框中，可以编辑新创建的cohesive材料的名称，然后点击Mechanical→Elasticity→Elastic→Traction，在空格中输入相应的刚度。

6 损伤准则

6.1 初始损伤准则

        初始损伤对应于材料开始退化，当应力或应变满足于定义的初始临界损伤准则，则此时退化开始。Abaqus 的Damage for traction separation law中包括: Quade Damage（二次应变）、Maxe Damage（最大应变）、Quads Damage（二次应力-常用）、Maxs Damage（最大应力）、Maxpe Damage、Maxps Damage六种初始损伤准则，其中前四种用于一般复合材料分层模拟，后两种主要是在扩展有限元法模拟不连续体(比如crack问题)问题时使用。

使用图2所示的双线本构模型，其中：tn0、ts0、tt0分别代表纯Ⅰ型、纯Ⅱ型或纯Ⅲ型破坏的最大名义应力，εn0、εs0、εt0代表相应的最大名义应变，当定义界面单元的初始厚度为1时，则名义应变等于与之相对应的相对位移δn、δs、δt  。

Quade Damage为二次名义应变准则:当名义应变比的平方和等于1时，损伤开始。

Maxe Damage为最大名义应变准则:当任何一个名义应变的比值达到1时，损伤开始。

Quads Damage为二次名义应力准则:当各个方向的名义应变比的平方和等于1时，损伤开始。

Maxs Damage为最大名义应力准则:当任何一个名义应力比值达到1时，损伤开始。

       Edit Material对话框中，点击Mechanical-→Damage for Traction Separation Laws然后根据自己的需要点击相应的损伤准则。其中最常用是Quads Damage.

6.2 损伤演化规律

        选择了初始损伤准则之后，然后点击Sub options - Damage Evolution窗口，其中Type包括Displacement和Energy , Displacement为基于位移的损伤演化规律,而Energy为基于能量的损伤演化规律。Softening 中包括Linear, Exponential 及Tabular三种刚度退化方式....Evolution中的所有的选项都是用来确定单元达到强度极限以后的刚度降阶方式。一般常用:以能量来控制单元的退化，即Type -> Energy；线性软化模型，即Softening → Linear，Degradation → Maximum ; Mixed mode behavior -> BK，Mode mix ratio ->Energy，并选中Power。