动态断裂韧性测试方法大致可以分为经验公式法、应变片法和光学方法，3种方法各有优劣。早期动态断裂实验研究的思路是将宏观量(加载点位移、载荷等)与微观量(动态应力强度因子)建立起联系。对于不同构型的试样，提出了相对应的经验公式^[1-2]，但经验公式法的局限性日益突出。随着光学测量技术的发展，裂纹尖端应变场可以通过光学测量手段测得，典型的测量手段有光弹法^[3]、焦散线法^[4]、数字图像相关法^[5]、相干梯度敏感干涉法^[6]等。结合高速摄影技术，裂纹扩展过程也可以被观察到^[7]。目前各种光学方法在材料断裂领域得到了广泛的应用。应变片法是将应变片粘贴在裂纹尖端附近，通过应变片信号计算应力强度因子。由于应变片可以直接得到裂纹尖端的信息，一般认为该方法较为准确，可以用来标定其他方法^[8]。应变片法相对于光学方法成本更低，但在高温或低温等特殊条件下不适用。Dally等^[9]最早采用应变片法对材料准静态与动态应力强度因子进行了测定，采用单应变片测量I型应力强度因子，给出了应变片粘贴角度满足一定条件时，应力强度因子与应变片信号之间的关系。Rittel^[10]将应变片法应用于Ⅱ型问题，采用双应变片，通过不同方向的应变来求解应力强度因子。目前虽然有学者讨论了应变片的角度^[11]、位置^[12]如何影响起裂时间测定，但应变片位置对应力强度因子测定精度的影响仅有定性分析，尚未细致讨论^[13-14]。文献[10]提出的双应变片法的计算误差取决于系数矩阵条件数。

本文分析应变片位置对系数矩阵条件数的影响，指出条件数最小的位置，同时采用数值模拟方法分析了不同条件数的情况下随机误差对动态应力强度因子计算解的影响，给出合适的应变片粘贴位置。

1 应变片法计算动态应力强度因子

如图 1所示的坐标系中，应变可以表示为^[10]：

$ {\varepsilon _{ij}} = \frac{1}{{\sqrt {2{\rm{ \mathsf{ π} }}{r_k}} }}\left[ {{K_Ⅰ}g_{ij}^{\rm{Ⅰ}}\left( {{\theta _k}} \right) + {K_Ⅱ}g_{ij}^Ⅱ\left( {{\theta _k}} \right)} \right] $

(1)

式中：下标i、j为应变分量方向，取值范围为1、2；下标k为应变片的编号，图 1中应变片1测量ε₁₁，应变片2测量ε₂₂。

应变分量与应力强度因子为：

