环形缺口小冲杆试样结合内聚力模型提取断裂韧性参数

张 宇， 刘海亭， 翁 琳， 沈 耀

(上海交通大学 材料科学与工程学院, 上海 201100)

在高温高压或者辐照的恶劣环境下，材料长期服役会导致其力学性能退化，容易引发事故[1].对于核电设备材料而言，辐照后组织和性能的均匀性显著降低，如果韧性不足则会造成断裂[2].所以从安全以及经济效益方面考虑，评价在役设备和辐照材料的断裂韧性至关重要.常见的断裂韧性指标为临界应力强度因子KIC或断裂韧度JIC，可通过常规力学性能试验获得.然而，某些在役设备和辐照材料取样条件有限，无法满足常规试验样品的尺寸要求[3]，因此需要采用小样品外推获取.小冲杆测试作为小样品技术的一种，最早由Baik等[4]提出，它有着近乎无损取样、耗材少等优点[5]，可以在任何承压设备的表面取样，且无需对设备进行修补，弥补了常规力学性能试验因材料不足而无法进行的缺陷.起初，小冲杆测试技术用于核反应堆壳体材料辐照后的脆化评定，大幅度节省了辐照样品的取样量.到20世纪90年代，小冲杆测试技术开始应用到各个工程领域，可以用来评估材料的屈服强度、抗拉强度以及断裂韧性等力学性能，这对在役设备的剩余寿命预测有着重要意义.

基于小冲杆试验提取材料断裂韧性的主要方法主要有经验关联法[6]、应变能密度法[7]、人工神经网络法[8]以及反向有限元法[9].经验关联法对一系列小冲杆试验和标准断裂韧性试验的数据进行关联，建立两者之间的关系.但关联公式因人而异，受试验设备、试样类型和尺寸的影响很大.应变能密度法最早由Petersen等[10]提出，他们将应变能密度定义为应力应变曲线下的面积，启裂点处的值对应为临界应变能密度.当裂尖附近最小应变能密度超过临界值时，裂纹发生扩展[11]，但是该方法的启裂点确定比较困难.人工神经网络法需预先构造一个包含很多材料参数及其对应有限元模拟结果的数据库，再利用人工神经网络通过试验结果从数据库中找出相匹配的材料参数.该方法虽然准确度高，但前期准备工作量大，训练过程复杂.反向有限元法基于某种回归算法调整模型材料参数，使模拟曲线逐渐逼近试验结果，当两者之间的差别满足设置的精度要求时，便可获取材料参数.与前3种方法相比，反向有限元法有着更高的精度，且实现过程相对容易.Yang等[12]采用基于黄金分割搜索算法的反向有限元方法提取材料参数，很好地预测了材料的屈服应力和硬化指数.但是想要预测材料的断裂性能，还需对算法进行改进和优化.

采用反向有限元法提取材料的断裂韧性需要引入具体的损伤模型.目前比较常见的损伤模型有GTN模型和内聚力模型.GTN模型是最初由Gurson提出，而后由Tvergaard和Needleman发展和完善的细观韧性断裂模型[13].GTN模型能够显式考虑材料的塑性损伤，包含9个模型参数，其中需要标定的参数有3个，每个参数影响载荷位移曲线的不同阶段[14].内聚力模型将材料分为无损伤连续基体和内聚力单元两部分[15]，能够显式描述材料内部的裂纹扩展行为.相比于GTN模型，内聚力模型中与断裂相关的控制性参数少，只有极限分离应力σc与断裂能Γ0两个参数，有利于提高反向有限元法的准确性.

常规的小冲杆样品在加载过程中不会产生应力集中，适用于提取弹塑性本构关系.若要获取材料的断裂韧性，需要采用预制缺口的小冲杆试样.常见的预制缺口有中心贯穿孔型、直线型和环形.本文采用环形缺口试样，因其最接近于小冲杆试样的轴对称特征[16].为确保反向有限元法的精度，需对样品的几何尺寸和缺口尺寸进行优化，以提升模拟载荷位移曲线对模型参数的敏感程度.

为了从小冲杆测试中准确获取被测材料的断裂韧性，本文基于内聚力损伤模型，对缺口样品的尺寸优化进行了研究，对反向有限元法的算法进行改进.首先，建立环形缺口小冲杆试验的有限元模型，从样品的径厚比、裂纹深度和有无预制裂纹3个方面分析内聚力模型参数的敏感性.然后，采用基于遗传算法的反向有限元法，验证样品优化设计方案的反算准确性.最后，结合随机游走模型对收敛算法进行优化，进一步提高反向有限元法的效率和精度.

