梯形应力脉冲在弹性杆中的传播过程和几何弥散1)

图1柱杆及其横截面位移分量

Fig.1Cylindrical rod and the displacement components in its cross section

在此基础上 Rayleigh计算了有限长柱发生纵向谐振时横向运动给柱杆带来的附加动能.谐振波数越大,附加动能的影响越大,因此导致谐振周期降低. Rayleigh给出了谐振周期与谐振波数之间的解析表达式,通常被称为 Rayleigh横向惯性修正^[6].

事实上,Rayleigh在其著作中并没有推导考虑横向惯性效应的波动方程形式. Love^[14]从 Hamilton 原理出发给出了基于 Rayleigh横向运动规律的一维纵波波动方程. 在弹性动力学框架,系统的 Hamilton方程可表达为

(3)$\delta \int_{t_0 }^{t_1 } {\left( {T - V} \right) d t} = 0 $

其中系统动能 $T$ 为

(4)$ T = \int_V {\dfrac{\rho }{2}\left( {\dot {u}^2 + \dot {v}^2 + \dot {w}^2} \right)} d V = \\\ \qquad \int_L {d x} \int_A {\dfrac{\rho }{2}} \left[ {\left({\dfrac{\partial u}{\partial t}} \right)^2 + \nu ^2\left( {y^2 + z^2}\right)\left( {\dfrac{\partial ^2u}{\partial x\partial t}} \right)^2}\right]d A =\\\ \qquad \int_L {\dfrac{\rho A}{2}\left[ {\left( {\dfrac{\partial u}{\partial t}} \right)^2 + \nu ^2\kappa ^2\left( {\dfrac{\partial ^2u}{\partial x\partial t}} \right)^2} \right]d x}$

式中,$A$ 是柱杆的横截面积, $\rho $是材料质量密度, $\kappa $ 是截面对$x$轴的回转半径, $\kappa ^2 =\dfrac{1}{A}\int_A {\left( {y^2 + z^2} \right) d y d z} $.另一方面,系统势能 $V$ 为

(5)$V = \int_V {\dfrac{\sigma \varepsilon }{2}d V} = \int_L{d x\int_A {\left( {\dfrac{E\varepsilon ^2}{2}}\right)d A} } = \int_L {\dfrac{EA}{2}\left( {\dfrac{\partial u}{\partial x}} \right)^2 d x} $

式中, $\sigma $ 为轴向应力,$\sigma = E\varepsilon $;$E$为杨氏模量. 把式 (4) 和式 (5) 代入式 (3) 可得

(6)$\delta \int_{t_0 }^{t_1 } {d t} \int_L A\left\{ {\dfrac{\rho }{2}\left[ {\left( {\dfrac{\partial u}{\partial t}} \right)^2 + \nu ^2\kappa ^2\left( {\dfrac{\partial ^2u}{\partial x\partial t}} \right)^2} \right] +}\right. \\\ \qquad \left. { \dfrac{E}{2}\left( {\dfrac{\partial u}{\partial x}} \right)^2} \right\}d x = 0 $

对式 (6) 进行分部积分并应用初始条件和边界条件,整理得

(7)$\int_t {\int_L {\left( { - c_0^2 \dfrac{\partial ^2u}{\partial x^2} + \dfrac{\partial ^2u}{\partial t^2} - \nu ^2\kappa ^2\dfrac{\partial ^4u}{\partial x^2\partial t^2}} \right)} } \delta u d x d t = 0$

其中,$c_0 = \sqrt {E / \rho } $ 是初等理论下的杆纵波波速.注意到 $\delta u $在积分空间和时间上任意变化,为使积分为零,被积函数必须为零,因此得到

(8)$\dfrac{\partial ^2u}{\partial t^2} = c_0^2 \dfrac{\partial ^2u}{\partial x^2} + \nu ^2\kappa ^2\dfrac{\partial ^4u}{\partial x^2\partial t^2}$

上式即为 Rayleigh-Love横向惯性修正后的波动控制方程,与初等理论对照可知,右边第二项即代表横向惯性效应.

1.2 Rayleigh-Love 方程的弥散关系

由于方程 (8) 中横向惯性修正项 $\nu ^2\kappa ^2\dfrac{\partial ^4u}{\partial x^2\partial t^2}$ 的存在,使得杆中弹性纵波将不再以恒速 $c_{0}$传播,不同频率 $f$ (或波长 $\lambda $) 的谐波将以不同的相速度 $c$ 传播. 以行波解

(9)$u\left( {x,t} \right) = u_0 \mbox{e}^{\mbox{i}\gamma \left( {x - ct} \right)} $

其中 $\gamma \left( { = {2\pi } / \lambda } \right)$

为波数,代入式 (8),则得 $ - \gamma ^2c^2 = - c_0^2 \gamma ^2 + \nu ^2\kappa ^2\gamma ^4c^2$,整理后得到

(10)$\dfrac{c}{c_0 } = \dfrac{1}{\sqrt {1 + \nu ^2\kappa ^2\gamma ^2} } $

对于半径为 $a$ 的圆柱杆,$\kappa = a / {\sqrt 2 }$,将式 (10)右端做泰勒展开后,略去高于二次的幂项,则得

(11)$\dfrac{c}{c_0 } \approx 1 - \pi^2\nu ^2\left( {\dfrac{a}{\lambda }} \right)^2 $

