二維聲波方程的Crank-Nicolson無條件穩(wěn)定方法研究

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

(解放軍理工大學(xué) 電磁環(huán)境效應(yīng)及光電工程國家級重點實驗室，南京 210007)

二維聲波方程的Crank-Nicolson無條件穩(wěn)定方法研究

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

(解放軍理工大學(xué) 電磁環(huán)境效應(yīng)及光電工程國家級重點實驗室，南京 210007)

一階速度-壓力聲波方程的有限差分?jǐn)?shù)值模擬中，由于受到Courant-Friedrich-Levy (CFL)穩(wěn)定性條件的限制，在分析精細(xì)結(jié)構(gòu)問題時往往效率低下。將Crank-Nicolson(CN)方法引入到聲波方程的有限差分模擬中，給出了聲波方程的CN差分格式。通過Von Neuman 方法推導(dǎo)分析了CN方法的穩(wěn)定性條件，證明了該方法的無條件穩(wěn)定性。同時，采用非均勻網(wǎng)格技術(shù)進(jìn)行網(wǎng)格剖分，進(jìn)一步提高了仿真效率，減少了內(nèi)存消耗。仿真實驗中，建立了二維多層精細(xì)結(jié)構(gòu)的聲傳播模型，通過與傳統(tǒng)時域有限差分的仿真結(jié)果進(jìn)行對比分析，驗證了該方法的有效性。

聲波方程；Crank-Nicolson方法；無條件穩(wěn)定；非均勻網(wǎng)格

聲波方程的時域有限差分模擬因其實現(xiàn)簡單，計算效率高，被廣泛應(yīng)用于地震勘探[1-2]及復(fù)雜環(huán)境聲傳播[3]等領(lǐng)域。但是由于受到CFL穩(wěn)定性條件限制，在處理某些精細(xì)結(jié)構(gòu)時，空間網(wǎng)格必須劃分得非常小，從而使得時間步長也必須非常小，大大降低了計算效率，有時這種限制對于整個時間區(qū)間的仿真是不可承受的。實際情況中，如對極薄地質(zhì)層，極薄隔音墻[4]及聲音消除器[5]等具有精細(xì)結(jié)構(gòu)特征的物體進(jìn)行有限差分模擬時，由于CFL穩(wěn)定性條件限制，仿真將非常耗時，時間上將出現(xiàn)嚴(yán)重的過采樣。因此，對聲波方程開展無條件穩(wěn)定研究，消除空間步長與時間步長之間的限制，提高計算效率，具有十分重要的意義。

聲波方程的有限差分模擬中，國內(nèi)外學(xué)者對二階聲波方程[6-9]及一階速度-壓力聲波方程[10-11]進(jìn)行大量仿真研究。一階速度-壓力聲波方程雖然在空間網(wǎng)格剖分方面較二階聲波方程復(fù)雜，但因其能實時觀測速度及壓力兩個場量，對實際物理問題的理解有著特殊的優(yōu)勢。在提高精度方面，Dablain等[12]將高階有限差分方法運(yùn)用到聲波方程的正演模擬中，在相同空間網(wǎng)格長度的情況下大大提高了數(shù)值模擬精度；Cohen等[13]在Dablain的基礎(chǔ)上，將高階有限差分方法運(yùn)用到非均勻介質(zhì)中并取得較好的結(jié)果；Bayliss等[14]將高階方法運(yùn)用于彈性波方程中，得到(τ2+h4)的精度，其中τ與h分別代表時間步長與空間步長；Liu等[6,7]考慮到時間域使用高階差分方法會引起不穩(wěn)定，對色散關(guān)系進(jìn)行泰勒級數(shù)展開，在時空聯(lián)合域推導(dǎo)了新的差分格式，使得數(shù)值模擬的精度得到進(jìn)一步提高。在提升效率方面，Peaceman等[15-16]首先引入交替隱式差分(ADI)方法解決熱擴(kuò)散問題并減少了CPU時間；Li等[17]首先提出了高階緊支交替隱式差分方法，并用于解決拋物線方程的數(shù)值模擬問題；Qin[18]將此方法擴(kuò)展到波動方程的數(shù)值模擬中，提高了計算效率。

在電磁波計算領(lǐng)域，許多國內(nèi)外學(xué)者對精細(xì)結(jié)構(gòu)的仿真效率問題進(jìn)行深入研究，并提出了一系列無條件穩(wěn)定的有限差分方法[19-21]，大大提高了仿真效率。在此基礎(chǔ)上，本文將電磁波計算領(lǐng)域的CN-FDTD無條件穩(wěn)定方法擴(kuò)展到一階速度-壓力聲波方程中，消除了時間步長與空間步長之間的穩(wěn)定性條件限制，提高了仿真效率。同時，為了減少內(nèi)存及消除不同空間步長導(dǎo)致的額外反射，本文在仿真中采用了漸變型的非均勻網(wǎng)格剖分技術(shù)。仿真實驗中，建立了具有精細(xì)結(jié)構(gòu)的2維聲傳播模型，激勵源作用于2維空間中一點，精細(xì)結(jié)構(gòu)兩側(cè)各有一點作為觀測點。通過對比分析觀測點的波形，驗證了本文所提方法的有效性。

