有限振幅声波的传播问题是非线性声学领域里的经典问题。随着大功率声源在医疗超声、无损检测以及国防工业等领域的广泛应用，针对有限振幅声波传播问题的研究显得越来越重要^[1-2]。在黏滞热传导介质中，Burgers方程与Westervelt方程是描述有限振幅声波传播的基本模型，Blackstock^[3-4]针对Burgers方程有过详细的研究，并得到了该方程在不同实际情况下的解析解，IM Hallaj等^[5]基于Westervelt方程用FDTD方法研究了非线性声场的时空特性，张振福等^[6]用数值模拟的方法分析了非线性声场的演变规律，并验证了冲击波的形成过程。

尽管对黏滞热传导介质中有限振幅声波传播问题的研究已经十分详尽，但是考虑弛豫效应作用下有限振幅声波传播问题的研究却很少，Wilson^[7]研究了含气泡多孔介质中的线性声传播，并导出了适用于全频段的弛豫传播模型。而在实际的海水、大气介质中，除了黏滞热传导效应外，还存在化学弛豫、分子弛豫等过程带来的弛豫效应，海水介质中甚至还有硫酸镁、硼酸等物质带来的化学弛豫过程，弛豫效应的存在，一方面会对声波产生吸收作用，另一方面会使介质具备较弱的色散特性，因此在研究有限振幅声波在这类介质中的传播特性时，不应该忽略弛豫效应带来的影响。Polyakova^[8]在Burgers方程的基础上考虑了单个弛豫过程的情况，但是由于忽略了色散效应，在得到的方程中没有出现描述色散的项。为了得到弛豫流体介质中有限振幅声波的传播模型，需要对Burgers方程进行修正，并加入由弛豫过程带来的弛豫修正项，得到扩展的Burgers方程。对于扩展的Burgers方程进行解析求解是困难的，数值模拟方法不失为对其进行研究的好方法。张军^[9]在研究水下强声波脉冲传播时考虑了由硫酸镁带来的弛豫吸收作用，但没有单独研究弛豫色散带来的影响。

本文将建立考虑海水弛豫效应后的有限振幅传播模型——扩展的Burgers方程，在声波频率远小于介质的特征弛豫频率时，结合数值模拟的方法，研究声波传播规律，并对冲击波形成的条件及其演变规律做分析。

1 海水弛豫效应作用下有限振幅声波的传播模型及数值算法 1.1 海水弛豫效应作用下有限振幅声波的传播模型

只考虑二阶非线性效应时，黏滞流体介质中的有限振幅声波的传播模型可以由下述Burgers方程来表示^[10]

$ \frac{\partial v}{\partial x}-\frac{\beta }{{{c}_{0}^{2}}}v\frac{\partial v}{\partial \tau }=\frac{b}{2{{\rho }_{0}}{{c}_{0}^{3}}}~\frac{{{\partial }^{2}}v}{\partial {{\tau }^{2}}} $

(1)

式中：ρ₀为介质的密度，c₀为介质处于平衡态时的声速，v为质点振速，x为传播距离，τ为推迟时间，β为介质的非线性系数，b为声扩散率。

除了经典的黏滞热传导吸收以外，介质中的各种驰豫过程还会引起超吸收效应。在这类介质中，考虑了弛豫修正项的扩展的Burgers方程能准确地描述有限振幅声波的传播，方程的形式如下^[11]：

$ \begin{align} &\ \ \ \frac{\partial v}{\partial x}-\frac{\beta }{c_{0}^{2}}v\frac{\partial v}{\partial \tau }=\frac{b}{2{{\rho }_{0}}c_{0}^{3}}~\frac{{{\partial }^{2}}v}{\partial {{\tau }^{2}}}+ \\ &\frac{1}{2{{c}_{0}}}~\frac{\partial }{\partial \tau }\underset{q}{\mathop{\sum }}\, {{\chi }_{q}}\int_{-\infty }^{\tau }{\frac{\partial v}{\partial \tau ^\prime }\text{exp}~\left( -\frac{\tau -\tau ^\prime }{{{\tau }_{q}}} \right)\text{d}\tau ^\prime }~ \\ \end{align} $

(2)

