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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2a)

(2b)

(2c)

(3)

(4)

(5a)

(5b)

(5c)

通過(guò)對(duì)式(5)中的三個(gè)式子的調(diào)整，可以簡(jiǎn)化得到一個(gè)隱式方程

(6)

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

[M][P]=[S]

(7)

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

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

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

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

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

(8)

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

1.3吸收邊界條件

為了截?cái)嘤?jì)算區(qū)域，我們使用一階Mur吸收邊界來(lái)模擬無(wú)窮大區(qū)域[22]。當(dāng)計(jì)算點(diǎn)在邊界時(shí)，我們可以得到其吸收邊界條件為

(9)

(10)

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

2 穩(wěn)定性分析

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

(11a)

(11b)

(11c)

(12)

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

(13)

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

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

(14a)

(14b)

(14c)

(15)

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

(16)

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

3 仿真結(jié)果

仿真實(shí)驗(yàn)中，我們建立了一個(gè)具有精細(xì)結(jié)構(gòu)的二維聲傳播模型，并采用了非均勻網(wǎng)格對(duì)其剖分。模型中將一極薄均勻介質(zhì)層(d=0.01 m,c=6.0 km/s)夾在其他兩層介質(zhì)層中間，具體示意圖如圖2所示。圖中，整個(gè)計(jì)算區(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為二維空間中的點(diǎn)激勵(lì)源，在點(diǎn)(28 m,30 m)處；p1與p2為二維空間中的觀測(cè)點(diǎn)，分別在點(diǎn)(36 m,30 m)與(82.01 m,30 m)處。為充分驗(yàn)證本方法可行性及有效性，仿真實(shí)驗(yàn)中分別設(shè)計(jì)了雷克子波與高斯脈沖作為點(diǎn)激勵(lì)源。

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

3.1雷克子波源

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

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

3.2高斯脈沖源

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

(b) p2點(diǎn)記錄波形

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

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

4 結(jié) 論

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

(a) p1點(diǎn)記錄波形

(b) p2點(diǎn)記錄波形

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

[1] 梁文全，楊長(zhǎng)春，王彥飛，等． 用于聲波方程數(shù)值模擬的時(shí)間-空間域有限差分系數(shù)確定新方法[J]．地球物理學(xué)報(bào), 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] 劉少林，楊長(zhǎng)春，王彥飛，等. 求解聲波方程的辛RKN格式[J]. 地球物理學(xué)報(bào), 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ì)算中的偏微分方程有限差分法[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