基于扩展有限元法的重组竹Ⅰ型断裂数值模拟分析

刘明凯，周爱萍，刘燕燕，盛宝璐

(南京林业大学生物质材料国家地方联合工程研究中心，南京 210037)

重组竹是竹纤维束通过顺纹组胚、胶合热压制成的材料，其力学性能优于常用结构木材。但由于竹束疏解工艺简单，难以均匀浸胶，重组竹不可避免会带裂缝工作。当复合材料存在基体裂纹、孔洞和缺口时，层间断裂是最常见的失效形式。这种层间分层会导致纤维断裂，降低复合材料的使用寿命。目前，重组竹处于基础性研究阶段，由于完整标准体系和设计计算理论的欠缺，该材料还未大规模投入实际工程应用中。重组竹构件在工程应用时，体内裂纹一般以复合型裂纹的形式出现，应力状态复杂。单一的Ⅰ/Ⅱ型裂纹是研究复合型裂纹的前提，其中，Ⅰ型层间断裂一直是研究人员最为关注的问题，因为与剪切模式相比，拉伸分层所需的能量较低，裂纹更易在该模式下萌生。为防止断裂失效的发生，有必要对复合材料Ⅰ型裂纹的启裂和扩展规律进行研究，这对工程设计和应用指导具有重要意义。

除解析法和试验法是研究断裂问题的重要方法，数值模拟法也是研究裂纹扩展问题的有效手段。随着计算机技术的迅猛发展，运算能力极大提高，断裂力学在有限元软件方面的应用更加广泛。传统有限元法存在局限性，预制裂纹面要与划分的网格重合，并且均要预先确定裂纹扩展的方向。扩展有限元法(extended finite element method，XFEM)克服了上述缺陷，给裂纹扩展在实际应用中带来了便利。虚拟裂纹闭合技术(VCCT)和内聚区模型(CZM)是模拟裂纹扩展的2种常用方法。VCCT基于线弹性理论对裂缝进行分析，忽略了塑性区能量耗散内聚力。CZM的软化阶段有效弥补了VCCT中缺少的塑性效应，该方法已被证明是模拟双材料界面[1-2]和层压复合材料[3-4]断裂的有力工具。Mo⊇s等[5]最初提出了XFEM与CZM相结合的方法，得到了与网格无关的裂纹表示和基于CZM的裂纹扩展。Heidari-Rarani等[6]将XFEM与VCCT、CZM相结合对复合材料进行模拟对比，与XFEM-VCCT相比，XFEM-CZM提供了更准确的结果，可以很好地预测裂纹的萌生和扩展。

目前，有关重组竹的断裂研究主要集中在断裂破坏机理分析和断裂韧度测定，而重组竹的断裂仿真模拟以及断裂韧性与试件尺寸相关的研究还很少。笔者以ABAQUS有限元软件为平台，基于扩展有限元法对不同尺寸的重组竹双悬臂试件进行裂纹扩展数值模拟，将模拟结果和试验结果进行对比来验证扩展有限元法的准确性，该内聚模型可用于建立构件计算模型和确定构件承载力；再通过荷载-位移曲线得到裂纹扩展阻力曲线(R曲线)，并依此进行重组竹断裂韧度的尺寸效应分析。

1 扩展有限元法和内聚区模型

扩展有限元法和常规有限元法的区别在于，扩展有限元法在单位分解法的基础上增加了阶跃函数和渐进场函数来对位移进行具体描述，实现了裂纹和网格独立存在。因此，在模拟裂纹扩展时无须重新划分网格，能方便地分析不连续问题，特别是不连续边界演化问题。求解域上的近似位移场[uh(x)]为：

(1)

在黏性裂纹中，裂纹的扩展受尖端附近裂纹面上牵引位移关系的控制。该模型是由Dugdale[7]和Barenblatt[8]引入金属领域的，后来演化为内聚区模型，能够消除裂纹尖端的应力奇异性，可将裂纹尖端处的应力幅度限制在有物理意义的水平。

图1 内聚区模型Fig. 1 Cohesive zone model

2 试验概述和仿真模拟

2.1 试验概述

断裂韧性测定试验对试件的形状和尺寸有严格的要求，双悬臂梁(DCB)是测试复合材料断裂性能最常用的结构形式。梁高1/2处的预制裂纹将整梁划分为上下两部分，通过拉伸试验研究材料的界面强度。Ⅰ型层间韧性能量释放率试验基于ASTM D5528-13“Standard test method for Mode I interlaminar fracture toughness of unidirectional fiber-reinforced polymer matrix composites” 测定，试件尺寸如图2所示。