1 聲波方程的無條件穩(wěn)定算法

1.1Crank-Nicolson無條件穩(wěn)定算法

為了方便分析，我們考慮密度恒定的二維一階速度-壓力聲波方程，其形式可以表示為

(1)

式中：v,u分別是y與x方向的速度；p是聲壓；c是聲速；Ω=[l1×l2]是計算空間范圍；T是時間長度。將式(1)中的一階空間導(dǎo)及一階時間導(dǎo)用泰勒級數(shù)展開，得到傳統(tǒng)的差分迭代格式為

(2a)

(2b)

(2c)

(3)

(4)

(5a)

(5b)

(5c)

通過對式(5)中的三個式子的調(diào)整，可以簡化得到一個隱式方程

(6)

式(6)中，用中心差分格式對其離散化，我們可以發(fā)現(xiàn)每個點與其相鄰的點有關(guān)系。將計算區(qū)域內(nèi)所有的點按一定順序集合起來，可寫成矩陣形式

[M][P]=[S]

(7)

1.2非均勻網(wǎng)格剖分技術(shù)

考慮到進(jìn)行大范圍仿真模擬時，如果空間剖分過密，會嚴(yán)重影響計算效率及內(nèi)存。為此，我們采用非均勻網(wǎng)格技術(shù)對空間進(jìn)行剖分，在精細(xì)結(jié)構(gòu)區(qū)域剖分得較為密集，在非精細(xì)結(jié)構(gòu)區(qū)域區(qū)域剖分得稀疏些。但是，這種剖分同樣會產(chǎn)生一個問題：當(dāng)精細(xì)區(qū)域網(wǎng)格與其他區(qū)域網(wǎng)格大小差距較大時，會在交界處產(chǎn)生明顯的反射。所以，在精細(xì)結(jié)構(gòu)區(qū)域與非精細(xì)結(jié)構(gòu)區(qū)域區(qū)域的中間需建立一塊漸變過渡區(qū)域(如圖1中網(wǎng)格漸變區(qū)域所示)，以避免網(wǎng)格大小差別太大引起反射。圖1為二維空間的非均勻網(wǎng)格剖分示意圖。

圖1 非均勻網(wǎng)格剖分示意圖

過渡區(qū)域內(nèi)的空間網(wǎng)格步長Δx過渡滿足如下變化規(guī)律：

Δx過渡=Δx精細(xì)·αn，α>1

(8)

式中：α為尺度變化因子；n表示過渡區(qū)域中空間網(wǎng)格總數(shù)。對尺度變化因子α來說，其取值越接近1，則由網(wǎng)格突變產(chǎn)生的額外反射就會越小，但是過小的α?xí)沟眠^渡區(qū)網(wǎng)格總數(shù)大幅增加，會嚴(yán)重影響計算效率與內(nèi)存消耗。因此，α的選取與實際需求有關(guān)，一般我們認(rèn)為當(dāng)α取值為1.25時產(chǎn)生的偏差在1%以內(nèi)。

1.3吸收邊界條件

為了截斷計算區(qū)域，我們使用一階Mur吸收邊界來模擬無窮大區(qū)域[22]。當(dāng)計算點在邊界時，我們可以得到其吸收邊界條件為

(9)

(10)

計算此方程得到的結(jié)果將具有邊界吸收效果。

2 穩(wěn)定性分析

2.1傳統(tǒng)FDTD的穩(wěn)定性條件

(11a)

(11b)

(11c)

(12)

為使式(12)中兩根模不大于1，按實系數(shù)多項式根的定理[23],可得到不等式

(13)

2.2CN-FDTD的穩(wěn)定性條件

按照以上的分析步驟，將帶有增長因子的各場量式子代入式(5),可得到

(14a)

(14b)

(14c)

(15)

按實系數(shù)多項式根的定理，式(15)中兩根不大于1的不等式為

(16)

由于不等式(16)恒成立，所以方程兩根不大于1，即一階速度-壓力聲波方程的CN差分格式是無條件穩(wěn)定的。

3 仿真結(jié)果

仿真實驗中，我們建立了一個具有精細(xì)結(jié)構(gòu)的二維聲傳播模型，并采用了非均勻網(wǎng)格對其剖分。模型中將一極薄均勻介質(zhì)層(d=0.01 m,c=6.0 km/s)夾在其他兩層介質(zhì)層中間，具體示意圖如圖2所示。圖中，整個計算區(qū)域大小為60 m×90.01 m，并且極薄介質(zhì)層位于x方向72 m處。介質(zhì)A,極薄介質(zhì)層和介質(zhì)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)處。為充分驗證本方法可行性及有效性，仿真實驗中分別設(shè)計了雷克子波與高斯脈沖作為點激勵源。

