修正内聚区长度计算公式在冰I 型裂纹扩展中的应用

倪宝玉，徐莹，黄其，尤嘉，薛彦卓

1 哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001

2 中国船舶及海洋工程设计研究院，上海 200011

3 中国船级社 江苏分社，江苏 南京 210011

0 引 言

结构物与冰的相互碰撞研究在冰区结构物的设计制造、极地地区资源的开发和利用，以及冰区船舶和结构物的安全航行及作业等领域均受到普遍重视。目前，在结构物−冰碰撞的研究中使用较广泛的方法包括离散元法[1-2]和有限元法[3]。其中，离散元法用于冰裂纹扩展的研究具有一定的优势，但对于海冰变形的研究不够充分；有限元法用于海冰与结构变形的模拟较成熟，但对于冰裂纹的扩展模拟不准确。基于此，有研究者将内聚力单元法（cohesive element method, CEM）应用于冰裂纹的模拟。CEM 法是数值模型中模拟裂纹萌生和扩展的专用有限元法，最早可以追溯到研究者对裂缝尖端Dugdale-Barenblatt（D-B）[4-5]模型的开发和研究。关于内聚力模型的最初思想是由Barenblatt[5]提出的原子牵引分离定律，应用该定律可以避免裂纹尖端处不切实际的应力奇异性。后来，Dugdale[4]使用内聚力模型来描述受到法向应力影响的理想弹塑性材料在裂纹尖端附近的塑性变形。而后，Needleman[6]将D-B 模型中的裂纹区域称为内聚力区（cohesive zone）。内聚力模型在许多研究者的研究中得以改进并应用于多个领域，与其他裂纹计算模型相比，内聚力模型具有能够模拟多个裂纹路径而无需采用裂纹路径跟踪的优点。此外，CEM 法无需预先知道裂纹扩展的方向，可通过任意放置内聚力单元来模拟裂纹路径的传播。本文将主要研究CEM 法在冰裂纹扩展方面的应用，这对于冰区结构物设计制造有着重要意义。

内聚力模型凭借在裂纹扩展和材料断裂方面的优势，近年来，被越来越多地应用于冰力学和结构物与冰碰撞的模拟中。例如，Mulmule[7]和Dempsey[8]采用I 型断裂模式的冰材料对内聚定律进行了反算；Kuutti 等[9]通过将二维内聚力单元插入三角形冰单元周围，模拟了垂直钢板与厚冰块碰撞的试验，结果发现密集的有限元模型网格可更好地模拟试验中冰块的破坏过程；Lu 等[10]和Wang 等[11]均采用CEM 法模拟了锥形体与层冰的碰撞过程；王峰等[12]、刘路平等[13]和蒋昱妍等[14]则采用CEM 法模拟了平整冰与竖直圆柱体、海洋平台等海洋结构物的碰撞过程，验证了CEM 法在处理此类碰撞问题时的可靠性。

目前，尽管CEM 法在研究结构物−冰相互作用方面已得到一定的应用，但关于冰的内聚区长度公式的研究仍十分缺乏。有部分研究者直接采用了其他材料的公式，例如Lu 等[10]采用弹性材料的公式结合冰裂纹试验，初步估算了内聚区长度。然而，大部分研究者并未关注内聚区长度的计算问题，例如Gürtner 等[15]、Kuutti 等[9]、Wang 等[11]在处理各种结构物−冰碰撞问题中均未提及内聚区长度与网格的关系，而在划分网格时多采用试错的方法。事实上，基于其他材料的研究成果在冰这种材料上不一定适用，且以往研究内聚区长度主要集中在低厚度模型[16-17]，例如复合材料的夹层，而对高厚度模型的研究很少。对于冰的断裂模式而言，冰通常都是高厚度模型，现有的适用于薄板类型的内聚区长度不再适用。

对于上述研究中存在的不足，本文拟针对冰材料的内聚区长度计算和网格划分问题进行新的适用性研究。首先，重现及推导现有的内聚区长度计算公式，并对公式进行修正；然后，基于有限元法建立双悬臂梁数值模型进行模拟，以研究单个内聚区长度内的有限元模型网格密度；最后，基于LS-DYNA 软件将相关结果应用于冰的三点弯曲试验，以验证CEM 法应用于模拟冰的I 型裂纹萌生和扩展时的有效性，期望能够促进CEM法在冰材料力学特性研究领域中的应用。