此即为 Rayleigh-Love 理论的行波弥散特性方程,当泊松比 $\nu =0$或者 $a/ \lambda $ 远小于 1时,相速 $c=c_{0}$,退化到了初等理论的情形.

2 Rayleigh-Love 波动方程初边值问题

2.1 采用 Laplace 变换法求解 Rayleigh-Love 波动方程

针对 Rayleigh-Love波动方程的研究大都限于讨论其弥散关系,对于 Hopkinson杆实验中特定的边界冲击问题,在文献中找到两种解析解：Davies^[4]边界载荷为阶跃应力脉冲 $\sigma_0 H\left( t \right)$的解析解,以及 Naitoh 和 Daimaruya^[19]考虑升时过程的边界速度载荷 $v_0 \left( {1-{\rm e}^{- t /{t_0 }}}\right)H\left( t \right)$ 的解析解.本文研究梯形边界载荷作用下应力波的传播问题.

首先考虑线性斜坡加载问题,对于一根长度为 $L$的圆杆,初始处于静止状态;一端施加斜率为 $\alpha$的斜坡应力脉冲,另一端自由,如图 2所示.则波动问题的边界条件和初始条件表示为

(12)$\left.\begin{array}{ll} \left( {E\dfrac{\partial u}{\partial x}} \right)_{x = 0} = \alpha tH\left( t \right)\,, & \left( {E\dfrac{\partial u}{\partial x}} \right)_{x = L} = 0 \\\ \left( u \right)_{t = 0} = 0\,, &\left( {\dfrac{\partial u}{\partial t}} \right)_{t = 0} = 0 \end{array}\right\} $

式中 $H(t)$ 为单位阶跃函数.

图2

图2波动问题的分析模型

Fig.2The analytical model of the wave propagation problem

以杆的几何参数和材料参数定义以下特征参量：杨氏模量 $E$为特征应力,杆波速 $c_{0}$ 为特征速度,杆直径 $d=2a$为特征长度,$d/c_{0}$ 为特征时间. 利用这些特征量将式 (8)和式 (12) 无量纲化,得到

(13)$\left. \begin{array} \frac{\partial ^2\bar {u}}{\partial \bar {t}^2} = \frac{\partial ^2\bar {u}}{\partial \bar {x}^2} + \nu ^2\bar {\kappa }^2\frac{\partial ^4\bar {u}}{\partial \bar {x}^2\partial \bar {t}^2} \\\ \left( {\frac{\partial \bar {u}}{\partial \bar {x} \ }} \right)_{\bar {x} = 0}= \bar {\alpha }\bar {t} H\left( \bar {t} \right)\,, \left( {\frac{\partial \bar {u}}{\partial \bar {x} \ }} \right)_{\bar {x} = \bar {L}=0} \\\ \left( \bar {u} \right)_{\bar {t} = 0} = 0\,, \ \ \left( {\frac{\partial \bar {u}}{\partial \bar {t} \ }} \right)_{\bar {t} = 0} = 0 \end{array} \right\}$

其中带有上划线的量为无量纲化物理量. 对位移函数 $\bar {u}\left({\bar {x},\bar {t}} \right)$ 做 Laplace 变换,从时间域 $\bar{t}$ 变为像域 $s$,变换后的像函数以大写字母表示,即

(14)$U\left( {\bar {x},s} \right) = \int_0^\infty {\rm e}^{ - s\bar {t}}\bar {u}\left( {\bar {x},\bar {t} \ }\right) d \bar {t} $

结合初始条件,则原来的偏微分初边值问题 (13)在变换域上表现为常微分边界值问题

(15)$\left. \begin{array} \bal \dfrac{d ^2U}{d \bar {x}^2} - \dfrac{s^2}{1 + \nu ^2\bar {\kappa }^2s^2}U = 0 \\\ \left( {\dfrac{d U}{d \bar {x}}} \right)_{\bar {x} = 0} =\dfrac{\bar {\alpha }}{s^2}\,, \ \ \left( {\dfrac{d U}{d \bar {x}}} \right)_{\bar {x} = \bar {L}} = 0 \end{array} \right\}$

式 (15) 是典型的二阶常系数线性齐次微分方程问题,容易得到其解

(16)$U\left( {\bar {x},s} \right) = \\\ \Bigg \{ \bar {\alpha }\sqrt {1 + \nu ^2\bar {\kappa }^2s^2} \cosh \left[ s\left( {\bar {L} - \bar {x}} \right) / \sqrt {1 + \nu ^2\bar {\kappa }^2s^2} \right] \Bigg \} \Bigg /\\\ \qquad \Bigg [s^3\sinh \left( s\bar {L} / \sqrt {1 + \nu ^2\bar {\kappa }^2s^2} \right) \Bigg]$

只要对式 (16) 进行 Laplace 反变换,即可求得位移的原函数.以下围绕位移像函数的反变换过程做具体阐述.

2.2 基于留数定理的 Laplace 反变换

根据 Bromwich 反变换公式

(17)$\bar {u}\left( {\bar {x},\bar {t}} \right) = \dfrac{1}{2 \pi {\rm i}}\int_{\beta - {\rm i}\infty }^{\beta + {\rm i}\infty } {U\left( {\bar {x},s} \right){\rm e}^{s\bar {t}}d s} $

