1.   引　言

非线性声学是声学的一个重要分支, 当系统的声强比较强时就会产生各种非线性效应, 如谐波产生、冲击波形成、声辐射力的出现等. 非线性声学在某些领域有着重要的运用, 如高声强聚焦超声^[1-4]、超声悬浮^[5-7]、声空化^[8]、声谐波成像^[9,10]、参量发射阵^[11,12]等. 在这些领域中声波的非线性方程的求解是非常重要的.

关于非线性声学系统波动方程的求解, 现阶段用到的方法一般包括: 1)完全的数值计算方法, 如有限元和有限差分法^[13-17]. 这类方法得到解的物理意义并不明确, 很难揭示非线性作用的物理本质, 而且在很多情况下还会引起数值发散问题, 并非适用于全部的非线性问题. 2)严格的解析方法^[18-21]. 这种方法只能处理极少数系统的非线性声学问题, 如理想流体中非线性声波的传播. 3)微扰法^[22]. 它的优点是方法简单和解的物理意义清晰, 但是只适合处理低声强时的非线性效应. 且它只考虑低阶谐波对高阶谐波的作用, 而忽略其反作用, 因此并不满足能量守恒定律.

对于声学非线性方程的求解问题, 本文提出了一种新的、半解析的数值迭代方法. 它是在频域内把声场展开为傅里叶级数的形式, 实现时间变量和空间坐标的分离. 然后根据计算精度的具体需求, 截断高频谐波而实现方程的求解. 它的解具有非常清晰的物理意义, 即是各阶谐波的组合. 经过研究发现, 在微扰法适用的声强范围内, 本文提出的方法也是适用, 且满足能量守恒定律(无耗散的系统). 在微扰法不适用的一个较宽的声强范围内, 迭代法依然适用且满足能量守恒定律(无耗散的系统). 只是在极高声强的情况下, 本文提出的方法才不适用.

2.   理论方法

2.1   非线性波动方程的迭代数值方法

在拉格朗日坐标系下, 黏性液体中一维非线性声波的位移满足下式^[23]:

$$\begin{split} & - \frac{{\rm{1}}}{{c_0^2}}{u_{,tt}} + {u_{,xx}} + \alpha {u_{,xxt}} \\=\; & 2\beta {u_{,x}}{u_{,xx}} + \alpha ({u_{,xxt}}{u_{,x}} + {u_{,xt}}{u_{,xx}}), \end{split}$$

(1)

式中$\beta$为液体的非线性系数; ${c_0}$是静态时(即不存在声波时)液体的声速; $\alpha = {{3\mu }}/({{{\rho _0}c_0^2}})$, 其中$\mu$为液体的体积黏滞系数, ${\rho _0}$是静态时(即不存在声波时)液体的密度. 下标中逗号后的坐标x和时间t表示对它们求偏导数. 令$\mu = 0$, 则方程(1)退化为理想液体中一维非线性声波的波动方程.

在许多情况下(如求解稳态问题), 把u中的变量t和x分离开来是有利的, 一般情况下u可以表示为(可以称为频域内的傅里叶级数展开):

$$\begin{split} u =\; & \dfrac{A_0}2 + \sum\limits_{n = 1}^\infty {\left( {1/2} \right){A_n}\exp \left( {{\rm{i}}n\omega t} \right)} \\ &+ \sum\limits_{n = 1}^\infty {\left( {1/2} \right)A_n^*\exp \left( { - {\rm{i}}n\omega t} \right)} , \end{split}$$

(2)

其中i是一个虚单位, $\omega$表示波的角频率; ${A_0}$为一个实数场变量, ${A_0}/2$表示声波的“直流”部分; ${A_n}(n \geqslant 1)$为第n阶谐波的复数场变量(即复振幅), ${A_n}\exp \left( {{\rm{i}}n\omega t} \right)$的实部是第n阶谐波真实的位移, $A_n^*(n \geqslant 1)$是${A_n}$的复数共轭场变量. 注意${A_n}$和$A_n^*$中已经不包含时间变量t, 它们只是空间坐标x的函数.

一般情况下高阶谐波是比较弱的, 根据计算精度的需要可以忽略掉某些高阶谐波. 为了简化理论的叙述, 本文只考虑阶数小于或者等于$N(N \leqslant 6)$的各阶谐波, 忽略掉其他的高阶谐波, 称之为N阶近似.

将方程(2)代入方程(1)中, 因为有相同时间因子($\exp \left( {{\rm{i}}n\omega t} \right)$, $n = 0, \pm 1, \pm 2, \cdots$)的项之和必须为零, 所以可以得到下面的方程:

