在爆炸冲击领域，材料的动态破坏问题一直是研究的热点和难点^[1]。相关研究一般需要从材料实验入手，常见的实验有：拉伸、压缩及剪切等准静态实验，Hopkinson压杆、泰勒杆及夏比冲击试验等动态实验^[2-6]。由于测量技术的限制，实验一般只能测定较宏观的物理量，而全面深刻地表征材料的力学特性，还需要结合高精度的有限元模拟。

材料失效模型是模拟结构破坏的关键，相关研究很多。通过对铝合金试件进行的大量研究，有了一个基于材料宏观破坏失效应变的BW(Bao Yingbin-Wierzbickis)失效模型，这种失效模型将断裂分为3种模型：韧性断裂模式、剪切断裂模式及混合断裂模式，适用范围较广^[7-8]。Y.W.Lee^[9]基于平面应力状态模型实验对BW失效准则进行了修正，通过引入断裂主应变空间概念，给出了基于失效主应变的FFLD(fracture forming limit diagram)失效模型。M.Luo等^[10]提出了一种基于拉伸与弯曲耦合作用下的MMC(modified Mohr-Coulomb)断裂失效模型。而工程上应用最广的是JC本构和失效模型^[11]，该模型是很多学者进行材料力学性能研究的基础^[12-13]。对于金属材料，JC本构模型能较好地反映材料塑性强化、应变率硬化和温度软化效应，而且形式相对简单。相应的JC失效模型则充分考虑了应力三轴度、应变率和温度对材料失效应变的影响，模型参数具有明确的物理意义。然而，针对工程上的具体材料，相应的模型参数一般需要进行一系列材料实验才能得到，尤其应变率相关的参数的实验，复杂且价格昂贵。

Q370d钢强度高、韧性好，在桥梁、船舶等焊接结构中应用广泛^[14]。本文中，结合夏比冲击试验^[15]和数值模拟，利用正交设计和回归分析对Q370d钢的JC失效模型参数进行研究，试图探索获取JC失效模型参数的方法。

1.   有限元模型

1.1   网格模型

图 1为V型缺口试样的有限元模型图，笛卡尔坐标系Oxy平面为冲击平面。试件横截面厚度为B，高度为W，长度为L，V型缺口试件的深度为a。标准试样无V型缺口。

图  1  夏比冲击试验V型缺口试件的几何、网格图

Figure  1.  Geometry and grid for Charpy V-notch test specimen

下载: 全尺寸图片 幻灯片

撞击平面内，区域1为缺口及刀刃撞击区，应力梯度变化大，单元需准确模拟材料的塑性及破坏，采用精细网格(0.1 mm)；区域2为弹性应力波传播区，由3级1:2过渡为粗网格；区域3为砧座接触区，为了准确模拟砧座接触力，由2级2:1过渡为较密网格。横向采用全均匀网格，为了保持断裂区单元较好的形状比，横向网格尺寸统一为0.1 mm。虽然这样处理会使弹性区域网格形状比过大，但由于横向应力梯度变化较小，这种网格不会导致较大误差。摆锤与砧座简化为解析刚体，与试件之间定义自动接触。

1.2   材料模型

材料采用JC本构模型和失效模型进行描述，表达式分别为：

式中：$\dot{\varepsilon}_{0}$为准静态实验应变率，一般取10^-4 s^-1，T^*=(T－T_r)/(T_m－T)，T_r时室温，T_m为材料熔点，$\overline{\varepsilon}_{D}^{\mathit{\boldsymbol{pl}}}$为材料的临界失效应变函数，σ^*为应力三轴度(平均应力σ_m与Mises应力σ的比值)，其他为材料参数。忽略温度因素，Q370d钢的JC本构模型参数分别为：d₁=350.40 MPa，d₂=565.16 MPa，n=0.483，d₃=0.019。其中，d₁、d₂、n通过圆棒准静态拉伸实验获得，而应变率参数d₃由于实验条件有限，参考了文献[16]。下面对3种厚度的无缺口试件进行实验和数值模拟。如图 2所示，试件厚度W分别为4、5和6 mm，每组厚度为两个试件，共6个，实验结果取平均值。在300 J的摆锤试验机上进行实验，摆锤初始冲击速度为5.2 m/s。定义实验后试件中横截面上的两个特征参数为d_t和d_c，试件长度方向上的特征参数为l(见表 1)。由表 1可以看到：3种厚度试件的冲击功值，数值模拟与实验的误差在3%之内，变形后特征参数的误差不超过7%，这表明数值模型和本构参数有较高精度。