其积分路径如图 3红线所示.

图3

图3反 Laplace 变换的 Bromwich 积分路径及被积函数奇点分布

Fig.3The Bromwich integral contour for the inverse Laplace transform and the singularity distributions of the integrand

注意到式 (17) 中的被积函数 $U\left( {\bar {x},s} \right) {\rm e}^{s\bar {t}}$ 为单值函数,在像平面无支点,对 Bromwich积分路径添加一个大圆弧 (其上的积分贡献等于 0),构成回路积分 (图 3).根据回路积分规则,式 (17) 等于积分路径左侧所有奇点的留数和.如图 3所示,$U\left( {\bar {x},s} \right){\rm e}^{s\bar {t}}$的所有奇点都在虚数轴上,包括四重极点 $s_0 = 0$ 和简单极点 $s_n =n\pi {\rm i}\sqrt {\bar {L}^2 + n^2\pi ^2\nu ^2\bar {\kappa }^2} $ \ $\left( {n = \pm 1, \pm 2, \pm 3, \cdots } \right)$.采用留数定理计算积分式 (17),首先：在 $s_0 = 0$ 点处的留数

(18)${\rm Re} s\left[ {\bar {u}\left( {\bar {x},s} \right) {\rm e}^{s\bar {t}},s_0 }\right] = \\\ \qquad \frac{1}{3!}\mathop {\lim }\limits_{s \to 0} \frac{d^3}{d s^3}\left[ {s^4\bar {u}\left( {\bar {x},s} \right) {\rm e}^{s\bar {t}}} \right] = \\\ \qquad - \bar {\alpha }\bar {t}(\frac{\bar{t}^{2}}{6 \bar{L}}+\frac{\bar{x}^{2}}{2 \bar{L}}+\frac{\bar{L}}{3}-\bar{x}+\frac{v^{2} \bar{\kappa}^{2}}{\bar{L}})$

在 $s = s_n $ 点处的留数