$${A_{0,xx}} = \beta {F_0},$$

(3)

$$(1 + {\rm{i}}\omega \alpha ){A_{1,xx}} + {\omega ^2}c_0^{ - 2}{A_1} = \Big(\beta + \frac{1}{2}{\rm{i}}\omega \alpha \Big){F_1},$$

(4)

$$(1 + 2{\rm{i}}\omega \alpha ){A_{2,xx}} + 4{\omega ^2}c_0^{ - 2}{A_2} = \Big(\beta + {\rm{i}}\omega \alpha \Big){F_2},$$

(5)

$$(1 + 3{\rm{i}}\omega \alpha ){A_{3,xx}} + 9{\omega ^2}c_0^{ - 2}{A_3} =\Big (\beta + \frac{3}{2}{\rm{i}}\omega \alpha \Big){F_3},$$

(6)

$$(1 + 4{\rm{i}}\omega \alpha ){A_{4,xx}} + 16{\omega ^2}c_0^{ - 2}{A_4} =\Big (\beta + 2{\rm{i}}\omega \alpha \Big){F_4},$$

(7)

$$(1 + 5{\rm{i}}\omega \alpha ){A_{5,xx}} + 25{\omega ^2}c_0^{ - 2}{A_5} =\Big (\beta + \frac{5}{2}{\rm{i}}\omega \alpha \Big){F_5},$$

(8)

$$(1 + 6{\rm{i}}\omega \alpha ){A_{6,xx}} + 36{\omega ^2}c_0^{ - 2}{A_6} =\Big(\beta + \frac{6}{2}{\rm{i}}\omega \alpha \Big){F_6},$$

(9)

其中

$$\begin{split} {F_0} =\; & {A_{0,x}}{A_{0,xx}} + A_{1,x}^*{A_{1,xx}} + A_{2,x}^*{A_{2,xx}}\\ &+ A_{3,x}^*{A_{3,xx}} + A_{4,x}^*{A_{4,xx}} + A_{5,x}^*{A_{5,xx}} \\ & + A_{6,x}^*{A_{6,xx}} + {A_{1,x}}A_{1,xx}^* + {A_{2,x}}A_{2,xx}^*\\ & + {A_{3,x}}A_{3,xx}^* + {A_{4,x}}A_{4,xx}^* + {A_{5,x}}A_{5,xx}^*\\ & + {A_{6,x}}A_{6,xx}^*, \end{split}$$

(10)

$$\begin{split} {F_1} =\; & {A_{0,x}}{A_{1,xx}} + {A_{1,x}}{A_{0,xx}} + A_{1,x}^*{A_{2,xx}} + {A_{2,x}}A_{1,xx}^*\\ & + A_{2,x}^*{A_{3,xx}} + {A_{3,x}}A_{2,xx}^* + A_{3,x}^*{A_{4,xx}} \\ & + {A_{4,x}}A_{3,xx}^* + A_{4,x}^*{A_{5,xx}} + A_{5,x}^*{A_{6,xx}} \\ &+ {A_{5,x}}A_{4,xx}^* + {A_{6,x}}A_{5,xx}^*,\\[-10pt] \end{split}$$

(11)

$$\begin{split} {F_2} = \; &{A_{1,x}}{A_{1,xx}} + {A_{2,x}}{A_{0,xx}} + {A_{0,x}}{A_{2,xx}} \\ &+ {A_{3,x}}A_{1,xx}^* + A_{1,x}^*{A_{3,xx}} + {A_{4,x}}A_{2,xx}^* \\ &+ A_{2,x}^*{A_{4,xx}} + {A_{5,x}}A_{3,xx}^* + {A_{6,x}}A_{4,xx}^* \\ &+ A_{3,x}^*{A_{5,xx}}+ A_{4,x}^*{A_{6,xx}}, \end{split}$$

(12)

$$\begin{split} {F_3} =\; & {A_{0,x}}{A_{3,xx}} + {A_{3,x}}{A_{0,xx}} + {A_{1,x}}{A_{2,xx}}\\ & + {A_{2,x}}{A_{1,xx}} + A_{1,x}^*{A_{4,xx}} + {A_{4,x}}A_{1,xx}^* \\ & + A_{2,x}^*{A_{5,xx}} + A_{3,x}^*{A_{6,xx}} + {A_{5,x}}A_{2,xx}^* \\ &+ {A_{6,x}}A_{3,xx}^*, \end{split}$$

(13)