图  2  标准试件的最终变形

Figure  2.  Final deformation of the standard specimens

下载: 全尺寸图片 幻灯片

表  1  标准试件的实验和数值模拟结果

Table  1.  Experimental and simulated results of the standard specimens

W/mm	Ak, exp/J	Ak, num/J	εAk/%	dt, exp/mm	dt, num/mm	εdt/%	dc, exp/mm	dc, num/mm	εdc/%	lexp/mm	lnum/mm	εl/%
4	181.7	180.3	-0.8	6.86	7.26	5.8	2.52	2.69	6.7	24.93	23.89	-4.2
5	226.6	224.7	-0.8	8.17	8.71	6.6	3.22	3.39	5.3	25.37	24.04	-5.5
6	277.8	270.6	-2.6	9.40	9.94	5.7	3.91	4.09	4.6	25.33	23.87	-6.1

下载: 导出CSV 

| 显示表格

1.3   断裂区网格尺寸讨论

为了考察不同网格尺寸对数值模拟结果的影响，取一组典型的失效模型参数(D₁=0.177, D₂=3.825, D₃=1.845，D₄=0)，采用二维模型，对试件的断裂区分别划分5种不同尺寸网格，讨论模型对网格的敏感性和依赖性，结果对比见表 2。表中，Δd为最大网格尺寸，A_k1为不考虑失效模型的冲击功，A_k2为考虑失效模型的冲击功，F_m为摆锤的最大冲击力。

表  2  不同断裂区网格尺寸下的冲击功、摆锤冲击力

Table  2.  Impact energy and pendulum forcein plane models with different mesh sizes

Δd/mm	Ak1/J	Ak2/j	Fm/kN
0.300	14.69	9.45	0.417
0.200	14.69	8.12	0.334
0.100	14.78	7.14	0.277
0.075	14.81	7.36	0.291
0.050	14.81	6.99	0.275

下载: 导出CSV 

| 显示表格

可以看到：在不引入失效模型，冲击功对网格敏感性较小，断裂区网格最大尺寸为0.3 mm就能达到较高精度，冲击功变化值不到1%。而在考虑材料失效的情况下，最大网格尺寸在小于0.1 mm后趋于稳定，冲击功值变化小于5%，同时摆锤的最大冲击力也随之趋于稳定。因此，综合考虑计算资源以及计算精度，断裂区的最大网格尺寸选为0.1 mm。

2.   失效模型参数拟合

2.1   参数分析

材料的失效过程一般包括损伤起始和损伤演化两个过程。如图 3所示，在F点及F点之前D=0，材料单元没有发生任何损伤，在G点D=1，材料单元完全破坏，在FG段之间的点，则将材料刚度乘以缩减系数(1~D)，模拟材料的塑性软化。

图  3  材料破坏的典型应力应变曲线

Figure  3.  Typical stress-strain curvesduring the material damage progress

下载: 全尺寸图片 幻灯片

对于JC失效模型，损伤起始方程(式(2))为3个表达式相乘，分别表征应力三轴度、应变率和温度的影响，各表达式之间参数不相互影响，共包含D₁~D₅等5个参数。对于应变率参数D₄，一般需要进行不同应变率的实验，而夏比冲击试验应变率水平一般在10³~5×10³ s^-1范围，不能反映材料在不同应变率下的特性，因此难以单纯通过夏比冲击试验对D₄进行优化拟合。而温度参数D₅对实验条件要求较高，所以在忽略温度影响的情况下，本文中只对D₁、D₂、D₃等3个参数进行优化拟合。

对于损伤演化参数，一般采用等效塑性变形u_pl或断裂能G_f来定义，对于金属材料一般建议采用u_pl，其定义如下：$\dot{\overline{u}}_{\mathrm{pl}}=L_{\mathrm{e}} \dot{\overline{\varepsilon}}_{\mathrm{pl}}$，其中L_e为单元的特征长度。当单元的等效塑性变形u_pl到达u_{pl, f}时，单元材料完全破坏。但很难通过实验直接得到u_{pl, f}具体值，一般取单向拉伸实验得到的延伸率乘以单元长度^[17]。本文中单元长度为0.1 mm，材料的延伸率通过拉伸实验测得为0.21，因此u_{pl, f}取0.021。