(19)${\rm Re} s\left[ {\bar {u}\left( {\bar {x},s} \right){\rm e}^{s\bar {t}},s_n } \right] = \mathop {\lim}\limits_{s \to s_n } \left[ {\left( {s - s_n } \right)\bar {u}\left( {\bar {x},s} \right){\rm e}^{s\bar{t}}} \right] =\\\ \qquad \dfrac{\bar {\alpha }\bar {L}}{\pi ^3}\dfrac{\left( { - 1}\right)^n}{n^3}\dfrac{\cos \left[ {n\pi \left( {1 - \dfrac{\bar {x}}{\mathop{L}\limits^{\_\!\_}}} \right)}\right]}{\sqrt {{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi ^2\nu ^2\bar {\kappa }^{\,\,2}} } \cdot \\\ \left( {\sin \dfrac{n\pi \bar {t}}{\sqrt {{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi ^2\nu ^2\bar{\kappa }^{\,\,2}} }-{\rm i}\cos \dfrac{n\pi \bar {t}}{\sqrt {{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi^2\nu ^2\bar {\kappa }^{\,\,2}} }} \right)$

上式的留数带有虚部,但是当 ${\rm Re} s\left[ {\bar {u}\left({\bar {x},s} \right) {\rm e}^{s\bar {t}},s_n }\right]$ 和${\rm Re} s\left[ {\bar {u}\left( {\bar{x},s} \right) {\rm e}^{s\bar {t}},s_{ - n} } \right]$相加时,虚部相互抵消. 因此,时域位移场函数为

(20)$\bar {u}\left( {\bar {x},\bar {t}} \right) = {\rm Re} s\left[ {\bar {u}\left( {\bar {x},s} \right){\rm e}^{s\bar {t}},s_0 } \right]+ \\\ \qquad \sum_{n = 1}^\infty {\left( {{\rm Re} s\left[ {\bar {u}\left( {\bar {x},s} \right) {\rm e}^{s\bar {t}},s_n } \right] + {\rm Re} s\left[ {\bar {u}\left( {\bar {x},s} \right) {\rm e}^{s\bar {t}},s_{- n} } \right]} \right)} = \\\ \qquad - \bar {\alpha }\bar {t}\left( {\dfrac{\bar{t}^{\,\,2}}{6\,{\mathop{L}\limits^{\_\!\_}}} + \dfrac{\bar {x}^{\,\,2}}{2\,{\mathop{L}\limits^{\_\!\_}}} +\dfrac{{\mathop{L}\limits^{\_\!\_}}}{3} - \bar {x} + \dfrac{v^2\bar {\kappa}^{\,\,2}}{{\mathop{L}\limits^{\_\!\_}}}} \right) + \\\ \qquad \dfrac{2\bar {\alpha }\bar {L}}{\pi ^3}\sum_{n = 1}^\infty{\dfrac{\left( { - 1} \right)^n}{n^3}\dfrac{\cos \left[ {n\pi \left( {1 - \dfrac{\bar{x}}{{\mathop{L}\limits^{\_\!\_}}}} \right)} \right]}{\sqrt {{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi ^2\nu^2\bar {\kappa }^{\,\,2}} }}\cdot \\\ \qquad \sin \dfrac{n\pi \bar {t}}{\sqrt{{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi ^2\nu ^2\bar {\kappa }^{\,\,2}} } $

一旦得到杆沿轴向的无量纲位移解,其轴向应力以及应变可以相应得到,无量纲应力和应变为

(21)$\bar {\sigma }\left( {\bar {x},\bar {t}} \right) = \varepsilon \left( {\bar{x},\bar {t}} \right) = \dfrac{\partial \bar {u}\left( {\bar {x},\bar {t}} \right)}{\partial \bar {x}} = \bar{\alpha }\bar {t}\left( {1 - \dfrac{\bar {x}}{{\mathop{L}\limits^{\_\!\_}}}} \right) - \dfrac{2\bar {\alpha}}{\pi ^2}\cdot \\\ \qquad \sum_{n = 1}^\infty {\dfrac{\left( { - 1} \right)^n}{n^2}\dfrac{\sin\left[ {n\pi \left( {1 - \dfrac{\bar {x}}{{\mathop{L}\limits^{\_\!\_}}}} \right)} \right]}{\sqrt{{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi ^2\nu ^2\bar {\kappa }^{\,\,2}} }} \sin \dfrac{n\pi \bar{t}}{\sqrt {{\mathop{L}\limits^{\_\!\_}}^2 + n^2\pi ^2\nu ^2\bar {\kappa }^{\,\,2}} } $

2.3 梯形入射脉冲在杆中产生的位移、应变和应力解

对于总脉宽为 $\bar {T}$,上升沿和下降沿时间均为 $\bar {t}_r$,斜率均为 $\bar {\alpha }$,平台应力为 $\bar {\sigma}_0\,\,(\,=\!\!\!\bar{\alpha }\bar {t}_r)$ 的梯形应力脉冲,由于原问题 (13) 是线性齐次的,可应用叠加原理得到梯形脉冲加载下的位移 ${\hat{u}\left( {\bar {x},\bar {t}} \right)}/{\bar {\sigma }_0 }$,应变 ${\hat {\varepsilon }\left( {\bar {x},\bar {t}} \right)} /{\bar {\sigma }_0 }$,应力 ${\hat {\sigma }\left( {\bar {x},\bar{t}} \right)}/{\bar {\sigma }_0 }$ 的解析解,即

(22)$\left. \begin{array} \bal {\hat {u}\left( {\bar {x},\bar {t}} \right)} /{\bar {\sigma }_0 } =\\\ \qquad \bar {u}\left( {\bar {x},\bar {t}} \right)H\left( \bar {t} \right) - \bar{u}\left( {\bar {x},\bar {t} - \bar {t}_r } \right)H\left( {\bar {t} - \bar {t}_r } \right) - \\\ \qquad \bar {u}\left( {\bar {x},\bar {t} - \bar {T} + \bar {t}_r } \right)H\left( {\bar {t} - \bar {T} + \bar {t}_r} \right) + \\\ \qquad \bar {u}\left( {\bar {x},\bar {t} - \bar {T}} \right)H\left( {\bar {t} - \bar {T}} \right) \\\ {\hat {\varepsilon }\left( {\bar {x},\bar {t}} \right)}/{\bar {\sigma }_0 } = \\\ \qquad \bar {\varepsilon }\left( {\bar {x},\bar {t}} \right)H\left( \bar {t}\right) - \bar {\varepsilon }\left( {\bar {x},\bar {t} - \bar {t}_r }\right)H\left( {\bar {t} - \bar {t}_r } \right) - \\\ \qquad\bar {\varepsilon }\left( {\bar {x},\bar {t} - \bar {T} + \bar {t}_r }\right)H\left( {\bar {t} - \bar {T} + \bar {t}_r } \right) +\\\ \qquad \bar {\varepsilon }\left( {\bar{x},\bar {t} - \bar {T}} \right)H\left( {\bar {t} - \bar {T}} \right) \\\ {\hat {\sigma }\left( {\bar {x},\bar {t}} \right)}/{\bar {\sigma }_0 } = \\\ \qquad \bar {\sigma }\left( {\bar {x},\bar {t}} \right)H\left( \bar {t} \right)- \bar {\sigma }\left( {\bar {x},\bar {t} - \bar {t}_r } \right)H\left({\bar {t} - \bar {t}_r } \right) - \\\ \qquad \bar {\sigma }\left( {\bar {x},\bar {t} - \bar {T} + \bar {t}_r }\right)H\left( {\bar {t} - \bar {T} + \bar {t}_r } \right) +\\\ \qquad \bar {\sigma }\left( {\bar{x},\bar {t} - \bar {T}} \right)H\left( {\bar {t} - \bar {T}} \right) \end{array} \right\} $

3 算例

3.1 典型的弥散波形

设梯形应力脉冲脉宽 $\bar {T} = 12$,上升沿和下降沿时间 $\bar{t}_r = 1$,平台应力为 $\bar {\sigma }_0 = 0.004$,斜率 $\bar{\alpha } = 0.004$,弹性杆泊松比 $\nu = 0.3$,离加载端 $\bar {x}= 5$ 处的典型应力弥散波形如图 4所示,该结果为对式 (21)级数前 500 项求和得到. 与杆端输入的梯形脉冲 (在图中假设 $\nu=0.0$,求和 1000项得到)相比,横向惯性带来的应力脉冲几何弥散现象清晰可见,弥散波形与实验观察结果也一致.

图4

图4距离杆端 $\bar {x} = 5$ 处的应力时程曲线(红色虚线为初始的梯形脉冲)

Fig.4Stress time history curve at $\bar {x} = 5$ from the rod end (The red dashed curve is the original trapezoidal pulse)

值得一提的是,解析解 (式 (20) $\sim$ 式 (22))是一个无限 Fourier级数表达式,在具体的数值计算中级数求和的有限项数 $N$对计算结果有影响.图 5给出了上述算例中采用不同项数所得到的结果对比,可以看出,在 $N=150$ 以内,计算得到的脉冲波形有较明显不同,但是 $N>200$ 时,计算所得到波形形状接近收敛.在后续所有分析中,所采用的求和项数均大于或者等于 500,以确保收敛.

图5

图5有限求和项数 $N$ 对波形计算结果的影响

Fig.5The influence of finite summation terms $N$ on the calculation results of waveforms

3.2 与有限元计算结果的比较

考虑一个与实际测试问题相近,带有物理单位的具体问题：杆的弹性模量 $E= 200$ GPa, 密度 $\rho = 7800$ kg/m$^{3}$ (波速 $c_{0 } =5 064$ m/s),泊松比 $\nu =0.3$,加载端作用的梯形脉冲幅值的平台应力为 800 MPa,总加载历时 120 $\mu$s,包括上升沿和下降沿各 10 $\mu$s. 对于杆径为 5, 14.5, 37 和 74 mm 四种情况,分别计算杆端、离杆端 100, 200, 300, 400 和 500 mm处的应力脉冲波形,并与二维轴对称有限元数值模拟结果^[17]进行比较,如图 6所示. 可以看出,解析计算结果与有限元数值模拟的计算结果吻合良好.注意到文献的数值计算结果 (图 6左侧曲线)包含了完整的三维惯性效应,这表明在本算例所涉及的参数范围,运用 Rayleigh-Love一维波动近似理论描述的弹性杆中应力脉冲传播现象,能够足够精确地再现波动过程和几何弥散效应.

图6

图6不同直径的杆离杆端不同距离处的应力波形之比较(左图为有限元模拟结果^[17],右图为解析计算结果)

Fig.6Comparisons of stress waveforms at different distances from the rod end of rods with different diameters (The left curves are finite element simulation results^[17]and the right curves are analytical calculation results

4 参数研究

围绕几何弥散现象,我们比较了典型的弥散波形与输入波形,如图 7所示.可以看出,弥散造成的梯形入射波的形状变化主要体现在两个方面：弥散波的峰值应力 $\bar {\sigma }_{\max } $超过梯形波平台 $\bar {\sigma }_0 $,以及上升时间 $\bar {t}'_r $长于梯形波上升沿 $\bar {t}_r $.为定量描述梯形脉冲在圆杆中传播时的弥散效应,采用以下两个相对变化量(弥散参数)作为衡量指标,即：$\dfrac{\Delta \bar {\sigma }}{\bar {\sigma }_0 } = \dfrac{\bar {\sigma }_{\max }- \bar {\sigma }_0 }{\bar {\sigma }_0 }$,$\dfrac{\Delta \bar{t}_r }{\bar {t}_r } = \dfrac{\bar {t}'_r - \bar {t}_r }{\bar{t}_r }$.计算表明,这两个弥散参数将随着传播距离、泊松比、以及梯形波上升时间的变化而发生变化.

图7

图7脉冲波形弥散前后幅值和上升沿时间的变化(橙色虚线为弥散前的初始波形)

Fig.7Variations of amplitude and wave front rising time before and after pulse waveform dispersion (The orange dashed line is the initial waveformbefore dispersion

4.1 传播距离$\bar {x}$的影响

图 8给出 $\bar {\alpha } = 0.004$, $\nu = 0.3$时杆中不同位置的应力波随时间变化波形. 在输入脉冲的斜率和泊松比保持不变情况下,随着传播距离的增大,各频率成分的波之间会因几何弥散而产生相位差,从而导致脉冲的整体波形的幅值产生振荡,振荡幅值同时上升沿时间显著延长.随着传播距离的增大,应力振荡的幅值和上升沿时间延长量均持续增大,振荡幅值可达到应力平台的 20%.4.2节将给出弥散参数随空间位置的变化规律.

图8

图8传播距离的变化对脉冲波形的影响 ($\bar {\alpha } = 0.004$,$\nu = 0.3)$

Fig.8The effect of variation of propagation distance on pulse waveforms($\bar {\alpha } = 0.004$,$\nu = 0.3)$

4.2 泊松比 $\nu $ 的影响

图 9比较了不同泊松比情况下,$\bar {\alpha } = 0.004$的梯形脉冲在相同位置 ($\bar {x} = 5$) 上产生的弥散波形. 显然,随着泊松比的增大,杆中质点的横向运动的动能就会越大,从而导致脉冲波形的弥散更加显著,具体而言：波形幅值的振荡程度变大,上升时间延长.泊松比对弥散效应的影响,在 Rayleigh-Love横向惯性修正的波动方程 (8) 中直接体现.

图9

图9泊松比的变化对脉冲波形的影响 ($\bar {\alpha } = 0.004$,$\bar {x} = 5)$

Fig.9The effect of variation of Poisson's ratio on pulse waveforms ($\bar{\alpha } = 0.004$,$\bar {x} = 5)$

综合图 8和图 9的结果,给出两个弥散效应参数,$\left( {{\Delta \bar {\sigma }}/ {\bar {\sigma }_0 }} \right)$ 和 $\left({{\Delta \bar {t}_r } / {\bar {t}_r }} \right)$随泊松比和传播距离的变化关系,如图 10所示.可以明显看到泊松比和传播距离对弥散波形变形的定量增强效果.

图10

图10弥散波形的幅值振荡和上升沿时间随泊松比和传播距离的变化 ($\bar{\alpha }= 0.004)$

Fig.10The changes in amplitude oscillation and wave front rising time of the dispersed waveforms with Poisson ratio and propagation distance ($\bar {\alpha } = 0.004)$

4.3 脉冲斜率 $\bar {\alpha }$(或上升沿时间 $\bar {t}_r $) 的影响

固定梯形脉冲平台应力 $\bar {\sigma }_0 = 1$,泊松比 $\nu =0.3$,传播距离 $\bar {x} = 5$,改变梯形脉冲斜率 $\alpha=\infty$,10, 5, 2, 1, 0.5,上升沿时间依次为 $\bar {t}_r = 1/\bar {\alpha }=0.0$, 0.1, 0.2, 0.5, 1.0, 2.0,其中,$\bar \alpha=\infty$, $\bar t_r=0.0$ 是矩形脉冲的情形.图 11(a) 给出相同位置 ($\bar {x} = 5$)处记录到的弥散波形,对于梯形入射脉冲,随着脉冲上升斜率的增大,即上升沿时间的减小,波形振荡愈发剧烈.图 11(b) 给出 $\bar {x} = 5$处应力波形的相对振荡幅值随入射梯形脉冲升时的关系,当升时降为0时,振荡波形收敛于矩形脉冲的情形,此时振荡峰值超出平台应力约 25%,此结论与 Davies^[4]的结论一致.

图11

图11上升沿时间对几何弥散效应的影响 ($\nu = 0.3$,$\bar {x} = 5$)

Fig.11The effect of wave front rising time on the geometric dispersion effect ($\nu = 0.3$,$\bar {x} = 5)$

前面的算例 (图 5) 已证实：对于有限升时 ($\bar {t}_r = 1.0$)的边界冲击加载,当计算取足够多项数 ($N>200$) 时所得到的应力弥散波形趋于收敛. 对于更短促的升时,为达到计算收敛有必要取更多项数,特别是对于 $\bar{t}_r = 0.0$的矩形脉冲加载,通常认为至少在边界上载荷具有强间断,使用三角级数构筑解析解时不可避免地发生 Gibbs振荡;但是在本文所研究问题中,由于基本微分方程具有本征的几何弥散效应,矩形脉冲的入射波形在很短地传播距离上都将被拉缓,变为连续波形.图 12给出了 $\nu = 0.3$, $\bar {t}_r = 0.0$ 的阶跃载荷 (Davies问题) 在不同位置的弥散波形算例,计算分别取 $N =200$, 500 和 1000项进行计算,可见边界阶跃脉冲迅速弥散成为连续波,取 $N =500$或更高的计算项数能够准确给出弥散特征. 由此可见在实际计算中得到的振荡现象均属于物理本质现象,而非数值误差.

图12

图12阶跃脉冲在杆不同位置产生的弥散波形 (Davies 问题) ($\nu = 0.3$)

Fig.12Dispersed waveforms generated by step pulses at different positions of the rod (Davies problem) ($\nu = 0.3)$}

5 结论

本文基于考虑横向惯性效应的 Rayleigh-Love一维应力波理论,运用 Laplace积分变换-反变换方法分析杆承受边界冲击脉冲时的波动问题,获得如下结论：

(1) 给出了边界承受梯形加载脉冲时 Rayleigh-Love方程的解析解,表达为三角函数级数形式,在取 $N$=500$\sim $1000计算项数情况下,该解析解所给出的应力弥散波形完全收敛;

(2)通过将解析计算所得到的弥散波形与轴对称有限元模拟结果进行对比,证明 Rayleigh-Love理论能够较好地表征 SHPB实验所涉及的弹性圆杆中的应力波几何弥散效应,是一个实用的波动理论;

(3)基于梯形脉冲加载问题的解析解进行参数研究,揭示了传播距离、泊松比和脉冲升时变化对应力波弥散效应的影响.在 Hopkinson杆实验中,通过考虑杆的横向惯性效应对测试得到应力波形进行修正,是实现精确测试的重要步骤,大量文献对此技术进行了讨论^[21-25].作为一个三维问题的近似处理,Rayleigh-Love一维应力波方程形式简单,在绝大部分 Hopkinson杆实验中具有良好的精度,可以作为分析各种类型入射波的几何弥散效应的基础.本文给出 Rayleigh-Love杆在梯形加载条件下的级数解析解,有一定理论价值,在后续工作中,我们推导出用积分形式表达的任意载荷函数产生的弥散波的通解,可用于一般实验测试中的几何弥散修正.

致谢

作者感谢与宁波大学王礼立教授卓有成效的学术讨论.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Hudson

Dispersion of elastic waves in solid circular cylinders

Physical Review, 1943 , 63 : 46 - 51

Ultrasonic elastic modes in solid bars are investigated theoretically and experimentally using the plane wave expansion method to calculate the dispersion curves k=k(omega) for longitudinal, torsional, and flexural waves. The plane wave extension method allows to consider rods of circular and square cross sections. The technique, which has received attention in the study of photonic and phononic crystals, is adapted in order to identify the various types of modes. Results are compared with predictions from semi-analytical models. The numerical approximation is validated with the experimental determination of the time-frequency dispersion curves. The technique based on the plane wave expansion method presented here could be a numerical alternative used to determine the wave propagation and modal vibration with high precision in structures like bars and cylinders. Practical applications of this study could include the inspection of long-span engineering systems with bar or cylinder like characteristics.

[2]

Mindlin

, McNiven

Axially symmetric waves in elastic rods

Journal of Applied Mechanics, 1960 , 27 : 145 - 151

[3]

Hutchinson

Vibrations of solid cylinders

Journal of Applied Mechanics, 1980 , 47 ( 4 ): 901 - 907

Solid-state optical refrigeration uses anti-Stokes fluorescence to cool macroscopic objects to cryogenic temperatures without vibrations. Crystals such as Yb³⁺-doped YLiF₄(YLF:Yb) have previously been laser-cooled to 91 K. In this study, we show for the first time laser cooling of a payload connected to a cooling crystal. A YLF:Yb crystal was placed inside a Herriott cell and pumped with a 1020-nm laser (47 W) to cool a HgCdTe sensor that is part of a working Fourier Transform Infrared (FTIR) spectrometer to 135 K. This first demonstration of an all-solid-state optical cryocooler was enabled by careful control of the various desired and undesired heat flows. Fluorescence heating of the payload was minimized by using a single-kink YLF thermal link between the YLF:Yb cooling crystal and the copper coldfinger that held the HgCdTe sensor. The adhesive-free bond between YLF and YLF:Yb showed excellent thermal reliability. This laser-cooled assembly was then supported by silica aerogel cylinders inside a vacuum clamshell to minimize undesired conductive and radiative heat loads from the warm surroundings. Our structure can serve as a baseline for future optical cryocooler devices.

[4]

Davies

A critical study of Hopkinson pressure bar

Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1948 , 240 ( 821 ): 375 - 457

DOIURL[本文引用: 4]

[5]

Kolsky

An investigation of the mechanical properties of materials at very high rates of loading

Proceedings of the Physical Society B, 1949 , 62 ( 11 ): 676 - 700

Silk fibroin (SF) is a model candidate for use in tissue engineering and regenerative medicine owing to its bio-compatible mechanochemical properties. Despite numerous advances made in the fabrication of various biomimetic substrates using SF, relatively few clinical applications have been designed, primarily due to the lack of complete understanding of its constitutive properties. Here we fabricate microstructurally aligned SF sponge using the unidirectional freezing technique wherein a novel solvent-processing technique involving Acetic acid is employed, which obviates the post-treatment of the sponges to induce their water-stability. Subsequently, we quantify the anisotropic, viscoelastic response of the bulk SF sponge samples by performing a series of mechanical tests under uniaxial compression over a wide range of strain rates. Results for these uniaxial compression tests in the finite strain regime through ramp strain and ramp-relaxation loading histories applied over two orders of strain rate magnitude show that microstructural anisotropy is directly manifested in the bulk viscoelastic solid-like response. Furthermore, the experiments reveal a high degree of volume compressibility of the sponges during deformation, and also evince for their remarkable strain recovery capacity under large compressive strains during strain recovery tests. Finally, in order to predict the bulk viscoelastic material properties of the fabricated and pre-characterized SF sponges, a finite strain kinematics-based, nonlinear, continuum model developed within a thermodynamically-consistent framework in a parallel investigation, was successfully employed to capture the viscoelastic solid-like, transversely isotropic, and compressible response of the sponges macroscopically.

[6]

Kolsky

Stress Waves in Solids. Oxford : Clarendon Press , 1963

[7]

王礼立 , 王永刚 .

应力波在用SHPB研究材料动态本构特性中的重要作用

爆炸与冲击, 2005 ( 1 ): 17 - 25

URL[本文引用: 1]

在材料动态本构关系的研究中,不论是由波传播信息反求材料本构关系,即所谓解第二类反问题,还是利用应力波效应和应变率效应解耦的方法(如SHPB技术),应力波传播实际上都起着关键作用.在一般性讨论的基础上,就SHPB试验技术分析了应力波传播如何影响材料动态本构特性的有效确定.对于应力/应变沿试件长度均匀分布假定以及一维应力波假定,着重进行了分析.

( Wang

Lili

, Wang

Yonggang

The important role of stress waves in the study on dynamic constitutive behavior of materials by SHPB

Explosion and Shock Waves , 2005 ( 1 ): 17 - 25 (in Chinese))

URL[本文引用: 1]

[8]

刘孝敏 , 胡时胜 .

大直径SHPB弥散效应的二维数值分析

实验力学, 2000 , 15 ( 4 ): 371 - 376

( Liu

Xiaomin

, Hu

Shisheng

Two-dimensional numerical analysis for the dispersion of stress waves in large-diameter-SHPB

Journal of Experimental Mechanics , 2000 , 15 ( 4 ): 371 - 376 (in Chinese))

[9]

Shen

, Yin

Analysis of geometric dispersion effect of impact-induced transient waves in composite rod using dynamic substructure method

Applied Mathematical Modelling, 2016 , 40 ( 3 ): 1972 - 1988

[10]

Tyas

, Watson

An investigation of frequency domain dispersion correction of pressure bar signals

International Journal of Impact Engineering, 2001 , 25 ( 1 ): 87 - 101

[11]

Pochhammer

Uber die Fortpflanzungsgeschwindigkeiten kleiner Schwingungen in einem unbergrenzten isotropen Kreiszylinder

Journal Fur Die Reine Und Angewandte Mathematik, 1876 , 81 : 324 - 336

[12]

Chree

The equations of an isotropic elastic solid in polar and cylindrical coordinates, their solutions and applications

Transactions of the Cambridge Philosophical Society, 1889 , 14 : 250 - 369

[本文引用: 1]

[13]

Rayleigh

JWS

The Theory of Sound Vol 1, 2nd Edition

London: Macmillan and Co, 1894

[14]

Love

AEH

A Treatise of the Mathematical Theory of Elasticity

New York: Dover Publication, 1944

[15]

Krawczuk

, Grabowska

, Palacz

Longitudinal wave propagation. Part I---Comparison of rod theories

Journal of Sound and Vibration, 2006 , 295 ( 3-5 ): 461 - 478

[16]

Stephen

, Lai

, Young

et al .

Longitudinal vibrations in circular rods: A systematic approach

Journal of Sound and Vibration, 2012 , 331 ( 1 ): 107 - 116

A systematic method is developed for expressing the frequency squared omega(2) and the corresponding displacement fields of harmonic waves in a long thin rod as an even power series in qa, where q is the wavenumber along the rod and a is a representative transverse dimension. For longitudinal waves in a circular rod, the evaluation is reduced to algebraic recursion, giving coefficients analytically in terms of Poisson's ratio v, to many orders. The second nontrivial coefficient, corresponding to Rayleigh-Love theory in the present longitudinal case and Timoshenko theory in the flexural case, is thus put on a firm footing without reliance on ad hoc physical assumptions. The results are compared to available exact predictions, and shown to be accurate for moderate values of qa (5% accuracy for qa <= 1.5) with just two terms. Improvements based on the Rayleigh quotient guarantee positivity and the correct asymptotic power, and the variational principle further ensures that the accuracy improves monotonically with the order of approximation. With these features, accurate results are obtained for larger qa (5% accuracy for qa <= 3), so that results are valid for rods that are by no means thin. Application of these methods to the flexural case has been presented separately. (C) 2011 Elsevier Ltd.

[17]

王礼立

. 应力波基础, 第2版. 北京 : 国防工业出版社 , 2005

[本文引用: 4]

( Wang

Lili

. Foundation of Stress Waves , 2nd Edition. Beijing : National Defense Industry Press , 2005 (in Chinese))

[本文引用: 4]

[18]

Brizard

, Jacquelin

, Ronel

Polynomial mode approximation for longitudinal wave dispersion in circular rods

Journal of Sound and Vibration, 2019 , 439 : 388 - 397

[19]

Naitoh

, Daimaruya

The influence of a rise time of longitudinal impact on the propagation of elastic waves in a bar

Bulletin of JSME, 1985 , 28 : 20 - 25

DOIURL[本文引用: 2]

[20]

Yang

, Shim

VPW

An analysis of stress uniformity in split Hopkinson bar test specimens

International Journal of Impact Engineering, 2005 , 31 ( 2 ): 129 - 150

[21]

Follansbee

, Frantz

Wave propagation in the split Hopkinson pressure bar

Journal of Engineering Materials and Technology, 1983 , 105 : 61 - 66

The initiation and propagation of damage in pure ice specimens under high rate compressive loading at the strain rate range of 100 s^-1to 600 s^-1was studied by means of Split Hopkinson Pressure Bar measurements with incorporated high-speed videography. The results indicate that local cracks in specimens can form and propagate before the macroscopic stress maximum is reached. The estimated crack velocity was in the range of 500 m/s to 1300 m/s, i.e., lower than, but in similar order of magnitude as the elastic wave speed within ice. This gives reason to suspect that already at this strain rate the specimen is not deforming under perfect force equilibrium when the first cracks initiate and propagate. In addition, in contrast to quasi-static experiments, in the high rate experiments the specimens showed notable residual load carrying capacity after the maximum stress. This was related to dynamic effects in fractured ice particles, which allowed the specimen to carry compressive load even in a highly damaged state.

[22]

Gong

, Malvern

, Jenkins

Dispersion investigation in the split Hopkinson pressure bar

Journal of Engineering Materials and Technology, 1990 , 112 : 309 - 314

[23]

, Lambros

Determination of the dynamic response of brittle composites by the use of split Hopkinson pressure bar

Composites Science and Technology, 1999 , 59 : 1097 - 1107

[24]

Park

, Zhou

Separation of elastic waves in split Hopkinson bars using one-point strain measurements

Experimental Mechanics, 1999 , 39 ( 4 ): 287 - 294

[25]

Gama

, Lopatnikov

, Gillespie

Hopkinson bar experimental technique: A critical review

Applied Mechanics Reviews, 2004 , 57 ( 4 ): 223 - 250