富志凯， 石立华， 黄正宇， 付尚琛

(解放军理工大学 电磁环境效应及光电工程国家级重点实验室，南京 210007)

二维声波方程的Crank-Nicolson无条件稳定方法研究

富志凯， 石立华， 黄正宇， 付尚琛

(解放军理工大学 电磁环境效应及光电工程国家级重点实验室，南京 210007)

一阶速度-压力声波方程的有限差分数值模拟中，由于受到Courant-Friedrich-Levy (CFL)稳定性条件的限制，在分析精细结构问题时往往效率低下。将Crank-Nicolson(CN)方法引入到声波方程的有限差分模拟中，给出了声波方程的CN差分格式。通过Von Neuman 方法推导分析了CN方法的稳定性条件，证明了该方法的无条件稳定性。同时，采用非均匀网格技术进行网格剖分，进一步提高了仿真效率，减少了内存消耗。仿真实验中，建立了二维多层精细结构的声传播模型，通过与传统时域有限差分的仿真结果进行对比分析，验证了该方法的有效性。

声波方程；Crank-Nicolson方法；无条件稳定；非均匀网格

声波方程的时域有限差分模拟因其实现简单，计算效率高，被广泛应用于地震勘探[1-2]及复杂环境声传播[3]等领域。但是由于受到CFL稳定性条件限制，在处理某些精细结构时，空间网格必须划分得非常小，从而使得时间步长也必须非常小，大大降低了计算效率，有时这种限制对于整个时间区间的仿真是不可承受的。实际情况中，如对极薄地质层，极薄隔音墙[4]及声音消除器[5]等具有精细结构特征的物体进行有限差分模拟时，由于CFL稳定性条件限制，仿真将非常耗时，时间上将出现严重的过采样。因此，对声波方程开展无条件稳定研究，消除空间步长与时间步长之间的限制，提高计算效率，具有十分重要的意义。

声波方程的有限差分模拟中，国内外学者对二阶声波方程[6-9]及一阶速度-压力声波方程[10-11]进行大量仿真研究。一阶速度-压力声波方程虽然在空间网格剖分方面较二阶声波方程复杂，但因其能实时观测速度及压力两个场量，对实际物理问题的理解有着特殊的优势。在提高精度方面，Dablain等[12]将高阶有限差分方法运用到声波方程的正演模拟中，在相同空间网格长度的情况下大大提高了数值模拟精度；Cohen等[13]在Dablain的基础上，将高阶有限差分方法运用到非均匀介质中并取得较好的结果；Bayliss等[14]将高阶方法运用于弹性波方程中，得到(τ2+h4)的精度，其中τ与h分别代表时间步长与空间步长；Liu等[6,7]考虑到时间域使用高阶差分方法会引起不稳定，对色散关系进行泰勒级数展开，在时空联合域推导了新的差分格式，使得数值模拟的精度得到进一步提高。在提升效率方面，Peaceman等[15-16]首先引入交替隐式差分(ADI)方法解决热扩散问题并减少了CPU时间；Li等[17]首先提出了高阶紧支交替隐式差分方法，并用于解决抛物线方程的数值模拟问题；Qin[18]将此方法扩展到波动方程的数值模拟中，提高了计算效率。

在电磁波计算领域，许多国内外学者对精细结构的仿真效率问题进行深入研究，并提出了一系列无条件稳定的有限差分方法[19-21]，大大提高了仿真效率。在此基础上，本文将电磁波计算领域的CN-FDTD无条件稳定方法扩展到一阶速度-压力声波方程中，消除了时间步长与空间步长之间的稳定性条件限制，提高了仿真效率。同时，为了减少内存及消除不同空间步长导致的额外反射，本文在仿真中采用了渐变型的非均匀网格剖分技术。仿真实验中，建立了具有精细结构的2维声传播模型，激励源作用于2维空间中一点，精细结构两侧各有一点作为观测点。通过对比分析观测点的波形，验证了本文所提方法的有效性。

1 声波方程的无条件稳定算法

1.1Crank-Nicolson无条件稳定算法

为了方便分析，我们考虑密度恒定的二维一阶速度-压力声波方程，其形式可以表示为

(1)

式中：v,u分别是y与x方向的速度；p是声压；c是声速；Ω=[l1×l2]是计算空间范围；T是时间长度。将式(1)中的一阶空间导及一阶时间导用泰勒级数展开，得到传统的差分迭代格式为