$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _{11}}}\\ {{\varepsilon _{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {\frac{{g_{11}^{{\rm{Ⅰ}}}\left( {{\theta _1}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}{r_1}} }}}&{\frac{{g_{11}^{{\rm{Ⅱ}}}\left( {{\theta _1}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}{r_1}} }}}\\ {\frac{{g_{22}^{{\rm{Ⅰ}}}\left( {{\theta _2}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}{r_2}} }}}&{\frac{{g_{22}^{{\rm{Ⅱ}}}\left( {{\theta _2}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}{r_2}} }}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{K_{{\rm{Ⅰ}}}}}\\ {{K_{{\rm{Ⅱ}}}}} \end{array}} \right] $

(2)

式中系数矩阵为F，反解式(2)即可解出每一时刻的应力强度因子：

$ \mathit{\boldsymbol{F}} = {\mathit{\boldsymbol{K}}^{ - 1}}\mathit{\boldsymbol{\varepsilon }} $

(3)

平面应力状态下，式(2)中的g_ij(θ)表达式分别为：

$ \left\{ {\begin{array}{*{20}{l}} {g_{11}^{{\rm{Ⅰ}}}\left( {{\theta _1}} \right) = \frac{1}{E}\cos \frac{{{\theta _1}}}{2}\left[ {1 - \nu - (1 + \nu )\sin \frac{{{\theta _1}}}{2}\sin \frac{{3{\theta _1}}}{2}} \right]}\\ {g_{11}^{{\rm{Ⅱ}}}\left( {{\theta _1}} \right) = - \frac{1}{E}\sin \frac{{{\theta _1}}}{2}\left[ {2 + (1 + \nu )\cos \frac{{{\theta _1}}}{2}\cos \frac{{3{\theta _1}}}{2}} \right]}\\ {g_{22}^{{\rm{Ⅰ}}}\left( {{\theta _2}} \right) = \frac{1}{E}\cos \frac{{{\theta _2}}}{2}\left[ {1 - \nu + (1 + \nu )\sin \frac{{{\theta _2}}}{2}\sin \frac{{3{\theta _2}}}{2}} \right]}\\ {g_{22}^{{\rm{Ⅱ}}}\left( {{\theta _2}} \right) = \frac{1}{E}\sin \frac{{{\theta _2}}}{2}\left[ {2\nu + (1 + \nu )\cos \frac{{{\theta _2}}}{2}\cos \frac{{3{\theta _2}}}{2}} \right]} \end{array}} \right. $

(4)

式中:E为杨氏模量; ν为泊松比，将E替换为E/(1－ν²); ν替换为ν/(1－ν)即可得到平面应变状态下的表达式。

2 系数矩阵条件数的变化规律

由矩阵理论可知，式(3)解的精度取决于矩阵F^-1的条件数。矩阵条件数的定义为：

$ {\mathop{\rm cond}\nolimits} (F) = {\left\| \mathit{\boldsymbol{F}} \right\|_2} \cdot {\left\| {{\mathit{\boldsymbol{F}}^{ - 1}}} \right\|_2} $

(5)

式中为‖F‖₂矩阵F的2-范数，显然，矩阵F^-1的条件数等于F的条件数。当系数矩阵病态，即条件数很大时，应变分量的微小变化可能引起应力强度因子解的极大误差。因此讨论系数矩阵条件数的变化规律是十分必要的。

系数矩阵条件数的表达式十分复杂，其值取决于应变片的粘贴位置。为简便，本文从以下2种特殊情况分析：

1) 2个应变片中心点连线与裂尖共线，即θ₁=θ₂=θ，r₁/r₂=a(a≠1)时，矩阵的行列式为：

$ |F| = \frac{{1 - {\nu ^2}}}{{4{E^2}{\rm{ \mathsf{ π} }}r\sqrt a }}(\sin (2\theta ) + 2\sin \theta ) \approx 0 $

(6)

可见系数矩阵F接近奇异，这种情况是高度病态的，因此这种情况应当排除；

2) 2个应变片中心点到裂尖的距离相等，即r₁=r₂，此时条件数与θ₁和θ₂有关。对于如图 2所示的网格，θ₁和θ₂的取值范围是±168.75°，间隔11.25°。取r₁=r₂=5 mm，E=68.5 GPa，ν=0.33，采用Matlab编程计算系数矩阵的条件数。限于篇幅，本文列出±90°范围内系数矩阵条件数的值，计算结果分别如表 1、2所示。可以看出，系数矩阵条件数在θ₁与θ₂分别取±90°时取得极小值1，而当θ₁=0°，θ₂=±90°时，系数矩阵条件数取得极大值565.7。将表 1中大于10的数据用方括号表示，大于2的数据用圆括号表示，可见多数情况下系数矩阵条件数不大于2，但在极大值附近梯度很大。

推广到一般情况，即θ₁与θ₂的值确定，但r₁≠r₂的情形。设r₁ < r₂，且r₂/r₁=a，设此时系数矩阵为F^*，则有：

$ {\mathit{\boldsymbol{F}}^*} = \left[ {\begin{array}{*{20}{l}} {\frac{{g_{11}^{{\rm{Ⅰ}}}\left( {{\theta _1}} \right)}}{{\sqrt {2\pi {r_1}} }}}&{\frac{{g_{11}^{{\rm{Ⅱ}}}\left( {{\theta _1}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}{r_1}} }}}\\ {\frac{{g_{22}^{{\rm{Ⅰ}}}\left( {{\theta _2}} \right)}}{{\sqrt {2\pi a{r_1}} }}}&{\frac{{g_{22}^{{\rm{Ⅱ}}}\left( {{\theta _2}} \right)}}{{\sqrt {2{\rm{ \mathsf{ π} }}a{r_1}} }}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} 1&0\\ 0&{\frac{1}{{\sqrt a }}} \end{array}} \right]\mathit{\boldsymbol{F}} $

(7)

