考虑单元网格尺寸影响的RTCL准则修正及应用研究

刘赫崴，刘 昆，王秀飞，王正耀

(江苏科技大学 船舶与海洋工程学院，江苏 镇江 212003)

随着计算机硬件及有限元算法的不断发展，有限元数值仿真方法已成为评估工程结构安全性的重要手段。对于船舶碰撞搁浅过程中的结构破坏问题，选择合适的失效准则是保证仿真计算结果合理准确的关键。

国内外学者开展了相关研究工作，其中使用较多的临界等效塑性应变准则认为当等效塑性应变达到临界值时，材料出现失效破坏。该准则形式简单，输入参数少，但在应力状态对失效的影响方面考虑不足。Alsos等[1]结合Hill[2]的局部颈缩分析准则和Bressan等[3]的剪应力准则而提出Bressan & Wlliams-Hill(BWH)失效准则，并将局部颈缩的出现作为断裂产生的判定依据，避免了局部颈缩区内过大变形，因而具有网格尺度敏感度低的优点，但仅适用于平面应力单元且忽略了弯曲对失效的影响。Johnson等[4]提出的通用断裂应变模型同时考虑应力三轴度、应变率和温度效应对失效的影响，具有物理意义清晰、参数易标定等优点，但其不能准确描述较高应变率下的断裂行为。Tørnqvist[5]结合Rice-Tracey准则[6]和Cockcroft-Latham[7]准则提出RTCL准则，相较于上述准则，其定义较为简单，仅需要失效应变εf作为唯一参数，且考虑全局应力三轴度，对应力状态复杂的船体结构失效问题[8-9]模拟效果较好。然而相关研究[10-11]表明，在使用RTCL准则进行断裂失效模拟时，不同的网格尺寸下仿真结果存在较大差异，需要改变失效参数εf的设定，才能较好地模拟网格尺寸差异较大的船体碰撞仿真。

综上，本文考虑网格单元尺寸的影响，对现有RTCL失效准则进行修正，建立了修正的RTCL失效准则(I-RTCL)，并通过二次开发嵌入大型非线性有限元程序中。在此基础上，充分考虑船舶与海洋工程碰撞问题特点，设计开展了不同应力状态下的单轴拉伸试验与仿真分析，验证了提出的I-RTCL失效准则在不同应力状态、不同网格尺寸下的适用性。

1 修正的RTCL失效准则

1.1 RTCL失效准则

一般认为，以静水应力与等效应力的比值构造出的无量纲量参数应力三轴度η可以较好地表征塑性变形过程中应力状态对韧性断裂的影响，并经常被用于失效准则的构建[12]。Bao等[13]发现在负应力三轴度条件下，断裂由剪切机制掌控；在高应力三轴度范围内(η≥1/3))，断裂由孔长大机制所主导；而在介于二者之间的低应力三轴度条件下，断裂由两种机制共同作用产生。

(1)

此外，Cockcroft-Latham准则能够准确预测韧性剪切断裂(η≤1/3)，Rice-Tracey准则可以正确模拟孔洞生长导致的断裂(η≥1/3)，Tørnqvist将二者组合并结合Bao等[14]有关断裂截止值的研究(η≤-1/3时不再发生断裂)，得到适用于全局应力三轴度的RTCL准则表示为

(2)

其中，

(3)

式中：D为累计损伤，当一单元在其厚度方向上的所有积分点都满足D≥1时，该单元删除；εf为失效应变，是RTCL准则唯一输入参数；f(η)为关于应力三轴度的函数，在不同应力三轴度区间内累计损伤函数不同。

1.2 RTCL失效准则修正

由式(2)可知，RTCL准则定义的累计损伤与失效应变密切相关，而在有限元软件中网格尺寸L又影响失效应变的设定，因此为降低RTCL准则在仿真中因L产生的误差，增设L为影响失效的自变量并以此构建I-RTCL准则。参照文献[15]制造平板试件与光滑圆棒试件，开展单轴拉伸试验及仿真，按表1建立不同网格尺寸、类型的有限元模型，尺寸、试件及部分有限元模型如图1所示。