式中：τ_q是介质中各个弛豫过程对应的特征弛豫时间，χ_q=(c_q∞²-c₀²)/c₀²，c_q∞是各个弛豫过程对应的高频声速。

本文主要研究海水介质中硫酸镁弛豫效应对于有限振幅声波传播的影响，为了得到硫酸镁作用下有限振幅声波的传播模型，需要在式(2)中保留硫酸镁弛豫效应带来的修正项，此时式(2)可以写为

$ \begin{align} &\ \ \frac{\partial v}{\partial x}-\frac{\beta }{c_{0}^{2}}v\frac{\partial v}{\partial \tau }=\frac{b}{2{{\rho }_{0}}c_{0}^{3}}~\frac{{{\partial }^{2}}v}{\partial {{\tau }^{2}}}+ \\ &\frac{\chi }{2{{c}_{0}}}\frac{\partial }{\partial \tau }\int_{-\infty }^{\tau }{\frac{\partial v}{\partial \tau ^\prime }\text{exp}~\left( -\frac{\tau -\tau ^\prime }{{{\tau }_{q}}} \right)\text{d}\tau ^\prime } \\ \end{align} $

(3)

式中：χ=(c_∞²-c₀²)/c₀²为硫酸镁弛豫效应的特征参量。式(3)可以写成下述无量纲化的形式:

$ \begin{align} &\frac{\partial \phi }{\partial \sigma }-\phi \frac{\partial \phi }{\partial y}=\text{ }\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}}\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}+ \\ &\text{ }D\int_{-\infty }^{y}{\frac{\partial \phi }{\partial y'}\text{exp}~\left( -\frac{y-y^\prime }{\omega {{\tau }_{r}}} \right)\text{d}y^\prime } \\ \end{align} $

(4)

式中：?=v/v₀，σ=xωβv₀/c₀²=x/x₀，v₀是参考振速，x₀是理想流体介质中频率为ω初始振速为v₀的有限振幅声波的间断距离，y=ωτ，D=χc₀/(2βv₀)表征了介质的弛豫效应与非线性效应之比。

为了研究硫酸镁弛豫效应对有限振幅声波传播过程带来的影响，需要对式(4)进行定量分析。而式(4)是一个非线性积分-微分方程，通过解析求解的方法对其进行分析是困难的，数值计算不失为一种好的方法。为了方便数值计算，对积分项进行化简，当ωτ_r$\ll$ 1时，即声波频率远小于介质的特征弛豫频率时，利用泰勒展开对$\frac{\partial \phi }{\partial y'}$进行分解，并保留起主要作用的前两项，可以得到

$ \frac{\partial \phi }{\partial y'}\approx \frac{\partial \phi }{\partial y}+\text{ }\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}\left( y^\prime -y \right) $

(5)

式(5)代入式(4)可得

$ \frac{\partial \phi }{\partial \sigma }-\phi \frac{\partial \phi }{\partial y}=\text{ }\left( D\omega {{\tau }_{r}}+\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}} \right)\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}-D{{\omega }^{2}}\tau _{r}^{2}\frac{{{\partial }^{3}}\phi }{\partial {{y}^{3}}} $

(6)

式中：$\frac{\partial \phi }{\partial \sigma }$为D′Alembertian项，$\phi \frac{\partial \phi }{\partial y}$为二阶非线性项，$\left( D\omega {{\tau }_{r}}+\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}} \right)\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}$为声吸收项，其中, $D\omega {{\tau }_{r}}\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}$是由弛豫效应带来的吸收项，$\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}}\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}$则是由黏滞热传导效应带来的吸收项，在满足低频近似条件ωτ_r$\ll$1时，由弛豫吸收引起的色散主要体现在$D{{\omega }^{2}}\tau _{_{r}}^{2}\frac{{{\partial }^{3}}\phi }{\partial {{y}^{3}}}$项。对于无量纲化的Burgers方程，声雷诺数$\mathit{\Gamma }\text{=}\frac{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}}{b\omega }$可以表征非线性效应与黏滞热传导效应的相对大小。在弛豫介质中，由于介质对声波的耗散作用由弛豫吸收与黏滞吸收两部分组成，所以在这种介质中表征非线性效应与耗散效应相对大小的物理量声雷诺数的表达式应为Γ′=(Dωτ_r+Γ^-1)^-1。