a0为初始裂纹长度；a为裂纹长度；l为试件长度；b为试件厚度；h为试件1/2高度；P为荷载。图2 DCB测试试件尺寸Fig. 2 DCB test specimen dimension

分析不同厚度、不同初始裂纹长(同向)的几何相似试样，确定适当的厚度和初始裂纹长，从而评估有效的断裂能，试件具体尺寸见表1。其中，厚度的取值是满足平面应变的条件之一；裂纹长度的取值是满足小范围屈服(塑性区尺寸远小于裂纹尺寸)的条件，只有小范围屈服，才能按线弹性计算出满足工程精度要求的双悬臂梁能量释放率(GIC)。本试验分为A、B两组，其中，A组试件初始裂纹长度相同，厚度不同；B组试件厚度相同，初始裂纹长度不同。试验共7组，其中每组至少包含3个试件。准静态DCB试验在万能试验机上进行，采用0.5 mm/min 恒定速率加载，低速率可以避免裂纹的不稳定传播。在裂纹扩展过程中，载荷-位移(P-δ)曲线是唯一需要记录的数据，是获得R曲线的基础。最后，将R曲线的GIC与模拟得到的GIC进行比较，评估该方法的可行性。

表1 试件尺寸Table 1 Specimen sizes mm

2.2 仿真模拟

模型是仿真的基础，模型的还原度直接决定了分析结果的可靠度。基于有限元软件ABAQUS，建立带预制裂纹的双悬臂模型。在对模型进行荷载设定时，为消除耦合带来的局部扭转效应，建立刚体并将其贯穿于双悬臂预制的销孔内，设置2个刚体的相反位移。

网格划分时，在保证计算精度的前提下为减少相对所需的计算时间，采用网格局部加密。划分部件并在裂纹扩展区附近细化局部网格，得到精确的模拟结果(高度方向每1 mm划分1次)；而在远离裂纹的区域，由于对裂纹尖端附近的应力状态影响较小，使用相对较大的网格(高度方向每10 mm划分1次)，并设置顺纹方向为坐标轴的1方向，其余2个方向为2和3方向，如图3所示。进一步的网格加密研究表明，更细的网格划分对结果精度影响较小。

图3 有限元模型网格划分Fig. 3 Mesh generation of finite element model

在材料模块中定义重组竹的材料特性(E为弹性模量，E1=11 212 MPa、E2=E3=2 561 MPa；ν为泊松比，ν12=ν13=0.304、ν23=0.054；G为剪切模量，G12=G13=1 418 MPa、G23=749 MPa)。重组竹为纤维定向增强复合材料，顺纹与横纹的力学性能差异较大，顺纹方向(1方向)的强度远高于另外2个方向(2和3方向)。实体部件使用C3D8R单元，采用静态分析并打开几何大变形。在XFEM裂纹模块中，定义裂纹面所处位置并允许裂纹沿顺纹方向扩展；“场输出”勾选PHILSM、PSILSM可观察裂纹的扩展状态，勾选STATUSXFEM可在结果中查看扩展有限单元的失效状态。

损伤起始准则采用最大名义应力准则，当分层过程中最大接触应力比值达到1时，材料开始发生损伤，该准则可表示为：

(2)

损伤演化准则为基于位移的、线性软化的损伤演化规律，一旦满足损伤起始准则，材料将根据损伤演化规律发生损伤，该规律描述了内聚区模型刚度退化的速率。令内聚刚度损伤因子为D，在起始阶段D=0，表示材料没有发生损伤，仍处于弹性状态。损伤开始后，D逐渐增加，直到D=1时材料完全失效。线性软化的内聚刚度损伤因子可表示为：

(3)

3 模拟结果及分析

3.1 结果对比

模拟结果的应力分布如图4所示。随着加载点相对位移的增大，裂纹沿着纤维方向向后扩展，与试验观察到的现象一致。通过对模拟数据的后处理获得P-δ曲线，进而获得R曲线，将模拟曲线和试验曲线进行对比。

图4 双悬臂梁仿真模拟应力图Fig. 4 Simulated stress diagram of double cantilever beam

3.1.1 荷载-位移曲线对比与分析

数值模拟得到的P-δ曲线与试验曲线吻合良好，如图5和6所示。峰值荷载的误差也在容许范围内，具体误差值见表2。