$$\begin{split} {F_4} = \; & {A_{0,x}}{A_{4,xx}} + {A_{4,x}}{A_{0,xx}} + {A_{1,x}}{A_{3,xx}} \\ &+ {A_{3,x}}{A_{1,xx}} + {A_{2,x}}{A_{2,xx}} + A_{1,x}^*A{}_{5,xx}\\ & + A_{2,x}^*A{}_{6,xx} + {A_{6,x}}A_{2,xx}^*, \end{split}$$

(14)

$$\begin{split} {F_5} =\; & {A_{0,x}}{A_{5,xx}} + {A_{5,x}}{A_{0,xx}} + {A_{1,x}}{A_{4,xx}} \\ &+ {A_{4,x}}{A_{1,xx}} + {A_{2,x}}{A_{3,xx}} + {A_{3,x}}{A_{2,xx}} \\ &+ {A_{6,xx}}A_{1,x}^* + {A_{6,x}}A_{1,xx}^*, \end{split}$$

(15)

$$\begin{split}{F_6} = \; &{A_{0,x{{x}}}}{A_{6,x}} + {A_{1,xx}}{A_{5,x}} + {A_{2,xx}}{A_{4,x}} + {A_{3,xx}}{A_{3,x}}\\ & + {A_{4,xx}}{A_{2,x}} + {A_{5,xx}}{A_{1,x}} + {A_{6,xx}}{A_{0,x}}.\\[-10pt] \end{split}$$

(16)

注意此处只给出了场变量的方程, 共轭场的方程并没有列出来, 只要对方程(3)—(9)取复数共轭就可以得到共轭场的方程, 因此方程(3)—(9)是完备的.

方程(3)—(9)是一组耦合的非线性方程, 直接求解它们是很困难的. 本文提出求解它们的一种新的简单迭代方法. 用${A^{\left( m \right)}}$和${A^{*\left( m \right)}}$($m \geqslant 0$)表示第m次迭代计算得到的场量. 在第m次迭代计算中, 采用了如下方法: 方程(3)—(9)等号左边的场量取为${A^{\left( m \right)}}$, 右边的场量取为${A^{(m - 1)}}$和${A^{*(m - 1)}}$. 在第m次迭代计算中, 用到如下的方程:

$$A_{0,xx}^{(m)} = \beta F_0^{(m - 1)},$$

(17)

$$(1 + {\rm{i}}\omega \alpha )A_{1,xx}^{(m)} + {\omega ^2}c_0^{ - 2}A_1^{(m)} \!=\!\Big (\beta + \frac{1}{2}{\rm{i}}\omega \alpha \Big)F_1^{(m - 1)},$$

(18)

$$(1 + 2{\rm{i}}\omega \alpha )A_{2,xx}^{(m)} + 4{\omega ^2}c_0^{ - 2}A_2^{(m)} = (\beta + {\rm{i}}\omega \alpha )F_2^{(m - 1)},$$

(19)

$$(1 + 3{\rm{i}}\omega \alpha )A_{3,xx}^{(m)} + 9{\omega ^2}c_0^{ - 2}A_3^{(m)} = \Big(\beta \!+\! \frac{3}{2}{\rm{i}}\omega \alpha \Big)F_3^{(m - 1)},$$

(20)

$$(1 + 4{\rm{i}}\omega \alpha )A_{4,xx}^{(m)} + 16{\omega ^2}c_0^{ - 2}A_4^{(m)} = (\beta + 2{\rm{i}}\omega \alpha )F_4^{\left( {m - 1} \right)},$$

(21)

$$(1 + 5{\rm{i}}\omega \alpha )A_{5,xx}^{(m)} + 25{\omega ^2}c_0^{ - 2}A_5^{(m)} = \Big(\beta \!+\! \frac{5}{2}{\rm{i}}\omega \alpha \Big)F_5^{(m - 1)},$$

(22)

$$(1 + 6{\rm{i}}\mu \omega )A_{6,xx}^{(m)} + 36{\omega ^2}c_0^{ - 2}A_6^{(m)} = (\beta + 3{\rm{i}}\omega \alpha )F_6^{(m - 1)},$$

(23)

分别用${A^{\left( {m - 1} \right)}}$和${A^{*\left( {m - 1} \right)}}$替换F中的A和${A^*}$, 得到的结果就是${F^{\left( {m - 1} \right)}}$. 方程(17)—(23)是一组非耦合的方程, 因此可以分别独立地计算出$A_n^{\left( m \right)}$. 这意味着当涉及到更多的高阶谐波时, 计算量不会急剧地增加.

2.2   本文迭代方法的具体运用过程

用迭代方法研究非线性声波在黏性液体中的传播问题. 现设在$x = 0$处有一列平面波朝x正向传播, 其为入射声波, 且声场可以表示为