1.2 有限振幅声波传播模型的数值算法

为了结合数值模拟的方法对式(6)进行具体分析，需要在空间域上与时间域上对其离散化，本文主要采用有限差分法求解上述式^[12]。

由于对流项$\phi \frac{\partial \phi }{\partial y}$在数值计算中的处理比较复杂，为了方便，引入一个变量F，令F=ϕ²/2，那么式(6)可以写成不含对流项的形式：

$ \frac{\partial \phi }{\partial \sigma }-\frac{\partial F}{\partial y}=\text{ }\left( D\omega {{\tau }_{r}}+\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}} \right)\frac{{{\partial }^{2}}\phi }{\partial {{y}^{2}}}-D{{\omega }^{2}}\tau _{r}^{2}\frac{{{\partial }^{3}}\phi }{\partial {{y}^{3}}} $

(7)

首先引入以下记号表示各阶差分算子：

$ F_{j}^{n}\approx \frac{{\rm{1}}}{{\rm{2}}}(F_{j}^{n}+F_{j\text{+}1}^{n})=\text{ }\frac{{\rm{1}}}{{\rm{2}}}\phi _{j}^{n}\phi _{j\text{+}1}^{n} $

(8)

$ {{\delta }_{\sigma }}\phi _{j}^{n}=\frac{\phi _{j\text{+}1}^{n}-\phi _{j}^{n}}{\Delta \sigma } $

(9)

$ {{\delta }_{y}}F_{j}^{n}=\frac{F_{j}^{n+1}-F_{j}^{n-1}}{2\Delta y} $

(10)

$ \delta _{y}^{2}\phi _{j}^{n}=\frac{\phi _{j}^{n+1}-2\phi _{j}^{n}+\phi _{j}^{n-1}}{{{(\Delta y)}^{2}}} $

(11)

$ \delta _{y}^{3}\phi _{j}^{n}=\frac{\phi _{j}^{n+2}-2\phi _{j}^{n+1}+2\phi _{j}^{n-1}-\phi _{j}^{n-2}}{2{{(\Delta y)}^{3}}} $

(12)

式中：上角标n表示波形中的第n个时间离散点，下角标j表示波形中的第j个空间离散点。将式(7)在ϕ_jⁿ及F_jⁿ点处通过有限差分法进行展开，并结合引入的各阶差分算子，可以得到时间—空间域中任意一点(n, j)处的差分格式，具体形式如下

$ \begin{align} &\ \ \ \ \ \ \ \ \ \ \ {{\delta }_{\sigma }}\phi _{j}^{n}-\frac{1}{2}{{\delta }_{y}}F_{j}^{n}= \\ &\text{ }(D\omega {{\tau }_{r}}+\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}})\delta _{y}^{2}\phi _{j}^{n}-D{{\omega }^{2}}\tau _{r}^{2}\delta _{y}^{2}\phi _{j}^{n} \\ \end{align} $

(13)

对式(13)做适当的变换可以得到

$ \begin{align} &{{a}_{1}}{{\phi }_{j+1}^{n-2}}+{{a}_{2}}{{\phi }_{j+1}^{n-1}}+{{a}_{3}}{{\phi }_{j+1}^{n}}+{{a}_{4}}{{\phi }_{j+1}^{n+1}}+{{a}_{5}}{{\phi }_{j+1}^{n+2}}= \\ &\text{ }{{b}_{1}}{{\phi }_{j}^{n-2}}+{{b}_{2}}{{\phi }_{j}^{n-1}}+{{b}_{3}}{{\phi }_{j}^{n}}+{{b}_{4}}{{\phi }_{j}^{n+1}}+{{b}_{5}}{{\phi }_{j}^{n+2}} \\ \end{align} $

(14)

其中