1 内聚区长度计算公式及原理

CEM 法在提出后主要应用于混凝土结构和复合材料等领域，展现出了良好的裂纹模拟能力，可以形象地模拟材料的张开、滑开、撕开等开裂破坏过程，以及很好地模拟混凝土的碎裂和复合材料的分层失效过程。目前，控制内聚区的常见牵引力−分离关系包括4 种形式：双线性型[18]、多项式型[19]、梯型[20]及指数型[21]，如图1 所示。图中，σ0为裂纹形成的最大应力， δn为损伤起始对应的上下面位移值， δf为“梯形”应力开始线性下降对应的上下面位移值， δ0为裂纹形成时上下面的最大位移值。

图1 牵引力−分离曲线Fig. 1 Traction-separation curve

图2 所示为不同材料的内聚力模型方向示意图。图中，定义裂纹扩展方向为长度方向、裂纹向外张开方向为厚度方向。由两图的对比可见，海冰属于典型的高厚度模型。

图2 不同材料内聚力模型方向定义Fig. 2 Defined dircections of cohesion model of different materials

内聚区长度可通过J积分推导计算。J积分是弹塑性断裂力学中一个与路径无关的积分[22]，可作为裂纹或缺口顶端应变场的平均度量。本文应用张开破坏的双悬臂梁模型来简述裂纹扩展区域长度的推导过程，如图3 所示。图中，Lcz为裂纹扩展区长度（也即内聚区长度），L为无初始裂纹的双悬臂梁长度，Lpc为初始裂纹长度，L0为整个双悬臂梁长度，h为单个悬臂梁的厚度。双悬臂梁左端的初始裂纹用于约束裂纹的产生位置和扩展方向，外扭矩M施加于双悬臂梁左端。

图3 双悬臂梁撕裂示意图Fig. 3 Schematics of tearing of double cantilever beam

J积分的回路积分定义公式如下[22]：

为将所有常系数（含无量纲系数a的项及常数项）统一到一起并用字符代替，使公式具有一般性，式（3）可进一步变形为

式中：c为关于L和h的无量纲系数。将式（4）进一步简化为

其中，

Hillerborg 等[23]基于D-B 模型得出原始的裂纹扩展区长度与裂纹位移和应力间的关系为

式中，f为裂缝形成时的最大应力，与式（6）中的σ0等价，即f=σ0；G=kδ0f，为断裂能释放率，其中k为常数，G适用于I 和II 型断裂模式，可以理解为图1 中各类型曲线与坐标系包围的面积，该面积又可表达为最大应力、最大应变与一个常数的乘积。对比式（6）和式（7），可知文献[23]的与式（6）中的c2是一致的，同样也适用于撕裂（内聚力单元的I 型拉伸失效）破坏。为此，式（5）可改写为

上述相关推导起初源于复合材料的裂纹扩展和计算裂缝区域长度，但因推导过程中并未限制材料的类型及属性，所以上述描述裂纹的思想也适用于冰的I 型裂纹的萌生和扩展。文献[23]提及的双线性δ-σ 曲线适用于混凝土这类脆性或准脆性材料，本文应用CEM 法，在冰的力学模拟中也采用了双线性δ-σ 曲线。

应用式（8）时，系数c1的选取是一大难点。以往相关研究者主要针对薄板模型，根据材料性质、内聚力模型本构关系等，对c1的取值给出了不同的推荐值，例如0.5[24]，π/8，0.731 或2.92[25]等。但是，这些参数均采用的是简单的固定值，忽略了c1本身也会随着模型尺寸变化的特征，且它们只对于低厚度的薄板材料（图2（a））较为适用，而对于高厚度冰材料（图2（b））的适用性尚未见诸于相关文献报道。为此，本文将针对系数c1在冰材料中应用的取值进行研究。

2 应用于冰材料的内聚区长度