圖2 具有精細(xì)結(jié)構(gòu)的二維聲傳播模型

3.1雷克子波源

對于傳統(tǒng)FDTD方法，由于受到CFL穩(wěn)定性條件的限制，我們將時間步長Δt=1.33×10-7s。而對于我們所提出的無條件穩(wěn)定方法，雖然不再受到穩(wěn)定性條件限制，但是完整的波形仍需要一定數(shù)量的時間采樣點來支撐，即須滿足采樣定理。這里我們采取的原則是在滿足采樣定理的基礎(chǔ)上，仍保證一定的余量，我們將CN方法的時間步長設(shè)定為Δt=1.67×10-4s。仿真實驗中，為了對比的公平性，兩種方法都使用了非均勻網(wǎng)格技術(shù)進(jìn)行網(wǎng)格剖分。圖3(a)和(b)分別為二維空間中p1與p2兩點記錄到的波形，其中實線為傳統(tǒng)FDTD方法得到的結(jié)果，虛線為CN-FDTD方法得到的結(jié)果。從圖中可知，本文所提出方法的仿真結(jié)果與傳統(tǒng)FDTD方法有很好的一致性。

表1顯示了傳統(tǒng)FDTD方法與本文所提的CN-FDTD方法在消耗CPU時間方面的對比。當(dāng)最小空間步長為0.001 m時，本文采用的時間步長是傳統(tǒng)FDTD方法的1 250倍左右，并且總消耗CPU時間是傳統(tǒng)方法的3.11%左右。

3.2高斯脈沖源

當(dāng)激勵源波形為高斯脈沖源時，點激勵源同樣作用在聲壓場上，其波形函數(shù)為f(t)=exp[-8(0.6f0t-1)2]，其中f0=200 Hz，仿真時間長度T=0.06 s。依據(jù)時間步長的設(shè)定原則，我們將傳統(tǒng)FDTD方法與CN-FDTD方法的時間步長分別設(shè)為Δt=1.33×10-7s與Δt=1.67×10-4s。同樣，在仿真實驗中兩種方法都使用了非均勻網(wǎng)格技術(shù)進(jìn)行網(wǎng)格剖分。圖4(a)和(b)分別為二維空間中p1與p2兩點記錄到的波形，其中實線為傳統(tǒng)FDTD方法得到的結(jié)果，虛線為CN-FDTD方法得到的結(jié)果。從圖中可知，本文所提出方法與傳統(tǒng)FDTD方法得到的仿真波形吻合得非常好。

(b) p2點記錄波形

方法Δx=0.001mΔtCPUtime傳統(tǒng)FDTD方法1.33×10-7s787.8181sCN-FDTD方法1.67×10-4s24.4756s

表2顯示了傳統(tǒng)FDTD方法與本文所提的CN方法在消耗CPU時間方面的對比。當(dāng)最小空間步長為0.001 m時，本文采用的時間步長是傳統(tǒng)FDTD方法的1 250倍左右，并且總消耗CPU時間是傳統(tǒng)方法的2.24%左右。

4 結(jié) 論

針對精細(xì)結(jié)構(gòu)對聲傳播模擬耗時嚴(yán)重的問題，本文引入Crank-Nicolson方法，推導(dǎo)了二維聲波方程的CN-FDTD格式，并采用Von Neumann法證明了此方法的無條件穩(wěn)定性。仿真實驗表明：本文所提方法得到的聲傳播波形與傳統(tǒng)FDTD方法得到的結(jié)果吻合得非常好，當(dāng)空間中存在精細(xì)結(jié)構(gòu)時，本方法CPU耗時較傳統(tǒng)FDTD方法減少了97%左右，極大提高了聲波方程的仿真模擬效率。

(a) p1點記錄波形

(b) p2點記錄波形

方法Δx=0.001mΔtCPUtime傳統(tǒng)FDTD方法1.33×10-7s836.8793sCN-FDTD方法1.67×10-4s18.7242s

[1] 梁文全，楊長春，王彥飛，等． 用于聲波方程數(shù)值模擬的時間-空間域有限差分系數(shù)確定新方法[J]．地球物理學(xué)報, 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]. 地球物理學(xué)報, 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] 張紅靜，周輝，陳漢明 等.聲波方程數(shù)值模擬中的任意廣角單程波吸收邊界[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] 張文生. 科學(xué)計算中的偏微分方程有限差分法[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月生

石立華 男，教授，博士生導(dǎo)師.1969年7生生 E-mail:lihuashi@aliyun.com

P631.4

： A

10.13465/j.cnki.jvs.2017.17.013