$${u_{\rm{i}}} = {B_{\rm{i}}}\exp \left[ {{\rm{i}}\omega \left( {t - x/{c_0}} \right)} \right],$$

(24)

其中${B_{\rm{i}}}$是一个已知量. 入射波的能流密度${P_{\rm{i}}}$可以表示为

$${P_{\rm{i}}} = \frac{1}{2}{\omega ^2}{\rho _0}{c_0}{B_{\rm{i}}}B_{\rm{i}}^{\rm{*}}.$$

(25)

本文需要计算在$x = L$处出射的各阶谐波的声场. 它们可以表示为

$${u_{{\rm{on}}}} \!=\! {B_{{\rm{on}}}}\exp [{{\rm{i}}n\omega ( {t \!-\! x/{c_0}} )} ]\;( n \!=\! 1,2, \cdots ,6),$$

(26)

其中${B_{{\rm{on}}}}$就是要计算的量. 出射波的能流密度${P_{{\rm{on}}}}$可以表示为

$${P_{{\rm{on}}}} = \frac{1}{2}{n^2}{\omega ^2}{\rho _0}{c_0}{B_{{\rm{on}}}}B_{{\rm{on}}}^*.$$

(27)

为了求解在$\left[ {0, L} \right]$坐标间隔内的非线性声场, 用有限差分法来求解方程(17)—(23). 在迭代计算中, 令${A^{\left( 0 \right)}} = 0$和${A^{*\left( 0 \right)}} = 0$, 非零的${A^{\left( m \right)}}$和${A^{*\left( m \right)}}$由边界激励条件产生. 用到的边界条件是: 在$x = 0$和$x = L$两个端点处, 各阶谐波的位移和垂直应力都是连续的.

3.   数值计算和讨论

通过数值计算分析了非线性声波在液体中的传播性质, 得到了本文提出的新数值方法的适用范围, 并证明了其有效性. 在所有的计算中, 如果没有特别说明, 那么用到的参数是: c₀ = 1.5 × 10³ m/s, ${\rho _0} =$ 10³ kg/m³, P₀ = 1.01 × 10⁵ Pa, $\mu = 1$ × 10^–3 Pa/s, $L = 0.05$ m, $\omega =$ 5 × 10⁶ rad/s, $\beta = 3.5$. 如果文中或图中对某个参数有特别说明, 那么该参数就替换为特殊说明处的数据.

图1显示了非线性声波在理想液体($\mu = 0$)中传播时能量守恒的破坏程度${E_{\rm{d}}}$、二阶谐波的相对能流${P_2}/{P_{\rm{i}}}$和三次谐波的相对能流${P_3}/{P_{\rm{i}}}$随入射能流${P_{\rm{i}}}$的变化情况. 能量守恒的破坏程度${E_{\rm{d}}}$定义为所有出射能流密度与入射能流密度的相对差值, 即${E_{\rm{d}}} = \left| {\displaystyle\sum\nolimits_{n = 1}^N {{P_{{\rm{on}}}}} - {P_{\rm{i}}}} \right|/{P_{\rm{i}}}$.${E_{\rm{d}}}$越小能量守恒越能保证, 它越大能量守恒定律破坏程度越大. 实线对应着本文提出的迭代法的情况, 虚线对应着微扰法的情况.

由图1(a)可见, 能量守恒的破坏程度${E_{\rm{d}}}$随入射波能流${P_{\rm{i}}}$的增大而增大. 其他参数不变, 入射波能流相同时, 明显可以看出本文提出的新数值迭代方法得到的能量守恒的破坏程度比微扰法得到的能量守恒的破坏程度小很多. 当入射波能流是$1.5 \times {10^7}$ J/(m²·s)时, 迭代法的破坏程度是7.7 × 10^–3, 微扰法的破坏程度是$0.506$; 当入射能流是$5.5 \times {10^7}$ J/(m²·s)时, 迭代法的破坏程度1.92 × 10^–2, 微扰法的破坏程度是$2.268$.

由图1(b)和图1(c)可知, 二次谐波的相对能流${P_2}/{P_{\rm{i}}}$和三次谐波的相对能流${P_3}/{P_{\rm{i}}}$均随入射能流${P_{\rm{i}}}$的增加而增加. 当入射波能流小于$1.5 \times$ 10⁷ J/(m²·s)时, 迭代法和微扰法得到的${P_2}/{P_i}$和${P_3}/{P_{\rm{i}}}$的值几乎相同; 当入射波能流大于$5.5 \times$ 10⁷ J/(m²·s)时, 迭代法和微扰法得到的${P_2}/{P_{\rm{i}}}$的值分别等于$0.6295$和$1.4741$及${P_3}/{P_{\rm{i}}}$的值分别大于$0.0421$和$0.4847$.

