水中冲击波与弹性薄板耦合作用的研究

刘德新; 冯洪庆; 佟景伟; 李鸿琦

doi:10.11883/1001-1455(2004)04-0324-6

摘要: 对水中冲击波与弹性薄板的流固耦合作用进行了研究,建立了实验装置。利用全息干涉法定量测量了水中动态流场,并测出受水中冲击波作用的弹性薄板在不同时刻的变形,得出相应的全息干涉条纹图。同时利用有限元法对水中冲击波与弹性薄板的流固耦合问题进行了数值计算。研究结果表明,用实验方法能够较好地模拟水中冲击波与弹性结构的耦合作用,用有限元法计算水中冲击波传播问题是可行的。

在地震或地下爆炸激发的地震波传播研究中, 具有一定压力边界的有限球形空腔可以视为震源, 因此固体介质中的球面应力波传播理论对于研究地震或地下爆炸的震源特征、应力波的传播演化规律等具有重要作用。F.G.Blake^[1]和H.L.Selberg^[2]基于折合位移势(reduced displacement potential, RDP)给出了理想弹性球面波传播的处理方法。RDP是球面波波动方程的一个解, 球面波传播中涉及到的粒子速度、粒子位移、应力、应变等物理量均可通过RDP给出。H.C.Rodean^[3]研究了理想弹性介质中球形发散压缩波的传播, 并讨论了球形发散应力波的特征频率(周期、波长)及其衰减规律。H.C.Rodean的工作对理想弹性球面应力波的传播规律进行了完整的总结。由于波传播和材料的动态力学性能相关, 研究者们基于不同的本构方程(弹塑性^[4-8]、弹黏塑性^[9-11]、Maxwell体^[12]和Kelvin-Voigt体^[13]黏弹性、线性和非线性ZWT^[14-20]等)对球面应力波的传播特征进行了分析, 进一步完善了球面波传播理论。

在考虑介质的黏弹性条件下, 从震源辐射出的波频率特征和理想弹性条件下波的传播特征不同。为深入研究黏弹性球面应力波的频率响应特性随传播距离的变化特征, 本文中基于线黏弹性球面波Laplace域的理论解^[19], 得到不同传播距离处粒子速度、应力、应变等物理量的传递函数, 并以标准线性固体模型为例, 重点讨论粒子速度频率响应函数的传播特征。以弹性半径为0.025 m的空腔爆炸为例, 采用Laplace数值逆变换方法^[21]对粒子速度波形的演化进行分析, 给出粒子速度强间断幅值及粒子速度峰值随传播距离变化的衰减曲线。

1. 线黏弹性球面应力波的频率响应函数

无限介质中线黏弹性球面波的Laplace解可按如下方程计算^[19]：

${{{\tilde{v}}}_{r}}(r, s)=-\left[ \frac{1}{r}\beta (s)+\frac{1}{{{r}^{2}}} \right]s\tilde{\varphi }(r, s)$

(1)

${{{\tilde{u}}}_{r}}(r, s)=-\left[ \frac{1}{r}\beta (s)+\frac{1}{{{r}^{2}}} \right]\tilde{\varphi }(r, s)$

(2)

${{{\tilde{\sigma }}}_{r}}(r, s)=\frac{1}{r}\left[ \left( \tilde{K}(s)+\frac{4}{3}\tilde{G}(s) \right){{\beta }^{2}}(s)+4\tilde{G}(s)\left( \frac{1}{r}\beta (s)+\frac{1}{{{r}^{2}}} \right) \right]\tilde{\varphi }(r, s)$

(3)

${{{\tilde{\sigma }}}_{\theta }}(r, s)=\frac{1}{r}\left[ \left( \tilde{K}(s)+\frac{2}{3}\tilde{G}(s) \right){{\beta }^{2}}(s)-2\tilde{G}(s)\left( \frac{1}{r}\beta (s)+\frac{1}{{{r}^{2}}} \right) \right]\tilde{\varphi }(r, s)$

(4)

${{{\tilde{\varepsilon }}}_{r}}(r, s)=\left[ \frac{1}{r}{{\beta }^{2}}(s)+\frac{2}{{{r}^{2}}}\beta (s)+\frac{2}{{{r}^{3}}} \right]\tilde{\varphi }(r, s)$