(2a)

(2b)

(2c)

(3)

(4)

(5a)

(5b)

(5c)

通过对式(5)中的三个式子的调整，可以简化得到一个隐式方程

(6)

式(6)中，用中心差分格式对其离散化，我们可以发现每个点与其相邻的点有关系。将计算区域内所有的点按一定顺序集合起来，可写成矩阵形式

[M][P]=[S]

(7)

1.2非均匀网格剖分技术

考虑到进行大范围仿真模拟时，如果空间剖分过密，会严重影响计算效率及内存。为此，我们采用非均匀网格技术对空间进行剖分，在精细结构区域剖分得较为密集，在非精细结构区域区域剖分得稀疏些。但是，这种剖分同样会产生一个问题：当精细区域网格与其他区域网格大小差距较大时，会在交界处产生明显的反射。所以，在精细结构区域与非精细结构区域区域的中间需建立一块渐变过渡区域(如图1中网格渐变区域所示)，以避免网格大小差别太大引起反射。图1为二维空间的非均匀网格剖分示意图。

图1 非均匀网格剖分示意图

过渡区域内的空间网格步长Δx过渡满足如下变化规律：

Δx过渡=Δx精细·αn，α>1

(8)

式中：α为尺度变化因子；n表示过渡区域中空间网格总数。对尺度变化因子α来说，其取值越接近1，则由网格突变产生的额外反射就会越小，但是过小的α会使得过渡区网格总数大幅增加，会严重影响计算效率与内存消耗。因此，α的选取与实际需求有关，一般我们认为当α取值为1.25时产生的偏差在1%以内。

1.3吸收边界条件

为了截断计算区域，我们使用一阶Mur吸收边界来模拟无穷大区域[22]。当计算点在边界时，我们可以得到其吸收边界条件为

(9)

(10)

计算此方程得到的结果将具有边界吸收效果。

2 稳定性分析

2.1传统FDTD的稳定性条件

(11a)

(11b)

(11c)

(12)

为使式(12)中两根模不大于1，按实系数多项式根的定理[23],可得到不等式

(13)

2.2CN-FDTD的稳定性条件

按照以上的分析步骤，将带有增长因子的各场量式子代入式(5),可得到

(14a)

(14b)

(14c)

(15)

按实系数多项式根的定理，式(15)中两根不大于1的不等式为

(16)

由于不等式(16)恒成立，所以方程两根不大于1，即一阶速度-压力声波方程的CN差分格式是无条件稳定的。

3 仿真结果

仿真实验中，我们建立了一个具有精细结构的二维声传播模型，并采用了非均匀网格对其剖分。模型中将一极薄均匀介质层(d=0.01 m,c=6.0 km/s)夹在其他两层介质层中间，具体示意图如图2所示。图中，整个计算区域大小为60 m×90.01 m，并且极薄介质层位于x方向72 m处。介质A,极薄介质层和介质B的声传播速度分别为3.0 km/s、6.0 km/s和4.7 km/s。其中，Source为二维空间中的点激励源，在点(28 m,30 m)处；p1与p2为二维空间中的观测点，分别在点(36 m,30 m)与(82.01 m,30 m)处。为充分验证本方法可行性及有效性，仿真实验中分别设计了雷克子波与高斯脉冲作为点激励源。

图2 具有精细结构的二维声传播模型

3.1雷克子波源

对于传统FDTD方法，由于受到CFL稳定性条件的限制，我们将时间步长Δt=1.33×10-7s。而对于我们所提出的无条件稳定方法，虽然不再受到稳定性条件限制，但是完整的波形仍需要一定数量的时间采样点来支撑，即须满足采样定理。这里我们采取的原则是在满足采样定理的基础上，仍保证一定的余量，我们将CN方法的时间步长设定为Δt=1.67×10-4s。仿真实验中，为了对比的公平性，两种方法都使用了非均匀网格技术进行网格剖分。图3(a)和(b)分别为二维空间中p1与p2两点记录到的波形，其中实线为传统FDTD方法得到的结果，虚线为CN-FDTD方法得到的结果。从图中可知，本文所提出方法的仿真结果与传统FDTD方法有很好的一致性。

