徐 娜 李 洋 周正干

(北京航空航天大學(xué) 機械工程及自動化學(xué)院，北京100191)

于 光

(上海飛機制造有限公司 航空制造技術(shù)研究所，上海200436)

超聲檢測具有檢測范圍廣、檢測深度大、使用方便、速度快等特點，是國內(nèi)外應(yīng)用最廣泛的無損檢測技術(shù)之一.然而，由于聲學(xué)的實際應(yīng)用問題很少能有嚴(yán)密的理論解，因此，從超聲傳播理論入手，研究聲波在介質(zhì)中的傳播及與缺陷的相互作用規(guī)律，建立超聲檢測系統(tǒng)模型，是超聲檢測的重要發(fā)展方向之一［1］.

時域有限差分法 (FDTD，F(xiàn)inite-Difference Time Domain)最初是一種電磁場數(shù)值計算的方法，目前被廣泛用于模擬彈性波在復(fù)雜介質(zhì)中的傳播.其差分格式中包含有介質(zhì)的參量，只需賦予各網(wǎng)格相應(yīng)的參量，即可模擬各種復(fù)雜的結(jié)構(gòu).另外，該方法采用步進法進行計算，易于實現(xiàn)各種復(fù)雜時域?qū)拵盘柕哪M，可方便獲得超聲波聲場在空間分布的可視化及某一空間點的時域信號波形［2］.采用FDTD方法對超聲波聲場進行數(shù)值計算，目前已有許多報道，如文獻［2］應(yīng)用FDTD方法對彈性固體中的聲場方程進行了二維聲場特性研究，文獻［3］采用FDTD方法對相控陣超聲檢測中聲場與缺陷的相互作用進行了數(shù)值仿真.

通常情況下，為了保證算法穩(wěn)定性，F(xiàn)DTD方法計算參數(shù)的選取需要滿足兩個限制條件，即仿真空間間隔遠(yuǎn)小于超聲波在被檢試件中傳播的波長，仿真時間間隔應(yīng)根據(jù)Courant限制條件來選取［4-5］.但是，在實際聲場計算中發(fā)現(xiàn)，當(dāng)相鄰兩種介質(zhì)的聲學(xué)參數(shù)相差較大時，采用Courant限制條件來選取計算參數(shù)無法保證算法的穩(wěn)定性［5］，大大降低了FDTD方法在超聲波聲場計算中的實用性.并且，在實際工程檢測中，常會遇到被檢對象中相鄰兩種介質(zhì)的聲阻抗相差很大的情況，或者耦合介質(zhì)與被檢工件的聲阻抗相差很大的情況(例如采用空氣作為耦合介質(zhì)).因此，如何有效解決FDTD方法數(shù)值穩(wěn)定性的問題，顯得尤為重要.

為了解決FDTD方法數(shù)值穩(wěn)定性的問題，很多學(xué)者都做過探討研究，提出了一些改進方法［6-8］，但大多應(yīng)用在電磁場數(shù)值計算中，很少研究其在超聲波聲場計算中的應(yīng)用.因此，本文研究FDTD方法在超聲波聲場計算中的應(yīng)用，通過對FDTD方法的改進擬解決其數(shù)值穩(wěn)定性的問題，從而改善其在超聲波聲場計算中的實用性.

1 FDTD方法及其改進

1.1 常規(guī)二維彈性波動方程的差分計算

聲場變量常用質(zhì)點速度ν和應(yīng)力σ表示，二維x-z平面彈性波動方程［9］為

式中，ρ為質(zhì)點密度;νx和νz分別為x和z方向的質(zhì)點速度;σx和σz分別為空間各點x和z方向的正應(yīng)力;σxz為切應(yīng)力.式(2)中的系數(shù)矩陣由材料的彈性模量及泊松比確定，對于各向同性介質(zhì)，該系數(shù)矩陣滿足:

式中，λ和μ是拉梅常數(shù)，也稱為切變彈性系數(shù).

取Δx和Δz分別為x和z方向的空間采樣間隔，T為時間采樣間隔，使用常規(guī)的有限差分法推出的差分計算公式，如式(4)～式(8):

實際計算過程中，按照如下步驟進行計算:

1)利用式(4)和式(5)計算(n+1/2)T時刻空間處質(zhì)點速度νx和νz的值;

2)利用式(6)～式(8)計算(n+1)T時刻空間各處應(yīng)力σx，σz及切應(yīng)力σxz的值;