$ \begin{align} &{{a}_{1}}=-\frac{1}{4}Bq, \text{ }{{a}_{2}}=\frac{1}{4}s\varphi _{j}^{n-1}-\text{ }\frac{1}{2}Ar+\frac{1}{2}Bq, \\ &\text{ }{{a}_{3}}=1+Ar, {{a}_{4}}=-\frac{1}{4}s\varphi _{j}^{n+1}-\text{ }\frac{1}{2}Ar-\frac{1}{2}Bq, \\ &{{a}_{5}}=\frac{1}{4}Bq, \text{ }{{b}_{1}}=\frac{1}{4}Bq, \text{ }{{b}_{2}}=\frac{1}{2}Ar-\frac{1}{2}Bq, \text{ }{{b}_{3}}=1-Ar, \\ &{{b}_{4}}=\frac{1}{2}Ar+\text{ }\frac{1}{2}Bq, {{b}_{5}}=-\frac{1}{4}Aq, A=D\omega {{\tau }_{r}}+\frac{b\omega }{2{{\rho }_{0}}\beta {{v}_{0}}{{c}_{0}}}, \text{ } \\ &B=D{{\omega }^{2}}\tau _{r}^{2}, s=\frac{\Delta \sigma }{\Delta y}, r=\frac{\Delta \sigma }{{{(\Delta y)}^{2}}}, \text{ }q=\frac{\Delta \sigma }{{{(\Delta y)}^{3}}} \\ \end{align} $

式(14)表示在差分算法中空间第j+1层与空间第j层的关系，在基于(n, j)点的差分方程中，用到了(n-2, j)，(n-1, j)，(n+1, j)，(n+2, j)，(n-2, j+1)，(n-1, j+1)，(n, j+1)，(n-1, j+1)，(n+2, j+1)这九个点，具体关系如图 1所示。

这是一个依赖于10个点的双层差分格式，如果给出了空间第j层的值，那么根据式(14)通过求解线性方程组可以求出空间第j+1层的值，由此可见这是一个隐式的差分格式。根据式(14)，如果已知空间初始位置j=0处的值，则可以计算出空间任意位置处的值。

式(14)同时可以写成如下矩阵的形式:

$ {\mathit{\boldsymbol{A}}_n}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{{\rm{ }}_{j + 1}} = {\mathit{\boldsymbol{B}}_n}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{{\rm{ }}_j} + {\mathit{\boldsymbol{E}}_n} $

(15)

其中向量Φ_j表示空间第j层的信号波形

$ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}}_{j}}=\left[\begin{matrix} \phi _{j}^{1} \\ \phi _{j}^{1} \\ \vdots \\ \phi _{j}^{^{N-1}} \\ \phi _{j}^{N} \\ \end{matrix} \right] $

(16)

式中:N表示时间采样点的总数，在边界条件上，令ϕ_j+1¹=ϕ_j¹=0，ϕ_j+1²=ϕ_j²=0，ϕ_j+1^N-1=ϕ_j^N-1=0以及ϕ_j+1^N=ϕ_j^N=0，这样在已知空间初始位置的值时，可以通过矩阵分解法准确地求出解向量Φ_j的值，从而能够有效地预测有限振幅声波波形的演变规律。

2 海水弛豫效应作用下有限振幅声波传播的仿真分析

结合式(15)给出的空间递推关系及给定的空间初始位置处声信号波形来模拟海水介质中有限振幅声波的传播规律，通过数值计算可以得到仿真结果。

根据MA Ainslie和JG McColm的经验公式，10 ℃海水介质中硫酸镁的弛豫频率为75.6 kHz^[13]，海水的非线性系数为4.99^[14]。若声源在初始位置σ=0处发出频率为7 kHz，声压幅值为200 kPa的单频声信号，无量纲参量取χ=0.01，Γ=(bω/2ρ₀βv₀c₀)^-1=1^[9]。在仿真计算中取时间波形的采样频率为100π kHz，空间步长Δσ=0.002 2。