(5)

${{{\tilde{\varepsilon }}}_{\theta }}(r, s)=-\left[ \frac{1}{{{r}^{2}}}\beta (s)+\frac{1}{{{r}^{3}}} \right]\tilde{\varphi }(r, s)$

(6)

式中： ${\tilde{K}}$ (s)和 ${\tilde{G}}$ (s)分别为复体积模量和复剪切模量, ${\tilde{v}}$ _r和 ${\tilde{u}}$ _r分别为径向粒子速度和粒子位移, ${\tilde{\sigma }}$ _r和 ${\tilde{\sigma }}$ _θ分别为径向应力和切向应力, ${\tilde{\varepsilon }}$ _r和 ${\tilde{\varepsilon }}$ _θ分别为径向应变和切向应变, s为Laplace变量, β(s)= $\sqrt{\frac{{{\rho }_{0}}{{s}^{2}}}{\tilde{K}(s)+\frac{4}{3}\tilde{G}(s)}}$ 为Laplace域的波传播系数。基于公式(1)~(6), 能够得到不同物理量的传递函数。假设弹性半径为r₀, 则不同半径r(r≥r₀)处的传递函数H_g(r, r₀, s)(下标g代表公式(1)~(6)左端的各物理量)可以写为：

${{H}_{g}}(r, {{r}_{0}}, s)=\frac{g(r, s)}{g({{r}_{0}}, s)}$

(7)

令s=ωi(ω=2πf为圆频率, f为频率, i为虚数单位), 则公式(7)表征不同物理量传递函数对应的频率响应函数。需注意的是, 上面的分析中并没有基于特定的线黏弹性本构模型, 即上述公式具有分析线黏弹性球面发散应力波的普适性。

2. 基于标准线性固体模型的粒子速度频率响应函数传播特征

在地下爆炸或地震研究中, 地运动信号(粒子加速度、速度、位移)是常见的分析对象, 因此本文中仅对粒子速度的频率响应函数进行讨论(理论上粒子加速度、速度、位移的频率响应函数相同)。粒子速度的频率响应函数可以写为：

${{H}_{{{v}_{r}}}}(r, {{r}_{0}}, \omega \text{i)=}\frac{\left[ \frac{1}{r}\beta (\omega \text{i})+\frac{1}{{{r}^{2}}} \right]}{\left[ \frac{1}{{{r}_{0}}}\beta (\omega \text{i})+\frac{1}{r_{0}^{2}} \right]}{{\text{e}}^{-\beta (\omega \text{i})(r-{{r}_{0}})}}$

(8)

线黏弹性球面应力波的传播系数β(ω i)可以写为：

$\beta (\omega \text{i})=\alpha (\omega )+k(\omega )\text{i}$

(9)

式中:α(ω)为线黏弹性球面应力波衰减因子, k(ω)为线黏弹性球面应力波的波数(k(ω)=ω/c(ω), c(ω)为线黏弹性球面应力波的相速度)。把公式(9)代入公式(8), 则有：

${{H}_{{{v}_{r}}}}(r, {{r}_{0}}, \omega \text{i)=}\frac{{{r}_{0}}}{r}\text{=}\frac{\left[ \alpha (\omega )+k(\omega )\text{i}+\frac{1}{r} \right]}{\left[ \alpha (\omega )+k(\omega )\text{i}+\frac{1}{{{r}_{0}}} \right]}{{\text{e}}^{-\alpha (\omega )(r-{{r}_{0}})}}{{\text{e}}^{-k(\omega )\text{i}(r-{{r}_{0}})}}=\left| {{H}_{{{v}_{r}}}}(r, {{r}_{0}}, \omega \text{i)} \right|\\ {{\text{e}}^{\text{i}{{\varphi }_{{{v}_{r}}}}(r, {{r}_{0, }}\omega )}}$

(10)

式中:|H_{v_r}(r, r₀, ωi)和φ_{v_r}(r, r₀, ω)分别为粒子速度响应函数的幅频特性和相频特性, 其表达式可以写为：