表1 有限元模型网格尺寸与种类Tab.1 Mesh sizes and types of FEMs

图1 尺寸、试件及有限元模型(mm)Fig.1 Dimensions, specimens and FEMs(mm)

以平板试件5 mm模型及光滑圆棒1 mm模型为例，多次尝试设定失效应变直至仿真与试验的载荷-位移曲线吻合，仿真失效前的应力云图与试验同样出现颈缩现象，如图2所示，此时得到对应网格尺寸下的失效应变数值分别为0.62与1.05，同理获得多组数据并将结果汇总如表2所示。

图2 试验与仿真载荷-位移曲线Fig.2 Load-displacement curve in experiment and simulation

表2 不同网格尺寸下有限元模型对应失效应变Tab.2 Failure strain of FEMs with different mesh sizes

(4)

图3 失效应变与网格尺寸关系Fig.3 Relationship between failure strain and element size

对于平板试件，a=0.717 2，b=-1.219，c=0.534 7；对于光滑圆棒，a=0.738 7，b=-0.506，c=0.342 2。

(5)

(6)

此时，I-RTCL失效准则的累计损伤表达式由式(2)转变为

(7)

1.3 I-RTCL失效准则的VUMAT开发

利用Fortran语言编写用户子程序VUMAT以自定义I-RTCL失效准则，图4展示了本文编写子程序判定单元是否失效的计算流程，其在每个增量步内完成如下工作：

(1) 编写的VUMAT子程序可以从ABAQUS主程序中获得材料参数、状态变量等参数信息用于材料的失效判定。在获得参数后，首先假定应变增量为弹性应变，使用广义胡克定律计算试探应力后计算其等效应力，将等效应力与屈服应力对比判断此增量步下该积分点是否处于屈服状态。

(2) 如等效应力小于屈服应力则未屈服，积分点不产生累积损伤，所得试探应力即为此增量步结束时的应力，之后更新能量和状态变量，子程序结束，进入下一增量步；如等效应力大于屈服应力则发生屈服，表示材料产生塑性变形，此时需要将弹性试探应力返回到屈服面上来计算出真实的应力。由所得真实应力计算应力三轴度，根据此时的应力状态利用I-RTCL准则计算累计损伤D，判断积分点是否失效。

(3) 如D<1则未发生失效，子程序更新能量和状态变量后，进入下一增量步；如D≥1则发生失效，失效积分点删除，当单元厚度方向上所有积分点均被删除时，控制单元失效的状态变量赋值为0，单元删除，进入下一增量步。

图4 VUMAT子程序计算流程图Fig.4 VUMAT subroutine calculation flow chart

2 I-RTCL失效准则试验验证及分析

2.1 拉伸试件设计

为模拟不同应力状态下的金属失效，按图5尺寸制作多种单轴拉伸试件，按表3划分为不同网格尺寸与类型的有限元模型，试件与部分有限元模型如图6所示，包含平板剪切试件：SST-1，SST-2；缺口平板试件：PNT-Rn，n=6，n=10，n=14(n代表缺口半径，下同)；缺口圆棒试件：SNB-Rn，n=6，n=9，n=18。

由图5可知，SST-1在被拉伸时连接段(即试件最先失效处)的受力状态接近于纯剪切，而SST-2的连接段相较于SST-1存在2.5 mm的错位，因此连接段除受剪力外还会受到拉力。通过控制PNT-Rn与SNB-Rn的缺口半径来控制连接段受力状态，总体上，随着缺口半径n的增大，应力三轴度η不断增大。

表3 有限元模型网格尺寸与类型Tab.3 Mesh sizes and types of FEMs