1 理论模型

1.1 内聚力模型

内聚力模型最早是由Matos等[17]提出，它们将裂纹问题视为非线性边值问题，适用于大变形[18]过程，可以很好地描述韧性断裂和裂纹扩展过程[19].在有限元中，内聚力模型通过基于牵引分离模型的内聚力单元实现.牵引分离模型有很多种，包括双线性模型、指数模型、常数模型和梯形模型[20]等.本文采用双线性模型，该模型只有σc和极限分离位移δ0两个控制性参数，两者都与样品的几何形状及尺寸无关.考虑到双线性模型中，Γ0(曲线下围成的面积)可以用σc和δ0唯一确定，因此δ0可以和Γ0相互转换.本文选取σc和Γ0作为有限元中内聚力模型参数的输入，采用反向有限元法将其提取后，可用于标准断裂试样的模拟，进而获得材料的断裂韧性.

1.2 有限元模型

本文考虑的样品为直径3 mm的环形缺口小冲杆薄片，厚度d分别为0.1、0.2及0.3 mm.设计的缺口深度d′分别为各自厚度的0.4、0.5及0.6倍，缺口宽度为各自厚度的0.1倍.考虑到试样为环形缺口的薄圆片，采用二维轴对称的方式建模，有限元软件为Abaqus.内聚力模型采用零厚度的内聚力单元实现.以T91钢的J2本构模型作为有限元模型本构参数的输入.

在裂纹尖端附近，应力梯度较大，网格要适当密集，而远离裂纹尖端的部位，网格要适当稀疏，这样才能同时保证计算效率和精度.因此，本文采用过渡网格策略[21]进行网格划分.图1为小冲杆试验的有限元模型，其中红线为0厚度的内聚力单元，单元类型为COHAX4，绿色网格为基体单元，单元类型为CAX4R.上夹具、下夹具和冲杆都设为解析刚体，RP-1、RP-2、RP-3分别为各自刚体的参考点.冲杆的直径为1 mm，上卡具的内径为1.02 mm，下卡具的内径为1.5 mm.接触方式为面对面接触，滑移方式为有限滑移，摩擦因子为0.1.

图1 小冲杆试验二维轴对称模型

1.3 算法模型

普通的回归算法在处理最优化问题时，容易陷入局部最优解，为获取全局最优解，本文采用两种智能算法：遗传算法和随机游走算法.遗传算法在搜索过程中引入了随机因素，能够消除对初始值的依赖性.初始生成的参数组集合称为初始种群，种群中每个参数组被称为个体，每个个体都有自身的适应度，本文将适应度定义为残差的倒数.每一轮迭代根据适应度的大小以几何概率来选择遗传到下一代的个体，再根据概率进行交叉和变异，产生下一代的种群.最终，末代种群产生最优个体，得到全局最优解.随机游走模型需要一个参数组作为起点，通过随机生成的n维单位向量，按照下式进行游走：

(1)

图2 算法流程图

2 内聚力模型的参数敏感性分析

反向有限元法的核心思想为基于某种回归算法调整模型材料参数，使模拟曲线逐渐逼近试验结果.若改变模型参数对载荷位移曲线的影响较小，则反算结果难以收敛，进而导致较大的误差.因此，需要从样品的几何和缺口设计出发，对内聚力模型参数作敏感性分析.本文首先根据未开缺口的小冲杆试样的模拟结果优选环形缺口的位置，再从试样直径与厚度的比例、缺口深度以及有无预制裂纹3个方面来分析内聚力模型参数的敏感性.

2.1 小冲杆试样的缺口位置

对无缺口的试样进行有限元模拟，研究试样变形时最大主应力点的位置以及应力方向，以此优选缺口的位置.

3种厚度d的无缺口小冲杆试样有限元模拟结果如图3所示.在变形的过程中，将最大主应力Smax的位置标记出来，观察其移动轨迹.结果表明：当加载位移到一定值后，3种厚度试样最大主应力的位置都不再变化，始终停留在图中黄色圆点处.同时，从图中可以看出，最大主应力在径向的分量较大，轴向的分量较小，因此断裂以径向拉开为主，符合环形缺口小冲杆试样的断裂规律.经过测量可知，不同厚度试样的最大主应力的位置距离试样中心都是0.3 mm左右，因此在距离试样中心0.3 mm的圆周上开缺口较为合理.

图3 无缺口小冲杆试样最大主应力分布图

2.2 小冲杆试样直径与厚度的比例

