二維聲波方程的Crank-Nicolson無條件穩定方法研究

富志凱， 石立華， 黃正宇， 付尚琛

(解放軍理工大學 電磁環境效應及光電工程國家級重點實驗室，南京 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