图 2给出的四组对比图分别描述了间断距离σ=1以内四个不同距离处的声波波形，每一组对比图中给出的分别是同一时间段内三种不同情况下的声波波形。从四组对比图可以看出，在间断距离以内，黏滞热传导效应对声波有吸收作用，由此导致声波幅值在传播过程中逐渐衰减，而弛豫效应会使这一现象加剧，这是因为弛豫效应会使得介质对声波产生逾量吸收作用，这与客观事实是相符的^[15]，同时从式(6)中也可以看出弛豫效应对于介质的声吸收作用是有贡献的。在间断距离以内，理想介质中的声波波形随着传播距离的增加逐渐向冲击波演变，从图 2(d)中可知当σ=1时波形的一侧已经变得非常陡峭，而黏滞热传导效应与弛豫吸收效应具有推迟这一过程形成的作用，这是因为理想介质中声波非线性传播会产生高次谐波，在时域波形上体现为所有的高次谐波都叠加在基波上，而黏滞热传导效应及弛豫吸收效应使得介质对于基波以及各阶谐波都会产生吸收作用，且介质的吸收作用随着频率的提高更为明显，高次谐波的能量在传播过程中很快损耗殆尽，所以声波的时域波形也就更接近基波对应的时域波形。

弛豫效应的存在，一方面会对声波产生逾量吸收，另一方面会使介质具备较弱的色散特性，这不仅从式(6)中可以明显看出，而且根据Kramers-Kronig色散关系可知^[16]，吸收介质原则上是色散介质，所以存在吸收的同时必定存在色散，弛豫介质是弱色散介质，所以弛豫效应产生吸收的同时也带来了色散。为了单独研究由于弛豫效应带来的色散对声波传播的影响，在仿真分析中，不考虑声吸收项，得到的结果如图 3所示。通过两者的对比可知在间断距离内，弛豫效应带来的色散会使声波波形变得不对称，不对称性随着传播距离的增加而加剧。文献[14]的研究结果也体现了上述结论。

在理想介质中，随着传播距离的增加，有限振幅声波的波形逐渐发生畸变，冲击波的形成发生在间断距离σ=1处，而实际的海水介质中同时存在黏滞吸收与弛豫吸收效应，这两种效应均具有推迟冲击波形成的作用。在这种介质中冲击波形成的阈值条件为Γ~1^[17]，当声雷诺数超过这个阈值时，很有可能形成冲击波。

为了研究海水介质中冲击波的形成问题，不考虑由硫酸镁带来的弛豫色散效应，则在这种介质中冲击波形成的阈值条件应为Γ′~1。仿真结果给出了当Γ′取不同值时，声波传播过程中的波形变化图。为了方便对比，在图 4中加入了理想介质间断距离σ=1处的冲击波波形，从图 4可以看出，当Γ′=20时，声波在σ=1时处于非线性积累阶段，此时的声波传播处于第一阶段，当声波传播至σ=π/2处附近时的波形与虚线的重合度非常高，并且冲击波的幅值达到最大值，出现这种现象的原因是当Γ$\gg $1时，在离声源不太远的距离内，非线性效应对声波传播的作用远大于耗散效应的作用，通过非线性的积累过程可以逐渐形成冲击波。当声波传播至σ=π/2处时冲击波具有最大的幅值，由于介质耗散效应的存在，此时冲击波的幅值略小于理想介质中冲击波的幅值，且继续传播的声波其幅值将逐渐减小，并逐渐从第二阶段过渡到第三阶段。从图 5中可以看出，当Γ′没有达到冲击波形成的阈值条件时，在传播过程中声波的幅值迅速衰减，同时声波波形几乎保持了初始的正弦波形，在离声源不远处声波已经衰减至小振幅波。这是因为当Γ$\ll $1时，非线性效应对声波传播的作用远小于耗散效应的作用，声波非线性传播产生的各阶谐波在很短的距离内便损耗殆尽了。

3 结论

1) 黏滞热传导效应对于声波的能量有吸收作用，使得声波幅值在传播过程中逐渐衰减，由硫酸镁弛豫效应引起的逾量吸收作用会使这一现象加剧。

2) 理想介质中的声波波形随着传播距离的增加逐渐向冲击波演变，硫酸镁弛豫吸收效应具有推迟冲击波形成的作用。

3) 硫酸镁弛豫效应带来的色散会使得声波波形变得不对称，不对称性随着传播距离的增加而显著。

4) 海水介质中冲击波的形成具有一定的阈值条件，当介质中的声雷诺数超过冲击波形成的阈值条件时，会形成冲击波，冲击波的幅度随着传播距离的增大而衰减，而当介质中的声雷诺数没有达到阈值条件时不会形成冲击波。