$\left| {{H}_{{{v}_{r}}}}(r, {{r}_{0}}, \omega \text{i)} \right|=\frac{{{r}_{0}}}{r}=\frac{\sqrt{(\alpha {{(\omega) +{{r}^{-1}})}^{2}}+{{k}^{2}}(\omega )}}{\sqrt{(\alpha {{(\omega) +r_{0}^{-1})}^{2}}+{{k}^{2}}(\omega )}}{{\text{e}}^{-\alpha (\omega )(r-{{r}_{0}})}}$

(11)

${{\varphi }_{{{v}_{r}}}}(r, {{r}_{0}}, \omega \text{)=}-k(\omega )(r-{{r}_{0}})+\Phi (r, \omega )-\Phi ({{r}_{0}}, \omega )$

(12)

公式(12)中Φ(r, ω)满足下式：

$\text{tan(}\Phi (r, \omega ))=\frac{k(\omega )}{\alpha (\omega )+{{r}^{-1}}}$

(13)

本文中基于标准线性固体模型展开对粒子速度频率响应函数传播特征的讨论。标准线性固体模型由1个线性弹簧和1个Maxwell体并联而成, 如所示。此模型的复体积模量 ${\tilde{K}}$ (s)和复剪切模量 ${\tilde{G}}$ (s)可以写为：

图 1 标准线性固体模型

Figure 1. Standard linear solid model

下载: 全尺寸图片幻灯片

${{{\tilde{K}}}_{s}}=\frac{1}{3(1-2\mu )}\left( {{E}_{0}}+\frac{{{E}_{1}}{{\theta }_{1}}s}{1+{{\theta }_{1}}s} \right)={{K}_{0}}+\frac{{{K}_{1}}{{\theta }_{1}}s}{1+{{\theta }_{1}}s}$

(14)

${{{\tilde{G}}}_{s}}=\frac{1}{2\left( 1+\mu \right)}\left( {{E}_{0}}+\frac{{{E}_{1}}{{\theta }_{1}}s}{1+{{\theta }_{1}}s} \right)={{G}_{0}}+\frac{{{G}_{1}}{{\theta }_{1}}s}{1+{{\theta }_{1}}s}$

(15)

式中：μ为泊松比, E₀和E₁分别为线性弹簧和Maxwell体中弹簧的弹性模量, K₀和K₁分别为对应的体积模量, G₀和G₁分别为对应的剪切模量, θ₁为Maxwell体的松弛因子。令s=ωi=2πfi, 把公式(14)~(15)代入公式(9), 可以给出衰减因子α(f)和波数k(f)的解析式^[20]。由文献[20], 当f→∞时(对应高频状态)：

${{\alpha }_{\infty }}=\underset{f\to \infty }{\mathop{\text{lim}}}\, \alpha (f)=\frac{{{E}_{1}}}{2({{E}_{0}}+{{E}_{1}}){{c}_{\infty }}{{\theta }_{1}}}$

(16)

${{k}_{\infty }}=\underset{f\to \infty }{\mathop{\text{lim}}}\, k(f)=\frac{2\pi f}{{{c}_{\infty }}}$

(17)

式中：c_∞代表高频体波波速, 可以写为：

${{c}_{\infty }}=\sqrt{\frac{(1-\mu )({{E}_{0}}+{{E}_{1}})}{\left( 1-2\mu \right)\left( 1+\mu \right){{\rho }_{0}}}}=\sqrt{\frac{({{K}_{0}}+{{K}_{1}})+\frac{4}{3}({{G}_{0}}+{{G}_{1}})}{{{\rho }_{0}}}}$

(18)

当f→0时(对应低频状态)：

${{\alpha }_{0}}=\underset{f\to 0}{\mathop{\text{lim}}}\, \alpha (f)=\frac{2{{\pi }^{2}}}{~{{c}_{0}}}\ \frac{{{E}_{1}}{{\theta }_{1}}}{{{E}_{0}}}{{f}^{2}}$

(19)

${{k}_{0}}=\underset{f\to 0}{\mathop{\text{lim}}}\, k(f)=\frac{2\pi f}{{{c}_{0}}}$

(20)

式中：c₀代表低频体波波速, 可以写为：