对于弹性阶段，由A组曲线对比可知，试件厚度较大时，模拟曲线中弹性阶段的直线斜率和试验斜率吻合度较好。其中，b=100 mm的试件(B-1)模拟曲线中的直线段斜率和峰值荷载均与试验曲线吻合较好。相反，较薄的试件承载能力模拟值略低于试验值，模拟曲线弹性阶段斜率也小于试验曲线的斜率。

最薄试件荷载差异(绝对值)为12.39%，其他厚度试件模拟差异(绝对值)为1.37%～8.40%。该结果与Manshadi等[13]对双悬臂梁模拟的结果一致，认为该差异是由于薄试件缺乏大规模纤维桥连而导致的。模拟曲线的下降段均与试验曲线较吻合，体现出ABAQUS中的损伤演化准则有效地模拟了材料软化阶段。

B组曲线进一步验证了厚试件模拟结果的准确性，3组曲线弹性阶段均吻合，不同的初始裂纹长度对初始刚度的影响较小。初始裂纹越长，双悬臂梁刚度越小，能承受的荷载越小。

综上所述，ABAQUS在对重组竹仿真模拟时，薄试件DCB模拟结果的准确度不及厚试件，双悬臂梁厚度越大，模拟结果越吻合。在实际应用中，为节省财力和物力，厚度大于某一值时，任意尺寸双悬臂梁可以通过该仿真模拟方法得到其可靠刚度。

图5 A组试件荷载-位移曲线对比Fig. 5 Comparison of load-displacement curves of specimens in Group A

图6 B组试件荷载-位移曲线对比Fig. 6 Comparison of load-displacement curves of specimens in Group B

表2 双悬臂梁的峰值荷载误差Table 2 Peak load errors of double cantilever beams

3.1.2 模型准确性验证

以裂纹扩展量(Δa)为函数的裂纹扩展阻力曲线称为R曲线，可由荷载-位移曲线得到。R曲线是一种几何相关性质[14-15]，即初始裂纹长度、厚度、高度、曲率等几何参数均可能影响R曲线。本研究对初始裂纹长度和厚度2种参数进行分析。

对于含有穿透裂纹的DCB试件，本研究采用柔度法计算应变能释放率。在线弹性断裂力学中，柔度法是测定材料断裂韧度的方法之一。

由经典梁理论，加载点的位移δ，即双悬臂梁的挠度可表示为：

(4)

式中：ω为单悬臂的挠度；E为弹性模量，此处使用顺纹方向的弹性模量E1。则柔度C可表示为：

(5)

考虑DCB为线弹性的情况下，Irwin[16]给出了能量释放率：

(6)

式中，柔度的立方根(C1/3)作为a的函数。

在实际应用中，由于双悬臂梁端并非处于完全固定边界条件下，且裂纹前缘可能因旋转偏离中间平面，该方法计算得到的G偏大，因此，须对G进行修正。纠正裂纹尖端旋转的一种方法是引入校正参数Δ，即双悬臂梁的有效裂纹长度为a+Δ。对于半厚度的DCB试样，Williams[17]提出了Δ的计算公式：

(7)

以a+Δ代替式(6)中的a得到修正后的G。A组和B组试件R曲线对比分别如图7和8所示。由图7和8可知，修正后的R曲线结果较准确，a0符合实际尺寸，且模拟曲线的水平段和试验曲线相近。

图7 A组试件R曲线对比Fig. 7 Comparison of R curves of specimens in Group A

图8 B组试件R曲线对比Fig. 8 Comparison of R curves of specimens in Group B

R曲线的试验结果与模拟结果在趋势上存在差异，试验结果在转折后呈上升趋势，而模拟结果在转折后呈下降趋势。该偏差可从两方面进行分析：一方面，仿真模拟时采用理想化模型的裂纹面为理想平面，而实际复合材料的非均质特性导致裂纹面呈现曲面，纤维桥连现象的增多使得断裂所需的能量相应增多，因此试验曲线呈上升趋势，仿真模型和实际试验存在差异；另一方面，牵引参数通常被认为是材料参数，文献[12]中的内聚参数是基于一定尺寸构件通过试验反推得到的，跟尺寸有一定的关联。该文献中采用的构件厚度为80和100 mm，初始裂纹长度为100和134 mm，本研究使用该内聚参数得到的模拟结果中，厚度为100 mm且初始裂纹长度为102和136 mm的DCB模拟得到的R曲线吻合度较高，R曲线趋于水平。对于其他尺寸的试件，模拟结果趋势存在一定偏差。