表1显示了传统FDTD方法与本文所提的CN-FDTD方法在消耗CPU时间方面的对比。当最小空间步长为0.001 m时，本文采用的时间步长是传统FDTD方法的1 250倍左右，并且总消耗CPU时间是传统方法的3.11%左右。

3.2高斯脉冲源

当激励源波形为高斯脉冲源时，点激励源同样作用在声压场上，其波形函数为f(t)=exp[-8(0.6f0t-1)2]，其中f0=200 Hz，仿真时间长度T=0.06 s。依据时间步长的设定原则，我们将传统FDTD方法与CN-FDTD方法的时间步长分别设为Δt=1.33×10-7s与Δt=1.67×10-4s。同样，在仿真实验中两种方法都使用了非均匀网格技术进行网格剖分。图4(a)和(b)分别为二维空间中p1与p2两点记录到的波形，其中实线为传统FDTD方法得到的结果，虚线为CN-FDTD方法得到的结果。从图中可知，本文所提出方法与传统FDTD方法得到的仿真波形吻合得非常好。

(b) p2点记录波形

方法Δx=0.001mΔtCPUtime传统FDTD方法1.33×10-7s787.8181sCN-FDTD方法1.67×10-4s24.4756s

表2显示了传统FDTD方法与本文所提的CN方法在消耗CPU时间方面的对比。当最小空间步长为0.001 m时，本文采用的时间步长是传统FDTD方法的1 250倍左右，并且总消耗CPU时间是传统方法的2.24%左右。

4 结 论

针对精细结构对声传播模拟耗时严重的问题，本文引入Crank-Nicolson方法，推导了二维声波方程的CN-FDTD格式，并采用Von Neumann法证明了此方法的无条件稳定性。仿真实验表明：本文所提方法得到的声传播波形与传统FDTD方法得到的结果吻合得非常好，当空间中存在精细结构时，本方法CPU耗时较传统FDTD方法减少了97%左右，极大提高了声波方程的仿真模拟效率。

(a) p1点记录波形

(b) p2点记录波形

方法Δx=0.001mΔtCPUtime传统FDTD方法1.33×10-7s836.8793sCN-FDTD方法1.67×10-4s18.7242s

[1] 梁文全，杨长春，王彦飞，等． 用于声波方程数值模拟的时间-空间域有限差分系数确定新方法[J]．地球物理学报, 2013, 56(10): 3497-3506.

LIANG Wenquan, YANG Changchun, WANG Yanfei, et al. Acoustic wave equation modeling with new time-space domain finite difference operators[J]. Chinese Journal of Geophysics, 2013, 56(10):3497-3506.

[2] 刘少林，杨长春，王彦飞，等. 求解声波方程的辛RKN格式[J]. 地球物理学报, 2013, 56(12): 4197-4205.

LIU Shaolin, YANG Changchun, WANG Yanfei, et al. Symplectic RKN schemes for seismic scalar wave simulations[J]. Chinese Jounral of Geophysics, 2013, 56(12) 4197-4205.

[3] LIU L B, ALBERT D G. Acoustic pulse propagation near a right-angle wall[J]. Journal of Acoustic Society of America, 2006, 119(4):2073-2083.

[4] CHO Y S. NDT response of spectral analysis of surface wave method to multi-layer thin high-strength concrete structures[J]. Ultrasonic, 2002, 40:227-230.

[5] YANG M, MENG C, FU C X, et al. Subwavelength total acoustic absorption with degenerate resonators[J]. Applied Physics Letters, 2015, 107(10): 1-7.

[6] LIU Y, SEN M. K A new time-space domain high-order finite-difference method for the acoustic wave equation[J]. Journal of Computational physics, 2009, 228: 8779-8806.

[7] LIU Y, SEN M K. Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation[J]. Journal of Computational Physics, 2013, 232: 327-345.

[8] 张红静，周辉，陈汉明 等.声波方程数值模拟中的任意广角单程波吸收边界[J]. 石油地球物理勘探, 2013,48(4):556-582.

ZHANG Hongjing, ZHOU Hui, CHEN Hanming, et al. Arbitrarily wide-angle wave equation absorbing boundary condition in acoustic numerical simulation[J]. Oil Geophysical Prospecting, 2013,48(4):556-582.

[9] LIAO W Y. On the dispersion, stability and accuracy of a compact higher-order finite difference scheme for 3D acoustic wave equation[J]. Journal of Computational and Applied Mathematics, 2014, 270:571-583.