从图1还可以看出, 当入射波能流小于$~1.5 \times {~}$10⁷ J/(m²·s)时, 两种方法(本文提出的新数值迭代法和微扰法)均适用; 当入射波能流小于$5.5 \times$ 10⁷ J/(m²·s)且大于$1.5\; \times$10⁷ J/(m²·s)时, 迭代法适用, 微扰法不适用; 当入射波能流大于$5.5 \times {}$10⁷ J/(m²·s)时, 两种方法均不适用.

图2显示了声波在理想液体($\mu = 0$)中传播时各阶谐波的相对能流${{{P_n}}}/{{{P_{\rm{i}}}}}~(n = 1, 2, \cdots , 6)$随入射波能流${P_{\rm{i}}}$的变化情况. 从图2可见, 各阶谐波的能流随入射波能流的变化趋势相同, 均随入射波能流的增加而增加, 但基频波的能流明显大于其他高阶谐波的能流. 从图2还可以看出, 基频波的相对能流随入射波能流的增加而减小, 高阶谐波的相对能流随入射波能流的增加而增加.

图2显示的结果是与文献[20, 24]中基频波能流与各高阶谐波能流之间的关系相一致. 从图2可以明显得到, 声波在介质中传播时, 基频波能流向各高阶谐波传递, 这同时也解释了, 基频波的相对能流随入射波能流的增加而减少, 但各高阶谐波的相对能流却随入射波能流的增加而增加.

图3显示了声波在理想液体($\mu = 0$)中传播时入射波能流不同的情况下二阶谐波的相对能流${P_2}/{P_{\rm{i}}}$随迭代次数m的变化关系. m表示迭代次数. 由图3可知, 当入射波能流是$0.4687 \times {}$10⁷ J/(m²·s), 迭代次数大于等于3时, 二次谐波的相对能流收敛于$0.145$; 当入射波能流是$2.7 \times {}$ 10⁷ J/(m²·s), 迭代次数大于等于4时, 二次谐波的相对能流收敛于$0.515$; 当入射波能流是$5.418 \times { }$10⁷ J/(m²·s), 迭代次数大于等于7时, 二次谐波的相对能流收敛于$0.639$. 由图3可见, 本文提出的新数值迭代方法具有收敛性, 且入射声强越大, 达到收敛的迭代次数需要越多.

图4显示了声波在黏性液体($\mu \ne 0$)中传播时, 黏度不同的情况下各阶谐波的能流P_n (n = 1, 2, 3, 4)随入射波能流${P_{\rm{i}}}$的变化. 从图4(a)—(d)可见: 当其他参数不变时, 1)基频波和各高阶谐波的能流均随入射波能流的增加而增加, 基频波的能流明显大于各高阶谐波的能流; 2)基频波和各高阶谐波的能流均随黏度的增加而减小; 3)黏度对各高阶谐波能流的影响程度大于对基频波能流的影响程度.

图5为声波在黏性液体($\mu \ne 0$)中传播时, 黏度不同的情况下各阶谐波的相对能流${P_n}/{P_{\rm{i}}}$ (n = 1, 2)随角频率的变化. 从图5(a)和图5(b)可见: 1)基频波的相对能流随角频率的增加而减少, 二次谐波的相对能流随角频率的增加而增加; 2)随黏度的增加基频波和二次谐波的相对能流均减小, 但基频波的影响程度比二次谐波的影响程度更大一些.

4.   结　论

用一种新的数值计算方法研究了声波在液体中的传播特性. 得到的主要结论如下:

1)第1次迭代计算的输出结果$A_N^{(1)}$与线性理论的预测结果相同, 第2次迭代计算的输出结果$A_N^{(2)}\left( {N \ne 1} \right)$与微扰理论的预测结果相同. 在第2次计算中已经考虑了谐波的损耗问题, 然而在微扰理论中没有考虑到这一点;

2)声波在介质中传播时, 基频波的相对能流随入射波能流(或角频率)的增加而减少, 但各高阶谐波的相对能流却随入射波能流(或角频率)的增加而增加;

3)声波的各阶谐波能流(或相对能流)均随黏度的增加而减小;

4)在相同的参数下, 本文提出的研究计算方法比微扰法能更好地保证在研究声传播的过程中能量守恒, 且当涉及到更多的高次谐波时, 用本文提出的数值方法计算时间不会急剧增加.