2.2   参数拟合

这里设计不同参数水平，通过数值模拟得到参数与冲击功的分析样本，再通过回归分析得到冲击功与模型参数之间的回归方程，最后结合夏比冲击试验冲击功来求解回归方程组。各参数取5个水平，D₁∈(0.15, 0.2, 0.25, 0.3, 0.35), D₂∈(3, 4, 5, 6, 7), D₃∈(－1.8, －2.2, －2.6, －3.0, －3.4)，共25种组合方案，对3种厚度(W=10, 6, 4 mm)的V型缺口试件，共进行了75次计算，试件的冲击吸收功A_k见表 3。

表  3  JC失效模型的正交设计试验表

Table  3.  Orthogonal design parameters for JC failure model

No.	D1	D2	D3	Ak/J
				W=10 mm	W=6 mm	W=4 mm
1	0.15	3	-1.8	151.57	89.04	57.61
2	0.15	4	-2.2	130.82	79.83	53.25
3	0.15	5	-2.6	113.26	67.86	46.08
4	0.15	6	-3.0	95.71	56.69	38.84
5	0.15	7	-3.4	81.50	47.67	32.54
6	0.20	3	-2.6	97.34	56.13	44.43
7	0.20	4	-3.0	92.88	54.81	34.06
8	0.20	5	-3.4	91.53	48.69	33.76
9	0.20	6	-1.8	287.76	185.80	107.57
10	0.20	7	-2.2	253.14	160.20	100.67
11	0.25	3	-3.4	87.83	48.38	31.22
12	0.25	4	-1.8	232.39	137.77	84.92
13	0.25	5	-2.2	203.56	122.17	77.73
14	0.25	6	-2.6	169.53	103.94	67.87
15	0.25	7	-3.0	150.13	87.96	57.72
16	0.30	3	-2.2	160.57	91.74	57.25
17	0.30	4	-2.6	152.92	87.35	54.97
18	0.30	5	-3.0	141.67	80.70	50.92
19	0.30	6	-3.4	129.61	73.81	46.59
20	0.30	7	-1.8	288.63	193.92	114.43
21	0.35	3	-3.0	130.85	73.26	44.74
22	0.35	4	-3.4	128.04	71.75	43.96
23	0.35	5	-1.8	286.80	179.69	103.27
24	0.35	6	-2.2	268.54	162.92	98.42
25	0.35	7	-2.6	230.74	138.54	87.47

下载: 导出CSV 

| 显示表格

采用线性变换对所有自变量及因变量进行归一化处理，分别将各参数取值的范围映射到区间[-1, 1]上，变换后的变量加上横线来区分。考虑参数之间的相互影响，引入D₁、D₂、D₃，二次非线性项D₁²、D₂²、D₃²和交叉乘积项D₁D₂、D₁D₃、D₂D₃。利用SPSS程序进行逐步回归，得到：

每种厚度试件为3个，共9个试件T1~T9，冲击实验结果如图 4~5所示。T1~T9的冲击功分别为208.3、203.1、199.4、104.9、106.7、112.0、55.5、55.7、56.7 J。

图  4  V型缺口试件的最终变形

Figure  4.  Final deformation of V-notch specimens

下载: 全尺寸图片 幻灯片

图  5  实验冲头力

Figure  5.  Experimental punch loads

下载: 全尺寸图片 幻灯片

将实验得到的冲击功各自归一化后代入回归方程的左边，则由回归方程(3)得到以D₁、D₂、D₃为变量的非线性方程，采用最小二乘法求解，得到近似解：D₁=(D₁=0.587, D₂=1.975, D₃=－3.586)。可以看到：近似解的值均不在参数设定的取值范围内，为外插型。为了进一步提高计算精度，在近似解附近增加了一个三水平D₁∈(0.5, 0.6, 0.7)、D₂∈(1.0, 2.0, 3.0)、D₃∈(－3.0, －3.5, －4.0)的正交设计试验，按照前面的方法可以求得近似解D₂=(D₁=0.788, D₂=1.835, D₃=－3.275)。