${{c}_{0}}=\sqrt{\frac{(1-\mu ){{E}_{0}}}{\left( 1-2\mu \right)\left( 1+\mu \right){{\rho }_{0}}}}=\sqrt{\frac{{{K}_{0}}+\frac{4}{3}{{G}_{0}}}{{{\rho }_{0}}}}$

(21)

结合公式(11)、(16)~(21), 当f→∞时, 有：

${{\left| {{H}_{{{v}_{r}}}}(r, {{r}_{0}}, 2\pi f\text{i)} \right|}_{f\to \infty }}=\frac{{{r}_{0}}}{r}{{\text{e}}^{-{{\alpha }_{\infty }}(r-{{r}_{0}})}}$

(22)

当f→0时, 有：

${{\left| {{H}_{{{v}_{r}}}}(r, {{r}_{0}}, 2\pi f\text{i)} \right|}_{f\to 0}}=\frac{r_{0}^{2}}{{{r}^{2}}}$

(23)

当公式(22)~(23)满足下式时：

${{\left| {{H}_{{{v}_{r}}}}(r, {{r}_{0}}, 2\pi f\text{i)} \right|}_{f\to \infty }}=\frac{{{r}_{0}}}{r}{{\text{e}}^{-{{\alpha }_{\infty }}(r-{{r}_{0}})}}={{\left| {{H}_{{{v}_{r}}}}(r, {{r}_{0}}, 2\pi f\text{i)} \right|}_{f\to 0}}=\frac{r_{0}^{2}}{{{r}^{2}}}$

(24)

速度响应函数对高频和低频的响应相同, 通过整理公式(24), 有：

$\text{ln}r-{{\alpha }_{\infty }}r=\text{ln}{{r}_{0}}-{{\alpha }_{\infty }}{{r}_{0}}$

(25)

定义函数χ(r), 有：

$\chi (r)=\text{ln}r-{{\alpha }_{\infty }}r$

(26)

如图 2所示, 函数χ(r)在r∈(0, $\alpha _{\infty }^{-1}$ ]区间内是递增函数, 在r∈[ $\alpha _{\infty }^{-1}$ , ∞)区间内是递减函数, 在r= $\alpha _{\infty }^{-1}$ 处取得极大值。存在下述两种情况：(1)当r₀≤ $\alpha _{\infty }^{-1}$ 时, 在r∈[ $\alpha _{\infty }^{-1}$ , ∞)区间内存在r₁满足式(25), 使得高频幅频响应和低频幅频响应相同。在r∈[r₀, r₁]区间内高频幅频响应高于低频幅频响应, 在r∈(r₁, ∞)区间内高频幅频响应低于低频幅度响应；(2)当r₀> $\alpha _{\infty }^{-1}$ 时, 在r∈[r₀, ∞)区间内当且仅当r₁=r₀时满足公式(25), 由此可知在r∈(r₀, ∞)区间内高频幅频响应均低于低频幅频响应。

图 2 函数χ(r)随r的变化

Figure 2. Function χ(r) vs. r

下载: 全尺寸图片幻灯片

以黄土标准线性固体模型参数(如表 1中所示)为例, 此时α_∞≈3.55 m^-1( $\alpha _{\infty }^{-1}$ ≈0.282m)。图 3给出了弹性半径r₀=0.025 m(r₀≤ $\alpha _{\infty }^{-1}$ )时粒子速度频率响应函数随传播距离的变化, 可以看出，图 3中曲线符合上述第一种情况的描述。当r>r₁=r₀> $\alpha _{\infty }^{-1}$ 时, 其粒子速度频率响应函数和图 3中r=2.5 m时的曲线类似, 高频幅频响应均低于低频幅频响应, 在频率f₀(定义为卓越频率)处使得幅频响应取得最大值, 频率f₁(定义为上限频率, 高于此频率的幅频响应均低于低频幅频响应, 即信号幅频响应主要集中在上限频率以下)处的幅频响应等于低频幅频响应。

表 1 黄土材料标准线性固体模型参数

Table 1. Parameters of the standard linear solid model for loess

密度ρ/(kg·m³)	弹性模量E₀/GPa	弹性模量E₁/GPa	松弛时间θ₁/μs	泊松比μ
1 800	1.60	0.33	21.0	0.25

下载: 导出CSV

| 显示表格

