波动方程非线性传播仿真

孟艳

【摘要】波动方程是波动传播过程中，介质中各个质点振动物理特性的数学描述。非线性作用是介质的一种固有属性，与衰减作用类似，在波动传播的过程中改变波源的振动波形。通过对波动方程的非线性仿真，为进一步的研究波的非线性特性提供了基础

【关键词】波动方程 非线性 仿真

一、引言

非线性现象是指波在传播过程中，产生波源整数倍的高次频率，非线性是传播介质的基本特性。 本文以三维Westervelt波动方程为基础， 利用其时间域有限微分解， 通过计算机仿真的方式，完成对波动场内任意点的波動仿真， 为仿真研究非线性特性提供基础。

二、非线性波动方程的数值解

常见波动方程都是非常复杂的偏微分方程，其中声场描述参数声压p是一个与时空都有关系的物理量，在笛卡尔坐标系下可以表示为：P（x，y，z，t）。时域有限微分法的基本思想是使用微分替换方程中的导数，将未知的时间空间变量用已知时间或者空间变量表达。通过不断的重复，计算时间和空间上未来的结果。

为了研究方便，我们首先将Westervelt方程在一维坐标下展开化简为仅含对p的偏导数，其中拉普拉斯算子定义为空间的偏导数：

对于 的导数，代入Westervelt方程得到：

上式中包含压强P在时间和空间上的偏导数，在时间t区间内，振动的函数。每个坐标点 代表了t时刻，x轴上的坐标点i的压强。 的一阶导数表示为：

使用上面的微分形式替换导数形式，可以获得 的精度，其中h为时间或空间的步长。上面两个公式中时间和空间步长都为1。从上面的两个公式可以化解出：

从上面两个公式可以看出，空间上点i+1在时间t的值可以通过上一空间点（i-1）以及点i处的导数获得。在给定初始条件 条件下，使用前一个公式在x轴方向上计算下一个坐标位置的值，使用后一个公式计算t+1时刻的坐标值，这样在x-t平面延伸，最后获得x轴方向上，t时间区间内，所有的点的压强。

三、仿真结果

从公式可以看出，非线性系数和液体黏滞系数分别作用到不同的压强分量上，我们可以将β与η分别设置为0来研究衰减和非线性对波动传播的影响。

首先我们以人体组织介质的物理参数值作为模拟的基础，假设x方向的长度为200 ，在100 处放置一个振源，振源的激励波为：

我们首先假设介质黏滞系数η及非线性系数β都为0，使用FDTD计算波动在x轴方向传播。经过400 时刻后，x轴方向上质点的振幅如图4.1所示：

对比图4.1中波形可以看出，在线性介质中，波动波形能够保持振源振动波形的光滑特性，但是，在非线性介质中，质点压缩周期的斜率高于线性作用下的传播。下面，我们对比非线性介质与线性介质中，距离振源长度为10 和20 处，质点在[0，400 ]时间周期内，振动的波形曲线和频谱曲线。

线性介质中传播的质点先到达最大距离。由于距离的导数是速度，那么，从振动曲线的斜率可以看出质点在对应时间区间内振动的速度。也就是说，在非线性介质中，时间区域[100，150]，[250，350]内，质点振动速度大于线性介质。这与非线性产生的物理解释完全吻合。

【参考文献】

【1】J.F.Greenleaf Measurement of the acoustic nonlinearity parameter B/A in human tissue by a thermodynamic method. The Journal of the Acoustical Society of America， vol. 76， no4， pp. 1023-1029

【2】V. P. Kuznetsov， “Equations of Nonlinear Acoustics，” Sov. Phys. Acoust. 16， 467-470.

【3】Y.-S. Lee， "Numerical solution of the KZK equation for pulsed finite amplitude sound beams in thermo viscous fluids，" Ph.D. dissertation， The University of Texas at Austin

【4】凌涛 ，郑海荣，“超声造影微泡非线性声学特性与成像研究进展”中国介入影像与治疗学2009年第6卷第4期 p392～p396

【5】TL SZABO Diagnostic Ultrasound Imaging Inside Out，Appendix B