将D₂代入数值模型计算，得到W为10、6、4 mm的冲击功分别为187.62、105.72、59.77 J，与实验平均值200.6、107.87、56.97 J相比，误差分别为-6.5%、-2.0%、4.9%，基本上吻合。而不同厚度试件的误差不同，可能源自其应力三轴度分布水平的差别。从冲击力曲线(见图 6)看，数值模拟结果比实验的大，原因可能为在断裂破坏时，截面附近材料瞬间温度升高，材料发生温度软化效应，导致冲头撞击力的峰值下降，冲头位移增大，而计算模拟中没考虑温度影响。另外，JC失效模型与材料实际破坏情况可能有一定差异。图 7为V型缺口试件计算的最终变形。

图  6  计算和实验的冲头力

Figure  6.  Computed and experimental punch loads

下载: 全尺寸图片 幻灯片

图  7  V型缺口试件计算的最终变形及观测点

Figure  7.  Computed final deformation of V-notch specimensand location of observation points

下载: 全尺寸图片 幻灯片

3.   断裂截面特性分析

以4 mm厚试件为例，图 8给出了试件断裂截面上部分观测点单元(见图 7)等效应力-等效塑性应变的关系，这些曲线与材料的JC本构方程(式(1))比较吻合，一些观测点在破坏前还经历了卸载。发生材料破坏单元的应力最终都减小到零，单元发生初始损伤的点是曲线出现折角的地方，对应曲线上标点的位置，随后材料损伤累积直到破坏。忽略弹性变形能，由${E_{{\text{pl}}}} = \int_0^{{{\overline \varepsilon }_{{\text{pl}},{\text{f}}}}} {\overline \sigma } {\text{d}}{\overline \varepsilon _{{\text{pl}}}}$和${G_{\text{f}}} = \int_{{{\bar \varepsilon }_{{\text{pl}}}},0}^{{{\bar \varepsilon }_{{\text{pl}},{\text{f}}}}} {\overline \sigma } {\text{d}}{\overline \varepsilon _{{\text{pl}}}}$，可以得到单元的塑性变形能密度E_pl和单元的断裂能密度G_f(损伤点后的曲线面积)，两者比$\varphi=\overline{G}_{\mathit{\boldsymbol{f}}} / \overline{E}_{\mathit{\boldsymbol{pl}}}$, 可以作为材料抗脆性断裂的衡量指标，比值越大表明材料抗脆性断裂能力越强。计算得到W为10、6、4 mm试件的φ值依次为0.159、0.157、0.148，可见随着截面尺寸的减小，其抗脆性断裂能力减小。

图  8  观测点的等效塑性应变与等效应力曲线

Figure  8.  Curves of equivalent stress vsequivalent plastic strain at observation points

下载: 全尺寸图片 幻灯片

应力三轴度是影响材料损伤失效的重要因素，截面上应力三轴度不同可能导致不同厚度试件误差大小有所差异。从试件变形图看到，试件缺口端附近观测点始终处于拉伸状态(左)，另一端由压缩状态进入拉伸状态(右)，分别选取两端破坏的单元作为观测点，给出了应力三轴度-等效塑性应变曲线，如图 9所示。可见由于单元失效，曲线末端近似垂直线段。不同截面尺寸上，各单元的应力三轴度变化趋势基本一致。同一观测点，试件截面越窄，应力三轴度水平越小，且曲线变化趋势有明显的滞后现象。

图  9  试件断裂截面上单元应力三轴度曲线

Figure  9.  Stress triaxiality curves of observation elements at the fracture cross-sections of specimens

下载: 全尺寸图片 幻灯片

4.   结语

材料宏观破坏时吸收的能量是表征材料破坏力学性能的一个重要指标，探究夏比冲击试验试件的冲击功与损伤失效模型参数之间的关系是研究材料动态力学性能的一种思路。本文中在对有限元模型及本构参数验证的基础上，结合正交设计和回归分析，通过夏比V型缺口试验及数值模拟，探讨了冲击吸收功与失效模型参量之间的关系，得到了Q370d钢一组具有工程适用价值的JC失效模型参数(D₁=0.788, D₂=1.835, D₃=－3.275)。这种结合正交设计和回归分析研究材料动态失效参数的思路是可行的。