图 3 幅频响应函数H_{v_r}(r, r₀, 2πfi)随频率f的变化

Figure 3. Amplitude-frequency response function H_{v_r}(r, r₀, 2πfi) vs. frequency f

下载: 全尺寸图片幻灯片

图 4给出了弹性半径r₀分别为0.025、0.25、2.5、25和250 m时卓越频率f₀和上限频率f₁随传播距离的变化。对于特定场地的地下填实爆炸来说, 弹性半径的大小和爆炸当量相关, 爆炸当量越大则弹性半径越大。从图 4中可以看出, f₀和f₁随着传播距离的增加而减小, 爆炸当量越大, 卓越频率f₀和上限频率f₁越小。为对比近区和远区粒子速度响应函数的频率特征, 表 2给出了不同弹性半径时0.5和1 000 km处的卓越频率f₀和上限频率f₁。弹性半径为0.025~25 m时, 在接近1 000 km传播距离处, 此4条曲线在图 4中基本重合, 表 2中的数据也证实了这一点, 其卓越频率f₀在2.33~2.48 Hz之间, 上限频率f₁在11.32~11.71 Hz之间。可以看出不同弹性半径下, 上限频率f₁的变化率比卓越频率f₀要小。弹性半径为25 m时对应的爆炸当量是弹性半径为0.025 m时爆炸当量的10⁹倍, 而在1 000 km处粒子速度频率响应函数的卓越频率f₀和上限频率f₁基本一致, 这说明当波的传播距离据爆心足够远时爆炸当量对地震信号的频谱特性没有显著影响。图 4中弹性半径为250 m时的频率曲线在1 000 km处卓越频率f₀为1.19 Hz、上限频率f₁为10.09 Hz。从变化趋势上看, 若传播距离继续增加, 弹性半径为250 m时的频率曲线会和弹性半径分别为0.025、0.25、2.5、25 m的4条频率曲线更加接近。

图 4 卓越频率f₀和上限频率f₁随传播距离r的变化

Figure 4. Predominant frequency f₀ and upper limit frequency f₁ varying with propagation distance r

下载: 全尺寸图片幻灯片

表 2 卓越频率f₀和上限频率f₁的变化

Table 2. Variations of predominant frequency f₀ and upper limit frequency f₁

弹性半径r₀/m	0.5 km处的卓越频率f₀/Hz	0.5 km处的上限频率f₁/Hz	1 000 km处的卓越频率f₀/Hz	1 000 km处的上限频率f₁/Hz
0.025	110.76	419.65	2.48	11.71
0.25	109.28	414.41	2.48	11.71
2.5	73.82	361.39	2.47	11.70
25.0	26.92	278.37	2.33	11.32
250.0	9.43	184.49	1.19	10.09

下载: 导出CSV

| 显示表格

3. 粒子速度波形的演化

前面对不同位置粒子速度的频率响应函数特征进行了讨论, 下面就粒子速度的波形演化特征进行数值模拟分析, 主要是研究强间断及粒子速度幅值随传播距离的变化规律。模拟中采用黄土作为示例材料(材料参数如表 1所示), 弹性半径取为0.025 m, 假设空腔边界的压力形式为：

${{\sigma }_{{{r}_{0}}}}={{\sigma }_{r}}({{r}_{0}}, t)=-({{\sigma }_{\text{d}}}{{\text{e}}^{-t/{{t}_{0}}}}+{{\sigma }_{\text{s}}})$

(27)

式中：σ_d+σ_s为空腔边界上施加压力的强间断幅值, σ_s为稳态空腔压力, t₀为压力的延迟时间。取σ_d=σ_s=1 MPa, t₀=21 μs。

图 5给出了0.025~4.775 m区域内不同位置处的粒子速度波形。由公式(18)和表 1, 可以给出c_∞=1 134.3 m/s。由于空腔边界处施加了强间断压力边界条件, 在空腔边界上粒子速度的强间断幅值 $v_{r\_\text{sd}}^{\text{visco}}$ (r₀)可以写为：

图 5 不同位置的粒子速度波形

Figure 5. Particle velocity histories at different locations

下载: 全尺寸图片幻灯片