内聚力模型的参数需要在合理区间内选取.参数值过高，内聚力模型难以收敛；参数值过低，则载荷位移曲线过早下降，与实际情况不符.经过尝试，对于T91钢，极限分离应力取值约在800～1 200 MPa，断裂能取值约在5～8 MPa·mm时，能够较好满足以上两点.因此，本文采用的5组内聚力模型参数分别为(σc=900 MPa，Γ0=5 MPa·mm)、(σc=800 MPa，Γ0=8 MPa·mm)、(σc=1 100 MPa，Γ0=7 MPa·mm)、(σc=1 000 MPa，Γ0=8 MPa·mm)及(σc=800 MPa，Γ0=5 MPa·mm).

图4为3种厚度试样在不同参数下模拟得到的结果，图中D为加载位移，F为加载反力.可以看到，其他条件相同时，厚度越小的试样，对内聚力模型参数的敏感性越强.对于d=0.3 mm的试样，改变参数对载荷位移曲线的影响不大，5条曲线几乎重合.对于d=0.2 mm的试样，5条曲线在断裂之后的区别较为明显，断裂之前区分度不高，且(σc=1 100 MPa，Γ0=7 MPa·mm)和(σc=1 000 MPa，Γ0=8 MPa·mm)对应的曲线较为接近.因为在断裂过程中网格畸变严重，所以载荷位移曲线后期不太光滑.对于d=0.1 mm的试样，5条曲线在断裂前后的区别都很明显，层次感分明，不同的参数组对应不同的最大载荷及断裂时刻的位移.由此可知，相比之下，d=0.1 mm的试样对内聚力模型的参数最为敏感.

图4 厚度对内聚力模型参数的影响

2.3 小冲杆试样的缺口深度

考虑到0.1 mm厚度的试样对内聚力模型参数最为敏感， 所以使用0.1 mm厚度的试样来分析缺口深度对参数敏感性的影响.缺口深度在厚度一半的附近选取，分别为0.04、0.05和0.06 mm.采用的内聚力模型参数与3.2节相同.

图5为3种缺口深度的试样在不同参数下模拟得到的结果.可以看到，在其他条件相同的情况下，缺口深度越大，对内聚力模型参数的敏感性越强，但是模型的收敛性会变差.对于缺口深度0.04 mm的试样，(σc=1 100 MPa，Γ0=7 MPa·mm)和(σc=1 000 MPa，Γ0=8 MPa·mm)得到的两条曲线几乎一致，说明该样品的参数敏感程度较低.对于缺口深度为0.05 mm和0.06 mm的试样，5条曲线差别都比较明显.然而，后者采用 (σc=800 MPa，Γ0=8 MPa·mm)和(σc=1 100 MPa，Γ0=7 MPa·mm)得到的曲线几乎没有下降的阶段，说明缺口深度为0.06 mm试样的内聚力模型收敛较为困难.基于敏感性和收敛性两个方面考虑， 择优选取的缺口深度为0.05 mm.

图5 厚度0.1 mm试样缺口深度对内聚力模型参数的影响

2.4 有无预制裂纹

由于实际加工的不稳定性，环形缺口小冲杆试样的裂纹尖端可能不够尖锐，应力集中也可能偏离裂尖的正上方.为了能够按照预期的路径发生断裂，可以考虑在缺口的正上方预制裂纹.但是，对于厚度较小的小冲杆试样，在开缺口的基础上预制裂纹难度很大.因此，有必要研究预制裂纹样品的参数敏感性，以确定引入预制裂纹的必要程度.

紧凑拉伸型(CT)试样中，预制裂纹的长度与剩余韧带长度的比例约为1∶10.基于相同比例，小冲杆样品的预制裂纹长度设为5 μm，如图6(a)所示，其中红色部分为内聚力单元.以d=0.1 mm的试样进行模拟，内聚力模型参数与3.3节相同，模拟结果如图6(b)所示. (σc=1 100 MPa，Γ0=7 MPa·mm)和(σc=1 000 MPa，Γ0=8 MPa·mm)得到的曲线重叠，参数组(σc=900 MPa，Γ0=5 MPa·mm ) 和(σc=800 MPa，Γ0=5 MPa·mm)的结果也很接近，载荷位移曲线的区分度不高.与无预制裂纹的结果相比，开预制裂纹后的试样对参数的敏感性有所降低.原因可能是有预制裂纹后，应力集中更明显，断裂更加容易，改变参数对断裂过程的影响变小.因此，预制裂纹对参数的敏感性的提升较小，且实际加工也较为困难，所以从经济效益和反算精度两个方面来看，都无预制裂纹的必要.

图6 厚度0.1 mm试样开预制裂纹的有限元模型及模拟结果