图5 拉伸试件尺寸(mm)Fig.5 Dimension of the specimens (mm)

图6 拉伸试件及有限元模型(mm)Fig.6 Tensile specimen and FEMs (mm)

2.2 拉伸试验与仿真概况

对相同网格尺寸的模型分别使用I-RTCL准则与RTCL准则开展拉伸仿真，对比修正前后仿真断裂时刻拉伸位移与试验值的误差以验证I-RTCL优势。考虑到平板与圆棒试件结构上差异较大，在仿真时分别选取用式(5)与式(6)修正的失效准则。

2.3 结果分析

2.3.1 失效单元应力状态

以NRB-Rn为例，其失效单元应力三轴度变化曲线如图7所示，在加载初期曲线存在振荡情况，随后曲线稳定在均值附近，失效单元应力状态较为稳定，说明利用单轴拉伸试件可有效模拟某一应力状态下的金属断裂失效行为。

图7 失效单元应力三轴度变化曲线Fig.7 Stress triaxiality of failure element

表4 不同试件失效单元平均应力三轴度Tab.4 Mean stress traxiality of failure element in different specimens

2.3.2 载荷-位移曲线

选取SST-1，SST-2，PNT-R6与NRB-R6拉伸试验与仿真的载荷-位移曲线如图8所示，可以看到：

(1) 修正前后的仿真与试验测量载荷-位移曲线吻合均较好，但使用I-RTCL失效准则的仿真断裂时刻拉伸位移分布更为集中，说明在多种应力状态、多种网格尺寸下，相比原RTCL准则本文提出的I-RTCL失效准则精度更高。

图8 试验与仿真对应载荷-位移曲线Fig.8 Force-displacement relation in test and simulation

(2) 纯剪切试件SST-1的拉伸曲线如图8(a)、图8(b)所示，其仿真曲线后半段载荷明显高于试验结果，结合0.4 mm工况下失效单元应力三轴度随拉伸产生的变化(图9)不难看到：在仿真后期失效单元的三轴度仍不断上升，单元呈现拉伸状态，因此载荷仍呈上升趋势。由于载荷随拉伸位移增大持续上升，极限强度也表现出与断裂位移相似的尺寸效应，断裂拉伸位移更集中的I-RTCL准则对应仿真的强度极限也更加接近试验值。

(3) 表5为RTCL准则与I-RTCL准则在不同网格尺寸下仿真断裂位移及其误差，可以看到对于中低应力三轴度状态(SST-1，SST-2，PNT-Rn)，修正后仿真的误差明显降低；而对于高应力状态，存在修正效果一般的情况(NRB-R9-0.5 mm，NRB-R18-1 mm)；若修正前误差较小，经过修正误差仍保持在较低水平(SST-1-0.8 mm，PNT-R6-2 mm，NRB-R9-1 mm)。

(4) 对于缺口试件，缺口半径越小，断裂发生越快，分析产生此现象的原因主要如下：一是缺口半径越小，试件失效区域应力三轴度越高，越容易发生断裂；二是缺口半径越大，试件发生变形的连接段越长，其产生相同应变需要拉伸的位移也越长。

表5 仿真断裂位移及误差Tab.5 Fracture displacements and errors in simulations

表5 (续)

图9 0.4 mm计算工况失效单元应力三轴度Fig.9 Stress triaxiality of failure element in 0.4 mm calculation condition

3 结 论

(1) 相同模型，采用不同的网格尺寸进行仿真计算得到的断裂位移相差较大，网格尺寸对失效分析准确性影响较大，失效应变随网格尺寸的增大而降低，并逐渐趋于定值。

(3) 使用I-RTCL准则进行多应力状态、多网格尺寸下的仿真，其计算精度高于原准则，其中对于中低应力三轴度，提高显著；对于高应力三轴度，精度有所提升但存在修正效果一般的情况；此外，若修正前误差较小，经过修正误差仍保持在较低水平。