代入式(4)可以证明${\mathop{\rm cond}\nolimits} \left( {{\mathit{\boldsymbol{F}}^*}} \right) = \sqrt a \cdot {\mathop{\rm cond}\nolimits} \left( \mathit{\boldsymbol{F}} \right)$。同理，当r₂ < r₁时，令r₁/r₂=a也有类似的结果。因此对于确定的θ₁与θ₂，当r₁=r₂时系数矩阵条件数最小。

图 3为文献[10]中应变片的粘贴位置，应变片的坐标为：r₁=2.65 mm，θ₁=1.4°；r₂=3.77 mm，θ₂=－74.7°。代入式(2)可以计算出系数矩阵的条件数为26.57，可见此时系数矩阵较为病态。考虑到测量误差，实际贴片位置的系数矩阵条件数值可能更大，因此按文献[10]的位置贴片是不合适的。

3 应变片法误差分析

应变片的原理是通过电信号求解应变，而电路中不可避免存在噪声干扰，带来随机误差。本文讨论条件数不同的情况下，随机误差对应力强度因子求解的影响。本文基于分离式Hopkinson拉杆实验装置，提出1种改进的紧凑拉伸试样(modified compact tension shear, MCTS)用于Ⅱ型及复合型动态断裂加载测试^[15]。本文采用ABAQUS软件对MCTS试样受到Ⅱ型加载的问题进行有限元分析，根据定义，采用应力外推法^[16]计算出动态应力强度因子值。MCTS试样有限元模型如图 4所示，入射杆端部施加的应力波形为参考文献[17]中实测波形。

为了便于讨论，本文中将数值模拟直接求得的解称为数值解，通过理论公式推导出的结果称为计算解。采用应力外推法求得MCTS试样动态应力强度因子数值解，如图 5所示。将应力强度因子值代入式(2)即可求得一点应变的计算解。在此基础上附加随机误差后，代入式(3)反解出应力强度因子，并与其数值解进行对比。

图 6显示了随机误差幅值为应变最大值的5%时，应力强度因子计算解与数值解的差异。图 6(a)、(b)分别选取系数矩阵条件数取极小值1与极大值565.7的情况。由图可见，当条件数取极大值时，5%的随机误差就足以使动态应力强度因子的解面目全非。

当条件数取极小值时，动态应力强度因子解的误差很小。增大随机误差的数量级，幅值为应变最大值的25%时，求解结果如图 7所示。作为对比，系数矩阵条件数取在相同数量级，分别为1、5.9。可见动态应力强度因子的计算误差与条件数的数量级有关。当条件数位于同一数量级时，同样幅值随机误差作用下解的计算误差也在同一数量级。

4 改进的紧凑拉抻试样动态应力强度因子解

选取不同应变片位置计算动态应力强度因子。节点的局部坐标如表 5所示。将动态应力强度因子数值解代入式2可以求得每一节点的应变分量的计算解。应变分量计算解与数值解的对比如图 8、9所示。

可以看出，第1组节点条件数取得极小值，但应变计算解与数值解误差较大；第2组节点条件数更大，应变计算解与数值解误差较小。节点应变的最大绝对误差与最大相对误差如表 6所示。应力强度因子计算解与数值解对比见图 10、11。解的误差见表 7。

二者比较可以看出，虽然在±90°方向上系数矩阵条件数取得极小值，但应变计算解与数值解的误差较大；而在-45°方向附近应变计算解与数值解的误差较小，故后者精度更高。由此可见，要提高动态应力强度因子解的精度不仅需要考虑系数矩阵条件数，还需要考虑应变计算解与数值解的误差。

5 结论

1) 应变片法求解动态应力强度因子的过程中，系数矩阵条件数取决于应变片位置。大部分位置满足条件数小于2，而某些特定位置下条件数可能超过500，粘贴应变片时需要避开这些位置。

2) 随机误差作用下，动态应力强度因子解计算误差与随机误差之比的数量级与系数矩阵条件数的数量级相同。系数矩阵条件数大于500时，5%随机误差即可对动态应力强度因子的解产生巨大影响；而条件数小于10时，条件数的影响并不显著。

3) 动态载荷作用下，试样裂纹尖端应变场是数值解与式(2)在θ=-45°或θ=±90°方向能够较好地符合。在这些方向上粘贴应变片能得到精度较高的解。