上述3方面的参数敏感性分析表明，通过对小冲杆样品进行几何尺寸的优化能够显著提升模拟结果对内聚力参数的敏感性.优化样品具有如下关键几何特征：厚度为0.1 mm，缺口深度为0.05 mm，无需预制疲劳裂纹.

3 反向有限元法提取模型参数

3.1 反向有限元法的准确性验证

理论上讲，可以基于被测材料的试验载荷位移曲线获取内聚力模型参数，通过对比模型参数的理论值，验证反向有限元法的准确性.然而，每种材料对应的内聚力模型参数具体数值未知，也无标准可查，难以判断反算结果是否准确.此外，某些在役设备和辐照材料的取样条件有限，进一步增加了试验验证的困难程度.因此，基于上述优化样品设计，本文以给定的 (σc=800 MPa，Γ0=5 MPa·mm，记为参数组1)和(σc=1 000 MPa，Γ0=6 MPa·mm，记为参数组2)进行有限元模拟，再将模拟所得载荷位移曲线作为逼近目标，采用反向有限元进行参数提取.通过比较反算参数与给定参数之间的差别，对反向有限元法的准确性进行验证.

残差函数定义为目标曲线和模拟曲线之间的面积差.基于遗传算法的反向有限元计算结果如图7和表1所示.可以看出，经过 1 500 次左右的迭代后，模拟和试验的载荷位移曲线差别已经很小，断裂之前的部分几乎重合，断裂后的部分稍有差异.迭代过程中，残差函数值整体呈下降趋势，最终反算获得的参数与给定参数较为接近，误差在6%以内.可以看出，基于遗传算法的反向有限元方法能够较为准确地提取内聚力模型参数.

图7 遗传算法的反算结果

表1 初始参数值与遗传算法的反算结果

3.2 搜索算法的改进

以上结果若是继续以遗传算法迭代，精度会有所提高，但是耗时过长，计算成本太高.所以将遗传算法和随机游走算法相结合，先将迭代解逼近到全局最优解附近，再采用收敛速度较快的随机游走算法，以达到缩短反算时间并提高精度的目的.采用种群内最多的重复个体数量来判断迭代解是否到达全局最优解附近.遗传算法是根据几何概率进行选择，个体适应度越大，被选择的次数也就越多，遗传到下一代就会出现重复的个体.经过多轮遗传，个体适应度会越来越大，重复个体的数量也随之增加.因此，当种群内最多的重复个体数量超过给定阈值后，便可认为迭代解逼近已到达全局最优解的附近.改进算法后的结果如表2和图8所示.可以看出，改进算法的迭代次数和迭代时间均得到了大幅降低，降幅约为2/3；另一方面，模拟和试验曲线几乎完全重合，反算所得参数与给定参数间的误差减小至1%以内，说明收敛精度得到了进一步提升.可以看出，改进搜索算法后的反向有限元法具有更高的计算效率和精度.

图8 遗传算法和随机游走算法相结合的反算结果

表2 初始参数值与遗传算法和随机游走算法相结合的反算结果

4 结论

本文以环形缺口小冲杆试样为研究对象，从试样的径厚比、缺口深度以及有无预制裂纹3个方面，系统地研究了内聚力参数对模拟结果的敏感性，得到了优化的缺口样品几何尺寸.在此基础上，采用预先选取的内聚力模型参数进行模拟，并将载荷位移曲线作为目标结果；利用改进算法后的反向有限元法来提取模型参数，验证了反向有限元法的准确性.主要结论如下：

(1) 对于直径3 mm的无缺口小冲杆试样，在加载后期，最大主应力始终停留在试样下表面的确定位置，距离圆心的距离为0.3 mm.将此处设为环形缺口位置较为合理.

(2) 试样直径与厚度的比例越大，缺口深度越大，对内聚力参数的敏感性越高，反向有限元法的结果就越精确.但是缺口深度过大，内聚力模型的收敛性会下降.从敏感性和收敛性两方面综合考虑，缺口深度选为厚度的一半.

(3) 预制裂纹对参数敏感性的提升较小，且实际加工困难，在样品设计时不作考虑.

(4) 基于遗传算法的反向有限元计算表明，反算所得参数与预先选取的参数之间误差在6%以内，模拟和目标试验的曲线基本重合，验证了反向有限元法的准确性.将遗传算法和随机游走算法相结合，对反向有限元法进行优化，误差则减小到1%以内，模拟和目标试验的曲线几乎完全重合.这说明，改进搜索算法后反向有限元法的效率与精度都得到了大幅提升，也进一步验证了样品设计的有效性.