3)回到步驟1)，將計算得出的(n+1/2)T及(n+1)T時刻的值分別作為(n-1/2)T及nT時刻的值，計算得出下一時刻的 νx，νz，σx，σz及σxz的值，不斷重復(fù)迭代從而實現(xiàn)數(shù)值模擬.

1.2 FDTD方法的改進

由于差分是利用空間某點及其周圍點對該點的描述，因此對微分方程進行差分離散后不可避免地會引入其他空間點的信息，而當(dāng)引入的其他空間點的屬性與原來該點的屬性相差較大時，差分方程將無法正確表述原微分方程.因此，當(dāng)相鄰兩種介質(zhì)的聲學(xué)參數(shù)相差較大時，在兩種介質(zhì)的分界面上經(jīng)過差分計算的方程組就無法正確表達(dá)原微分方程，從而引起差分算法的不穩(wěn)定.

針對該問題，提出將界面兩邊介質(zhì)的聲學(xué)參數(shù)在界面處進行平均處理，從而緩解兩者介質(zhì)聲學(xué)參數(shù)在界面處的突變問題，減小數(shù)值計算的誤差，提高算法的穩(wěn)定性.對FDTD方法進行改進，提出對波動方程式(1)和式(2)先進行積分處理，然后再進行差分離散處理，使由于引入其他空間點的信息而產(chǎn)生的誤差在很大程度上得到削弱.

以圖1為例介紹改進的FDTD方法.圖1是波動方程差分計算的離散網(wǎng)格示意圖，其中質(zhì)點速度νx，νz以及切應(yīng)力σxz處于介質(zhì)1和介質(zhì)2的界面上.首先對式(1)的方程1進行積分處理:

其積分范圍為

圖1 差分算法的離散網(wǎng)格示意圖

根據(jù)積分中值定理，對式(9)進行離散近似，可得

采用相同的積分計算方法對式(1)的方程2以及式(2)的3個方程進行計算，可得式(11)～式(14):

式(11)的積分范圍為

式(12)、式(13)積分范圍為

式(14)的積分范圍為

通過計算，使用改進的有限差分法推出的差分計算表達(dá)式如式(15)～式(19)所示.可看出式(17)和式(18)與利用常規(guī)FDTD方法獲得的差分表達(dá)式(式(6)、式(7))是一致的，這是因為正應(yīng)力σx和σz沒有位于兩種介質(zhì)的分界面上.但σxz位于其上，所以如式(19)所示，對式(2)中方程3的系數(shù)項C55進行了簡單平均處理.

通過該方法得到的差分表達(dá)式實質(zhì)上是對常規(guī)差分格式中相鄰兩種介質(zhì)分界面上的系數(shù)項進行了簡單平均處理，從而解決相鄰兩種介質(zhì)聲學(xué)參數(shù)相差太大而造成的差分計算不穩(wěn)定的問題.

2 FDTD方法改進前后的對比

建立圖2所示的仿真模型.試樣參數(shù):空氣層，厚度7.5 mm，密度 ρ=1.169 1 kg/m3，縱波聲速 νl=340m/s，聲阻抗 Z0=3.97 ×102Pa·s/m3;金屬層，厚度12.5 mm，密度ρ=7700 kg/m3，橫波聲速 νt=3230 m/s，縱波聲速 νl=5950 m/s，聲阻抗Z1=4.58×107Pa·s/m3.相鄰兩種介質(zhì)的聲阻抗相差5個數(shù)量級.

圖2 仿真模型示意圖

由于計算機容量限制，F(xiàn)DTD技術(shù)只能在有限區(qū)域進行，因此邊界將采用完全匹配層(PML，Perfectly Matched Layer)方法進行處理［2］.

按照仿真空間間隔遠(yuǎn)小于超聲波在被檢試件中傳播波長的要求以及仿真時間間隔根據(jù)Courant限制條件選取的方法，取空間間隔h=50μm，時間間隔T=5 ns.當(dāng)采用常規(guī)FDTD方法進行聲場模擬時，仿真程序運行一段時間后，便出現(xiàn)圖3所示現(xiàn)象，仿真程序無法繼續(xù)進行.

圖3 利用常規(guī)有限差分法計算的聲場圖

在相同空間和時間間隔下，采用改進的FDTD方法進行聲場模擬時，程序運行正常，如圖4所示.可見，改進的FDTD方法保證了差分計算的穩(wěn)定性.

圖4 利用改進有限差分法計算的聲場圖