3.2 尺寸效应

3.2.1 厚 度

汇总不同厚度试件的临界能量释放率，对比试验值和模拟值，如图9所示。当b≥40 mm时，临界能量释放率均保持在2.0 J/m2左右，即该厚度范围均可看成平面应变问题；但当厚度较小且b=20 mm时，试验测得的GIC高于模拟得到的GIC，较薄的试样具有较大的断裂韧度。材料断裂的临界能量释放率随厚度的增加而逐渐减小，最终趋于一个恒定的较低值。这一结果验证了断裂韧度随试样厚度的变化关系，如图10所示。

图9 不同厚度试件的临界能量释放率Fig. 9 Critical energy release rate of specimens with different thicknesses

图10 断裂韧度与厚度关系图Fig. 10 Relationship between fracture toughness and thickness

Hu等[18]认为与尺寸效应相关的准脆性断裂行为实际上是由于裂纹尖端处的FPZ与最近的结构边界相互作用的结果，断裂过程区(如桥连区)到最近边界的距离控制着断裂行为。FPZ在小试件中相对较大，因此，断裂以强度准则为主。在大试样中，FPZ相对于试样尺寸及其到最近试样边界的距离较小(FPZ的绝对尺寸没有意义)，因此，断裂以断裂韧性准则为主。

实际上，当构件厚度与塑性区尺寸的比值较小时，厚度方向的约束较弱，材料的塑性变形在该方向上不受限制，厚度方向可以自由地发生屈服，易产生45°斜断口；当构件厚度与塑性区尺寸的比值较大时，厚度方向不能自由产生屈服，周围弹性材料的限制使厚度方向应变为0，此时裂纹呈平断口。对于较厚的构件，两自由表面在厚度方向的应力为0而处于平面应力状态，但构件内部绝大部分处于平面应变状态，如图11所示。厚度越大，平面应力区域所占比例越小，因而可将构件整体视为完全处于平面应变区域。有限元分析结果中，双悬臂梁的剖面应力图验证了上述关于厚度的尺寸效应，如图12所示。该结果和陈涛等[19]研究铝合金材料R曲线的三维效应类似，薄试件有较大的断裂韧性，随着厚度的增加，阻力曲线趋于一条水平直线。

图11 裂纹尖端三维塑性区Fig. 11 3D plastic zone at the crack tip

图12 双悬臂梁模拟剖面应力图Fig. 12 Stress diagram of simulated section of double cantilever beam

3.2.2 初始裂纹长度

ASTM D5528-13建议裂纹分层前沿应距加载线50 mm，而ESIS TC4 Protocol “Determination of the mode I delamination resistance of unidirectional fibre-reinforced polymer laminates using the double cantilever beam specimen”建议从加载块(或钢铰链)的前边缘到薄膜前沿的距离应至少为45 mm。

随着初始裂纹的增长，裂纹启裂所需要的能量减小。初始裂纹较短时，横向剪切变形占主导地位，裂纹扩展所需的能量较大；初始裂纹较长时，双悬臂梁会产生大挠度变形，裂纹前端接近DCB试件的末端，临界能量释放率较小。此时为消除端部效应，应制备更长的试件。不同初始裂纹长度试件的临界能量释放率如图13所示，当试件处于平面应变条件时，初始裂纹长度对模拟结果的吻合度影响较小。

图13 不同初始裂纹长度试件的临界能量释放率Fig. 13 Critical energy release rate of specimens with different initial crack lengths

4 结 论

本研究运用ABAQUS软件对重组竹的裂纹扩展进行了三维仿真模拟，得到了和试验吻合度较好的模拟结果，主要结论如下：

1)本研究得到了基于内聚区模型的扩展有限元分析模型，在可靠内聚参数设定下，模拟裂纹扩展过程，得到较吻合的荷载-位移曲线和裂纹扩展阻力曲线。

2)试件厚度较大时，P-δ曲线吻合度较好；修正后的R曲线较准确，双悬臂梁的能量释放率随厚度的增加先下降后趋于平稳，并且当b≥40 mm时趋于稳定。因此，厚度达到40 mm后，模拟结果和试验结果吻合。

3)初始裂纹长度越大，双悬臂梁承载能力越低，能量释放率越小，且初始裂纹长度对模拟的吻合度影响较小。