本文在给出系数c1的修正公式之前，首先探讨内聚区长度Lcz随模型尺寸的变化规律。这里，将关注两个值：厚度值h和长厚比L/h。在有限元模型中，可以通过单元的应力或应变是否超过模型中的设定值来判断单元是否失效[26]。本文内聚力单元的变形位移量是通过追踪每个时刻单元节点坐标计算得到的。基于位移法则，通过节点坐标的变化来判断内聚力单元是否处于不可恢复的变形状态，结合牵引力−分离曲线（图1），来判定内聚力单元是否失效。为降低误差，本文尽可能减小了坐标数据输出的时间间隔，将输出间隔时间设置为5×10−5s。

构建双悬臂梁数值模型。双悬臂整个长度L0=150 mm，单个悬臂梁厚度h=1.55 mm，模型受力端有一个初始裂纹长度Lpc=35 mm 的切口，无初始裂纹的双悬臂梁长度L=115 mm。将裂纹扩展长度设为内聚区长度Lcz，双悬臂梁受方向相反、大小相等的拉力F作用。模型中的内聚力单元受力变形如图4 所示。由图可见，当加载时间由0 s 增加至0.35 s 时，内聚力单元随着力的加载逐渐被拉伸变形；当加载时间达到0.5 s 时，内聚力单元因受力达到极限并失效。

图4 内聚力单元拉伸变形Fig. 4 Tensile deformation of cohesive element

在保证初始裂纹长度Lpc=35 mm 不变，而改变无初始裂纹的双悬臂梁长度L（115，100 和85 mm）的情况下，绘制了如图5 所示双悬臂梁破坏过程中内聚区长度Lcz随单个悬臂梁厚度h的变化曲线，并进一步绘制出如图6 所示Lcz随Ln(L/h)的变化曲线。

图5 Lcz -h 曲线Fig. 5 Lcz-h curve