[10] TAN S R, HUANG L. J A staggered-grid finite-difference scheme optimized in the time-space domain for modeling scalar-wave propagation in geophysical problems[J]. Journal of Computational Physics, 2014, 276: 613-634.

[11] CHEN H M, ZHOU H, ZHANG Q C, et al. A k-space operator-based least-squares staggered-grid finite-difference method for modeling scalar wave propagation[J]. Geophysics, 2016, 81(2):T39-T55.

[12] DABLAIN M A. The application in heterogeneous media: velocity-stress finite-difference method[J]. Geophysics, 1986, 51(1):54-66.

[13] COHEN G, JOLY P. Construction and analysis of fourth-order finite difference schemes for the acoustic wave equation in non-homogeneous media[J]. SIAM Journal on Numerical Analysis, 1996, 33(4):1266-1302.

[14] BAYLISS A, JORDAN K E, LEMESURIER B J, et al. A fourth-order accurate finite-difference scheme for the computation of elastic waves[J]. Bulletin of the Seismological Society of America, 1986, 76: 1115-1132.

[15] PEACEMAN D W, RACHFORD H. The numerical solution of parabolic and elliptic differential equations[J]. Journal of the Society for Industrial and Applied Mathematics, 1959, 3(3): 28-41

[16] DOUGLAS J, PEACEMAN D W. Numerical solution for two-dimensional heat flow problems[J]. American Institute of Chemical Engineers Journal, 1955, 1(4):505-512.

[17] LI J C, CHEN Y T, LIU G Q. High-order compact ADI methods for parabolic equations[J]. Computers and Mathematics with Applications, 2006, 52(8):1343-1356.

[18] QIN J G. The new alternating direction implicit difference methods for the wave equations[J]. Journal of Computational and Applied Mathematics, 2009, 230(1): 213-223.

[19] SUN G, TRUEMAN C W. Unconditionally stable Crank-Nicolson scheme for solving two-dimensional Maxwell’s equations[J]. Electronics Letters, 2003, 39(7):595-597.

[20] CHUNG Y S, SARKAR T K, JUNG B H, et al. An unconditionally stable scheme for the finite-difference time-domain method[J]. Transactions on Microwave Theory and Techniques, 2003, 51(3): 697-704.

[21] HUANG Z Y, SHI L H, CHEN B, et al. A new unconditionally stable scheme for FDTD method using associated hermite orthogonal functions[J]. IEEE Transactions on Antennas and Propagation, 2014, 62(9): 4804-4809.

[22] BI Z, WU K, WU C, et al. A dispersive boundary condition for micro strip component analysis using the FD-TD method[J]. Transactions on Microwave Theory and Techniques, 1992, 40(4): 774-777.

[23] 张文生. 科学计算中的偏微分方程有限差分法[J]. 北京：高等教育出版社, 2006.

ACrank-Nicolsonunconditionallystablemethodforsolving2Dacousticwaveequations

FU Zhikai, SHI Lihua, HUANG Zhengyu, FU Shangchen

(National Key Laboratory on Electromagnetic Environmental Effects and Electro-optical Engineering, PLA University of Science and Technology, Nanjing 210007, China)

Considering the constraint of Courant-Friedrich-Levy (CFL) stability condition, it is time-consuming to solve the one-order velocity-pressure acoustic wave equation with the conventional finite-difference time domain (FDTD) method, especially, to analyze fine structure problems. Here, Crank-Nicolson (CN) method was introduced in finite difference simulation, the acoustic wave equation’s CN difference scheme was obtained. Based on Von Neuman method, the unconditional stability condition of CN method for acoustic wave equations was derived. With the proposed method, the time step was not restricted by the CFL stability condition any more. Meanwhile, the non-uniform grid technology was used to generate mesh grids to further save internal memory and improve simulation efficiency. In simulation tests, a 2-dimentional multi-layer fine structure’s sound propagation model was established. Through comparing the simulation test results with those using the traditional FDTD method, the effectiveness of the proposed method was verified.

acoustic wave equation; Crank-Nicolson method; unconditionally stable; non-uniform grid

2016-03-25 修改稿收到日期：2016-07-14

富志凯 男，博士生，1987年2月生

石立华 男，教授，博士生导师.1969年7生生 E-mail:lihuashi@aliyun.com

P631.4

： A

10.13465/j.cnki.jvs.2017.17.013