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

全文HTML

--> --> -->

1.引　言

非线性声学是声学的一个重要分支, 当系统的声强比较强时就会产生各种非线性效应, 如谐波产生、冲击波形成、声辐射力的出现等. 非线性声学在某些领域有着重要的运用, 如高声强聚焦超声^[1-4]、超声悬浮^[5-7]、声空化^[8]、声谐波成像^[9,10]、参量发射阵^[11,12]等. 在这些领域中声波的非线性方程的求解是非常重要的.
关于非线性声学系统波动方程的求解, 现阶段用到的方法一般包括: 1)完全的数值计算方法, 如有限元和有限差分法^[13-17]. 这类方法得到解的物理意义并不明确, 很难揭示非线性作用的物理本质, 而且在很多情况下还会引起数值发散问题, 并非适用于全部的非线性问题. 2)严格的解析方法^[18-21]. 这种方法只能处理极少数系统的非线性声学问题, 如理想流体中非线性声波的传播. 3)微扰法^[22]. 它的优点是方法简单和解的物理意义清晰, 但是只适合处理低声强时的非线性效应. 且它只考虑低阶谐波对高阶谐波的作用, 而忽略其反作用, 因此并不满足能量守恒定律.
对于声学非线性方程的求解问题, 本文提出了一种新的、半解析的数值迭代方法. 它是在频域内把声场展开为傅里叶级数的形式, 实现时间变量和空间坐标的分离. 然后根据计算精度的具体需求, 截断高频谐波而实现方程的求解. 它的解具有非常清晰的物理意义, 即是各阶谐波的组合. 经过研究发现, 在微扰法适用的声强范围内, 本文提出的方法也是适用, 且满足能量守恒定律(无耗散的系统). 在微扰法不适用的一个较宽的声强范围内, 迭代法依然适用且满足能量守恒定律(无耗散的系统). 只是在极高声强的情况下, 本文提出的方法才不适用.

2.理论方法

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

-->

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}$

式中$\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}$

其中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}, $

$(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}, $

$(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}, $

$(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}, $

$(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}, $

$(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}, $

$(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}, $

其中

$\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} $

$\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}$

$\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} $

$\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}$

$\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} $

$\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} $

$\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}$

注意此处只给出了场变量的方程, 共轭场的方程并没有列出来, 只要对方程(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)}, $

$(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)}, $

$(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)}, $

$(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)}, $

$(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)}, $

$(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)}, $

$(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)}, $

分别用${A^{\left( {m - 1} \right)}}$