图6 Lcz-Ln(L/h)曲线Fig. 6 Lcz-Ln((L/h) curve

由图5 可见，整体内聚力区长度Lcz会随无初始裂纹双悬臂梁长度L的增加而逐渐增加。在增加单个悬臂梁厚度h的情况下，Lcz相对于h增加的变化率逐渐减小。当h增加到一定数值后，Lcz基本不再发生变化，但对于低厚度条件下的Lcz变化影响较大。

由图6 可见，长厚比L/h的增加对Lcz的影响也呈现很明显的下降趋势。当L/h>30（Ln(30)≈3.40）时，Lcz基本不再发生变化；当L/h减小时，Lcz受L/h变化的影响逐渐敏感；当L/h进一步减小，达到L/h<1.25（Ln（1.25）≈0.223）时，Lcz再次在某个稳定值附近波动，本文将出现该现象的点称为拐点。根据3 种不同长度L的无初始裂纹双悬臂梁模型，拐点出现时，L/h的比值均在1.25 附近。由图5 和图6 可见，当L/h在1.25～30 的范围内取值时，Lcz受到L/h值变化的影响比较明显；当L/h不在该范围时，其对Lcz的影响比较小。

3 修正内聚区长度计算公式

目前，CEM 法主要应用于较薄，即L/h值偏大[16]的复合材料的力学特性研究中，其内聚区长度公式的修正系数以常数项修正为主，未考虑长厚比L/h值由1.25 趋近1 时内聚区长度不再随厚度变化的现象。而在实际的模拟中，应考虑存在内聚区长度变化的拐点，即在L/h到达一定数值后，内聚区长度变化不再明显的现象。本文拟在前人研究的基础上对系数c1的选取进行修正，将其与L/h值的变化联系起来，以提高内聚区长度公式的精确性。

根据不同长度L的无初始裂纹双悬臂梁模型，本文选取L/h值大于1.25 的部分，并假设c1为关于L与h的函数，再根据式（8）对模拟中的Lcz–h散点进行多次数据拟合，以建立c1与L/h的关系式。通过对数据的拟合整理，可得：

代入式（8），可得到新的内聚区长度计算公式为

上式中：x为当前L/h的取值；hmin为模拟的最小模型厚度。当hmin远小于模型长度时，hmin对计算结果影响很小，对于冰材料目前尚未见相关计算，本文参照文献[5] 和[16] 中对于复合材料的取值，选择hmin=1.5 mm。对于冰的I 型裂纹的萌生和扩展， σ0取为冰的抗拉强度，其值的选取参考了冰的抗拉试验[27]。此外，根据上节所述，拐点出现后内聚区长度基本不变，故该公式仅适用于L/h大于1.25 的情况，对于L/h小于1.25 的情况，均取L/h=1.25 对应的内聚区长度。

下文将验证系数c1对内聚区计算结果修正的有效性。重新构建一个高厚度的双悬臂梁数值模型，如图7 所示，其宽度B为20 mm，长度L分别为70 和210 mm，Lpc仍为35 mm，厚度2h不断变化。在1～150 的L/h范围内适当选取一定的散点绘制相应曲线，将模型计算结果与公式修正后的结果进行比较，如图8 所示。

图7 高厚度的双悬臂梁模型Fig. 7 Model of high-thickness double cantilever beam

由图8 可见，在2 种不同尺寸模型下，修正模型模拟结果与直接模拟结果在L/h的所有取值范围内都有较好一致性。由分析可知，除验证公式的有效性外，拐点处（L/h=1.25，即图8（a）中的虚线位置）的内聚力区长度与双臂梁模型整个长度也有一定的关系。随着整个长度和内聚区长度的增加，公式的计算结果与模拟结果间的误差逐渐变小，这点在L/h取值较小的范围内尤其明显。

图8 模拟结果和修正内聚区长度对比Fig. 8 Comparison of simulation results and modified cohesive zone length

4 修正公式在冰材料上的应用

网格密度决定了计算量大小。在内聚区长度Lcz中内聚力单元的密度越大，裂纹形成和力的曲线越光滑，精确性越高，反之，力的曲线更容易出现锯齿形波动。在单个内聚区范围内设置合适的内聚力单元个数可以很好地模拟裂纹形成和力的传递过程。基于上文Lcz计算方法，本节以内聚力单元的I 型拉伸失效为研究对象，研究一个Lcz中应布置的内聚力单元个数，即Lcz=mSc（其中，Sc为内聚力单元尺寸，m为单元数）。

本节狭长双悬臂梁模型的主尺度（L×B×2h）为150 mm×20 mm×3.1 mm，如图9 所示。冰参数如下：弹性模量5.72 GPa，剪切模量2.2 GPa，硬化模量4.26 GPa，体积模量5.26 GPa，上下牵引的运动速度为0.004 m/s。网格尺寸依次从2.0 mm 细化到0.5 mm，间隔0.05 mm。模拟在拉力带动下狭长双悬臂梁间的内聚力单元的I 型拉伸失效过程，并测试不同网格密度下网格密度对失效拉力和位移曲线的影响，如图10 所示。

图9 狭长双悬臂梁模型主尺度示意图Fig. 9 Schematics of main dimensions of narrow double cantilever beam model

由图10 可见，随着网格的细化，内聚力单元需传递的牵引力−时间曲线逐渐趋于光滑，整体波动趋于缓和，数值逐渐向中间靠拢。根据式（10）计算的内聚区长度Lcz=2 mm，结合上述模型数据，实际模拟得到的内聚区长度也在2 mm 附近。

根据图10 所示不同网格密度下牵引力−时间曲线分布情况，当网格尺寸小于等于0.5 mm 时，曲线的波动基本上可以接受，且0.5 mm 的网格密度所需计算量也可以接受。综上所述，内聚力单元在Ⅰ型拉伸失效模式下一个内聚区长度中至少存在4 个内聚力单元才能较为精确地描述材料的断裂过程，并可较好地保证载荷曲线的精确性和稳定性。

图10 不同网格密度下内聚力单元牵引力−时间曲线Fig. 10 Traction-time curves of cohesive element with different grid densities

在理论上推断了CEM 法网格长度与模型尺寸关系的精确性并进行数值验证后，本节将采用CEM 法对冰的三点弯曲试验进行数值模拟，以呈现冰由变形到破坏的过程。数值模型[28]主尺度（L×B×2h）为70 mm×70 mm×650 mm，支点间距离为600 mm，如图11 所示。相关参数如表1 所示。数值模型按照实际的试验模型建立，冰参数的取值根据试验结果推导得出。

表1 试验参数Table 1 Experimental parameters

图11 冰试样模型主尺度示意图Fig. 11 Schematics of main dimensions of ice specimen

为保证良好的网格稳定性和计算结果的精确性，根据式（10）计算所需内聚区长度Lcz，按照裂缝发展方向，以及上文中对于内聚区长度公式中厚度方向的定义，确定了三点弯曲试验模型的长厚比L/h约为0.215，这明显小于拐点出现时的L/h值1.25，故取拐点处的Lcz值。经计算，得到的拐点处Lcz=15.4 mm，并根据m=4，得到最大内聚力单元网格尺寸为3.85 mm。

为保证网格划分的取整性，计算中裂纹发展方向的内聚力单元长度取值为3.5 mm，垂直于裂纹发展方向的内聚力单元长度为4.5 mm，这样可以节省一定的计算量。在上述网格密度条件下，计算步长为3.68×10−7s，计算时长约4 h。

按照上述试验数据和冰参数构建有限元模型，在可能断裂的试样中间部分添加内聚力单元，添加方式如图12 所示。图中，红色为冰单元，灰白色为内聚力单元。冰力学模型仅提供变形和力的传递，冰的破坏和分离采用内聚力单元代替，故不考虑冰单元的破坏，而受破坏的仅有内聚力单元。冰力学模型为各向同性的线弹塑性模型，数值模拟过程如图13 所示，试验冰试样断裂形式如图14 所示。

图12 冰试样截面间内聚力单元Fig. 12 Cohesive element between ice sample sections

图13 有限元模拟三点弯曲过程Fig. 13 Finite element simulations of three-point bending process

图14 试验冰试样断裂[28]Fig. 14 Fracture of experimental ice specimen[28]

本节数值模拟所用冰模型是均匀的，内聚力单元牵引力法则选用的是双线性型，且在冰力学模型建立过程中不考虑冰的原始缺陷，因此，图13所示冰试样断裂面较光顺，导致整体模拟过程中冰的破坏较为平滑。根据数值模拟中物体的挠度与压力变化，绘制出如图15 所示压力−挠度曲线。

图15 压力−挠度曲线Fig. 15 Pressure-deflection curves

试验中，冰试样被破坏前的最大受力值为1 127.9 N，基于CEM 法模拟三点弯曲的极限载荷为1 161 N，断裂点的误差为2.9%，与试验结果较为接近。从图15 所示压力−挠度曲线的不同阶段来看，CEM 法模拟的曲线趋势和试验曲线基本一致。从最终断裂挠度来看，CEM 法对应的极限挠度值相较于试验数据有一定的误差，但较小。此外，由于该数值模拟选择的冰力学模型为线性模型，且内聚力单元建立过程中未考虑原始缺陷，使得模拟得到的曲线相较于试验数据得到的曲线更光滑。从两个曲线的对比可知，数值模拟过程中冰试样的弯曲变形、破坏过程与试验的基本一致，且峰值过后的载荷卸载过程很好地体现了冰的脆性和裂纹成型的连贯性。总之，本次模拟数据与试验数据的对比结果很好地验证了CEM 法及相关网格划分方式在模拟冰的力学特性方面的有效性。

5 结 语

本文从J积分出发，基于部分假设并结合前人的相关研究结果，重现了内聚区长度计算公式的推导过程，在原有的内聚区长度公式基础上，增加了关于L/h值的修正函数。通过建立新的高厚度双悬臂梁数值模型检验了内聚区长度公式的精确性，解决了该公式是否适用于冰力学模型的问题。通过构建双悬臂梁有限元模型进行数值模拟，结果发现，在一个内聚区长度中至少存在4 个内聚力单元才能够较为精确地描述材料的断裂过程。最后，将研究结果应用于三点弯曲试验模型的模拟中，结果显示断裂点的极限载荷误差为2.9%且在合理范围内，由此验证了修正后的内聚区长度公式在模拟冰的I 型裂纹萌生和扩展时的有效性。