在犧牲仿真時間的情況下，將計算參數(shù)的空間間隔和時間間隔均縮小5倍，即空間間隔h=10 μm，時間間隔T=1 ns時，采用常規(guī)FDTD方法再次進行聲場模擬，仍出現(xiàn)圖3所示發(fā)散現(xiàn)象.采用改進的FDTD方法進行聲場模擬，程序正常運行.

上述數(shù)值仿真實驗表明，當(dāng)相鄰兩種介質(zhì)的聲阻抗相差較大時，即使采用Courant條件選取適當(dāng)?shù)目臻g和時間間隔，甚至將空間和時間間隔縮小5倍的情況下，仍會引起差分計算的不穩(wěn)定.而采用改進的FDTD方法時可以有效解決發(fā)散現(xiàn)象，保證了差分計算的穩(wěn)定性.

3 FDTD改進方法的正確性驗證

3.1 仿真實驗

設(shè)計如圖5所示的檢測試樣，試塊中縱波聲速5 900 m/s，橫波聲速 3 230 m/s，密度 7 900 kg/m3.

按照圖5所示建立仿真模型，仿真參數(shù):空間間隔h=0.1 mm，時間間隔T=10 ns.激勵信號為1個周期2.5 MHz的正弦信號.

圖5 檢測試樣

在試樣上表面中心位置處(圖5中A點)激勵信號，圖6所示為某時刻的聲場空間分布情況.取應(yīng)力來表征回波信號幅值，并做歸一化處理，所獲得的信號A型圖如圖7所示.

圖6 超聲聲場分布圖

圖7 仿真信號A型圖

3.2 檢測實驗

采用單晶探頭檢測試樣，探頭中心頻率2.5 MHz，直徑20 mm，采用接觸法檢測.當(dāng)探頭放于試樣中心位置時，超聲信號A型圖見圖8.

對比圖7和圖8，仿真和實驗的結(jié)果均可識別出試樣上表面的界面波、缺陷波和圓弧面反射波，且時間位置與理論值一致.

仿真與檢測結(jié)果表明，改進的有限差分算法可以用于超聲檢測的數(shù)值仿真.

圖8 實驗信號A型圖

4 結(jié)論

1)FDTD的改進方法有效解決了相鄰兩種介質(zhì)聲阻抗相差較大時的算法穩(wěn)定性問題;

2)FDTD的改進方法能夠用于彈性固體中超聲波傳播的仿真計算，是超聲波聲場特性分析的一種有效數(shù)值方法.

References)

［1］丁輝.計算超聲學(xué)-聲場分析及應(yīng)用［M］.北京:科學(xué)出版社，2010:1-9 Ding Hui.Computational ultrasonics-analysis application of ultrasonic field［M］.Beijing:Science Press，2010:1-9(in Chinese)

［2］魏東，周正干.固體中脈沖超聲波傳播的有限差分模擬［J］.航空學(xué)報，2010，31(2):388-392 Wei Dong，Zhou Zhenggan.Finite difference simulation of pulsed ultrasonic propagation in solids［J］.Acta Aeronautica et Astronautica Sinica，2010，31(2):388-392(in Chinese)

［3］ Satyanarayan L，Sridhar C，Krishnamurthy C V，et al.Simulation of ultrasonic phased array technique for imaging and sizing of defects using longitudinal waves［J］.International Journal of Pressure Vessels and Piping，2007，84:716-729

［4］ Zhen Fenghua，Chen Zhizhang，Zhang Jiazong.Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method［J］.Microwave Theory and Techniques，2000，48(9):1550-1558

［5］ Schroder C T.On the stability of the FDTD algorithm for elastic media at a material interface［J］.Transactions on Geoscience and Remote Sensing，2002，40(2):474-481

［6］ Thoma P，Weiland T.Numerical stability of finite difference time domain methods［J］.Magnetics，1998，34(5):2740-2743

［7］ Remis R F.On the stability of the finite-difference time-domain method［J］.Journal of Computational Physics，2000，163(1):249-261

［8］ Pereda Jos'e A，Oscar Garc'a，Angel Vegas，et al.Numerical dispersion and stability analysis of the FDTD technique in lossy dielectrics［J］.Microwave and Guided Wave Letters，1998，8(7):245-247

［9］李太寶.計算聲學(xué)-聲場的方程和計算方法［M］.北京:科學(xué)出版社，2003:67-83 Li Taibao.Computational acoustics-sound field equation and calculate method［M］.Beijing:Science Press，2003:67-83(in Chinese)