$v_{r\_\text{sd}}^{\text{visco}}({{r}_{0}})=\frac{{{\sigma }_{\text{d}}}+{{\sigma }_{\text{s}}}}{{{\rho }_{0}}{{c}_{\infty }}}=0.980\ 8\ \text{m}/\text{s}$

(28)

粒子速度强间断幅值 $v_{r\_\text{sd}}^{\text{visco}}$ (r)可以写为^[19]：

$v_{r\_\text{sd}}^{\text{visco}}(r)=\frac{v_{r\_\text{sd}}^{\text{visco}}({{r}_{0}}){{r}_{0}}}{r}{{\text{e}}^{-{{\alpha }_{\infty }}(r-{{r}_{0}})}}$

(29)

通过图 5可以看出, 粒子速度强间断幅值衰减很快, 在图 5(c)中, r≥2.275 m处粒子速度强间断幅值已经几乎观察不到了, 这主要是高频波衰减快造成的。在本文中给出的计算区域内, 含有强间断的粒子速度波形逐渐演化为具有一定上升时间的波形。这里需要指出, 即使r≥2.275 m处的粒子速度强间断从视觉上观察不到, 但理论上是存在的。另外, 通过图 5还可以看出, 当r≥0.775 m时, 粒子速度强间断幅值不再是粒子速度波形的峰值, 这主要是因为不同位置处粒子速度的幅频响应函数对不同频率的波响应不同, 这些位置处高频波的幅频响应不是整个幅频响应曲线的峰值。

图 6给出了黏弹性黄土介质中粒子速度峰值、强间断幅值随传播距离的变化规律。作为对比, 图 6还给出了理想弹性材料中粒子速度峰值随传播距离的变化规律(理想弹性材料中粒子速度强间断幅值和粒子速度峰值相等)。从图 6可以看出, 黏弹性介质中粒子速度幅值的衰减曲线介于理想弹性介质中粒子速度幅值衰减曲线和黏弹性介质中粒子速度强间断幅值衰减曲线之间。在近区, 理想弹性材料中粒子速度强间断幅值衰减曲线(仅有几何扩散效应)和黏弹性介质中粒子速度峰值(或强间断幅值)的衰减曲线的偏离程度较小, 而随着波传播距离的增加, 两者的偏离程度加大。这说明：(1)在近区, 球面波的几何扩散效应对粒子速度峰值(或强间断幅值)的衰减起主导作用；(2)在远区, 介质的黏性对粒子速度峰值(或强间断幅值)的衰减起主导作用。

图 6 粒子速度峰值(强间断幅值)随传播距离r变化

Figure 6. Peak values of particle velocity vs. propagation distance r

下载: 全尺寸图片幻灯片

4. 结论

基于球面发散应力波在Laplace域的理论解, 得到了不同传播距离处粒子速度、应力、应变等物理量的传递函数, 并以标准线性固体模型为例, 重点讨论了粒子速度的频率响应函数的传播特征, 给出了卓越频率和上限频率随传播距离变化的曲线。以弹性半径为0.025 mm的空腔爆炸为例, 采用Laplace数值逆变换方法对粒子速度波形的演化进行了模拟。基于上述分析, 得出如下结论：

(1) 黏弹性球面波的幅频响应函数除了和弹性半径r₀及传播距离r有关之外, 还和黏弹性本构引入的本构衰减有关。理想弹性球面波的幅频响应函数仅和弹性半径r₀及传播距离r有关；

(2) 在线黏弹性条件下, 当弹性半径r₀≤ $\alpha _{\infty }^{-1}$ 时, 随着球面波传播距离的增加, 粒子速度幅频响应函数的高频响应会由高于低频响应状态过渡到低于低频响应状态。当弹性半径r₀> $\alpha _{\infty }^{-1}$ 时, 无论球面波传播距离如何变化, 粒子速度幅频响应函数的高频响应永远低于低频响应；

(3) 黏弹性介质中粒子速度幅值的衰减曲线介于理想弹性介质中粒子速度幅值衰减曲线和黏弹性介质中粒子速度强间断幅值衰减曲线之间。在近区, 球面波的几何扩散效应对粒子速度峰值(或强间断幅值)的衰减起主导作用。在远区, 介质的黏性对粒子速度峰值(或强间断幅值)的衰减起主导作用。