非線性聲波方程的幾種解法對(duì)比

張世功，丁 凱，張克聲，蘇向東

(1. 貴州理工學(xué)院，貴州 貴陽 550003；2. 貴州省醫(yī)工交叉工程研究中心，貴州 貴陽 550003；3. 近地面探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 無錫 214035)

0 引 言

測(cè)量介質(zhì)的非線性系數(shù)在醫(yī)學(xué)診療和材料的無損檢測(cè)上都有至關(guān)重要的作用。生物組織的聲阻抗在發(fā)生病變時(shí)會(huì)產(chǎn)生微小的力學(xué)性能變化，通過表征材料非線性系數(shù)可以找到不易(被線性超聲)發(fā)現(xiàn)的病變組織(診斷成像)并進(jìn)行后續(xù)治療[1-2]。而在固體介質(zhì)中，利用傳統(tǒng)的無損檢測(cè)手段很難發(fā)現(xiàn)微裂紋、熱損傷、應(yīng)力集中等早期危害。而表征非線性系數(shù)可及時(shí)發(fā)現(xiàn)這些早期損傷，對(duì)其干預(yù)可以減小服役過程可能會(huì)帶來的損失[3-5]。

眾所周知，大多非線性問題理論上難以解決，常用的攝動(dòng)理論只是一種近似方法。利用它圍繞非線性系數(shù)對(duì)非線性聲波方程進(jìn)行攝動(dòng)展開，能得到一系列解。原則上，所有階次解之和才是非線性聲波方程的解析解，但要得到高階次的攝動(dòng)解非常困難，而低階次的近似攝動(dòng)解的誤差較大，適用范圍也相對(duì)有限。數(shù)值仿真計(jì)算，比如有限元、有限差分等方法[6-7]也可以用來計(jì)算非線性方程。然而，數(shù)值解得到的都是數(shù)值，不易用來分析相關(guān)的物理規(guī)律。另外，非線性聲學(xué)中會(huì)引起高次諧波的產(chǎn)生，有限的空間步長(zhǎng)和時(shí)間步長(zhǎng)能造成計(jì)算時(shí)的發(fā)散問題。

目前，實(shí)驗(yàn)表征材料非線性系數(shù)時(shí)常用二階攝動(dòng)解進(jìn)行計(jì)算[5,8-9]，但二階攝動(dòng)近似解只能在近場(chǎng)或小振幅激勵(lì)信號(hào)才適用，當(dāng)樣品不太小或激勵(lì)源不太弱時(shí)，得到的實(shí)驗(yàn)結(jié)果會(huì)有較大誤差[9]。

在介紹氣、液、固三種狀態(tài)介質(zhì)中非線性聲波方程的基礎(chǔ)上，本文闡述了多種求解一維非線性聲波方程的方法。首先，通過理論分析，三種狀態(tài)中的非線性聲波方程具有相似的形式，這表明不同狀態(tài)下的非線性聲波應(yīng)有相似的傳播性質(zhì)。其次，利用符號(hào)計(jì)算工具可獲得8階攝動(dòng)解，并與有限元、有限差分兩種數(shù)值解比較，將得到的非線性聲波時(shí)域波形與實(shí)驗(yàn)進(jìn)行了對(duì)比。再次，通過比對(duì)線性方程和非線性方程的形式，采用線性方程的求解方法對(duì)非線性聲波方程求解，得到非線性方程的偽線性解。上述多種解得到的時(shí)域波形都顯示了波形畸變等非線性聲波的傳播性質(zhì)。為證實(shí)這些解的非線性傳播性質(zhì)，在水中開展非線性聲場(chǎng)的測(cè)量實(shí)驗(yàn)，驗(yàn)證了這些解的有效性。高階攝動(dòng)解能用來更精確地測(cè)量介質(zhì)的非線性系數(shù)。當(dāng)然，非線性聲波方程的多種解法也可為非線性聲學(xué)的其他相關(guān)深入研究提供理論依據(jù)。

1 氣液固介質(zhì)中的非線性聲波方程

固體中的非線性包括材料非線性和幾何非線性[10-11]。材料非線性源自原子的不規(guī)則排列，主要表現(xiàn)為應(yīng)力和應(yīng)變的非線性關(guān)系。而幾何非線性（也稱非經(jīng)典非線性[12]）主要考慮的是大的幾何變形引起的非線性問題。本文主要討論的是材料非線性問題。

利用應(yīng)變能公式[13-14]可以得到各向同性固體介質(zhì)中的一維非線性聲波方程[15]：

流體中的非線性聲波傳播問題研究過程要更復(fù)雜一些，Earnshaw等發(fā)展了初期的非線性聲傳播理論[17-20]，F(xiàn)ubini, Fay, Blackstock等對(duì)拉氏坐標(biāo)下的非線性聲波方程組式(2)進(jìn)行求解。

其中：ρ和ρ0，p和p0分別為有擾動(dòng)和無擾動(dòng)情況下的介質(zhì)密度和壓強(qiáng)。通過求解獲得了相關(guān)條件下的解析解。Beyer則改寫(2)式，得到了氣、液體中的非線性聲波方程[18]：

式中：γ是氣體的比熱容比；B/A為液體的非線性參數(shù)，與Landau描述固體介質(zhì)非線性性質(zhì)的A、B無關(guān)。設(shè)是流體中的非線性系數(shù)[19]。利用泰勒展開，可將式(3)、(4)改寫為

比較式(1)和式(5)，可以發(fā)現(xiàn)，氣液固三種介質(zhì)中的非線性聲波方程形式上完全一樣，只是固體中的非線性系數(shù)β是流體中非線性系數(shù)的–2 倍[6,9](文獻(xiàn)[19]中非線性參數(shù)B/A為5，對(duì)應(yīng)=–3.5)。這也是一些文獻(xiàn)中二階攝動(dòng)近似解的常系數(shù)是1/8[9]，然而另一些文獻(xiàn)中是1/4的原因[1,5]。

2 一維非線性聲波方程的五種解法

2.1 經(jīng)典解析解[19-20]

Fubini, Fay等對(duì)流體中的非線性聲波方程(2)求解，分別得到近距離和遠(yuǎn)距離的相對(duì)振速解：

式中：v和v0分別為諧波和激勵(lì)信號(hào)的質(zhì)點(diǎn)振動(dòng)速度，即v/ v0表示諧波的相對(duì)(激勵(lì)信號(hào))振動(dòng)速度。但式(6)無法在沖擊波產(chǎn)生位置附近取得較好的結(jié)果，Blackstock利用橋函數(shù)：

將式(6)橋接起來，相對(duì)完滿地解決了問題。式(6)在許多文獻(xiàn)中均有詳細(xì)闡述，因此不再過多展開說明。

2.2 攝動(dòng)方法

式(7)被稱為非線性聲波方程的高階攝動(dòng)展開方程，一階攝動(dòng)展開方程(7.1)的解為。其他方程均為非齊次方程，非齊次項(xiàng)由低階次諧波組成，即低階諧波可以認(rèn)為是高階諧波的源。低階(如i和j)次諧波耦合形成和(i+j)諧波，也可合成為差諧(i–j)波。

高次諧波解的形式會(huì)隨著攝動(dòng)展開階次的增加變得越來越復(fù)雜，為了得到n階的高次諧波解，所有比n小的諧波都需要已知。

利用符號(hào)計(jì)算工具，并進(jìn)行簡(jiǎn)單的人為干預(yù)，可以得到高達(dá)8階的非齊次方程特解：

式(8)中： θ =ωt- k x ，因更高階的諧波解過于復(fù)雜，只將偶數(shù)次諧波里的二次諧波列出，參見附錄。從式(8)中可以看出，所有偶數(shù)階次的諧波解中均包含二次諧波解。理論上分析可知，得到高次攝動(dòng)解的階次越高，二次諧波的最終表達(dá)式會(huì)越接近真實(shí)解。非線性聲波方程(1)的解析解應(yīng)該是所有攝動(dòng)展開方程的解之和。

2.3 有限差分方法

一維非線性聲波方程式(1)的差分形式可以表示為

這里不會(huì)產(chǎn)生橫波，求解過程相對(duì)橫波激勵(lì)的非線性聲波[21]更加容易，此處不再詳述計(jì)算細(xì)節(jié)，但需要注意的是它的空間和時(shí)間步長(zhǎng)不能采用通常差分方程中的步長(zhǎng)。為適應(yīng)高次諧波的短波長(zhǎng)，計(jì)算采用的空間和時(shí)間步長(zhǎng)要遠(yuǎn)遠(yuǎn)小于1/8的波長(zhǎng)和周期。文中采用了1/160波長(zhǎng)的空間步長(zhǎng)。另外，非線性聲波方程的邊界條件問題也較為復(fù)雜，可采用相對(duì)較長(zhǎng)的介質(zhì)，若只研究聲波反射之前的傳播性質(zhì)，邊界條件問題就可以不作考慮。

2.4 有限元方法[6]

利用有限元方法，方程(1)可寫為

其中：M、K、B、F分別表示質(zhì)量矩陣、剛度矩陣、非線性剛度矩陣和力源，其他相關(guān)參數(shù)的物理意義參見前期的研究工作[6]。

通過計(jì)算可以得到非線性聲波方程的有限元數(shù)值聲場(chǎng)。

2.5 偽線性解

即不同位置(或位移空間梯度)處的質(zhì)點(diǎn)具有不同的聲傳播速度，這樣，非線性聲波方程的解可化簡(jiǎn)為

式(13)即為非線性聲波方程的偽線性解，但它不能直接獲得，仍需要部分?jǐn)?shù)值計(jì)算。在計(jì)算下一時(shí)刻的聲場(chǎng)時(shí)，將當(dāng)前時(shí)刻下的每個(gè)位置的聲速利用式(12)進(jìn)行估計(jì)，代入(13)即可計(jì)算下個(gè)時(shí)刻的聲波時(shí)域波形。

3 理論計(jì)算與實(shí)驗(yàn)結(jié)果的對(duì)比

為驗(yàn)證上述解的正確性，同時(shí)為了測(cè)量不同位置的非線性聲場(chǎng)的便利性，在水中開展非線性聲學(xué)實(shí)驗(yàn)，實(shí)驗(yàn)在Ritek-SNAP-5000系統(tǒng)上開展。信號(hào)發(fā)生器發(fā)射 12個(gè)周期的脈沖射頻信號(hào)給功率放大器放大后再傳給水中主頻為2 MHz的超聲換能器，換能器將電信號(hào)轉(zhuǎn)化為超聲信號(hào)在水中傳播，并由水中另一個(gè)寬帶超聲換能器接收后再傳入系統(tǒng)進(jìn)行分析與存儲(chǔ)。水的密度為1 g·cm-3, 聲速經(jīng)測(cè)量為1491 m·s-1。取文獻(xiàn)[19]中= -3.5或 β=7進(jìn)行數(shù)值計(jì)算。

有限差分?jǐn)?shù)值解(Finite Difference Time Domain, FDTD)和偽線性解(Pseudo Linear Solution,PSEU)的時(shí)域波形與實(shí)驗(yàn)信號(hào)(Experimental Signal,EXP)的對(duì)比如圖 1所示，這些信號(hào)均是在距離聲源 6 cm處得到的。理論的偽線性解為無限長(zhǎng)的穩(wěn)態(tài)解，圖中的波形是用了梯形窗進(jìn)行時(shí)域上截?cái)嗟慕Y(jié)果。

圖1 有限差分法和偽線性解得到的距聲源6 cm處非線性聲波位移信號(hào)與實(shí)驗(yàn)結(jié)果比較Fig.1 Comparison of the nonlinear ultrasonic displacement signals at 6 cm from source obtained by FDTD, PSEU and EXP

Blackstock橋函數(shù)提供了不同傳播位置的(相對(duì)聲源的)相對(duì)振速解(或歸一化解)，但是通過式(6.3)并不能得到時(shí)域波形信號(hào)；同時(shí)，僅利用8階次的諧波(包括基頻)信號(hào)，得到的攝動(dòng)時(shí)域波形與實(shí)驗(yàn)信號(hào)相差仍然較大；而通過有限元方法計(jì)算的時(shí)域波形與有限差分的結(jié)果幾乎完全重疊。這三種結(jié)果均未在圖1中列出。

有限差分(或有限元)及偽線性解的時(shí)域信號(hào)均與實(shí)驗(yàn)信號(hào)相似，都能展現(xiàn)非線性聲波信號(hào)的波形畸變性質(zhì)。正弦波隨傳播距離逐漸演化為尖銳的三角波。通過對(duì)時(shí)間求導(dǎo)，也可以得到如圖2所示的質(zhì)點(diǎn)振速的鋸齒波。

圖2 有限差分法和偽線性解得到的距聲源6 cm處非線性聲波質(zhì)點(diǎn)振速信號(hào)與實(shí)驗(yàn)結(jié)果比較Fig.2 Comparison of the nonlinear particle velocity signals at 6 cm from source obtained by FDTD, PSEU and EXP

根據(jù)式(8)，諧波幅度隨傳播距離增加而增大，即非線性聲波的畸變程度會(huì)加劇，如圖3所示。實(shí)線從下向上為 2、4、6 cm 處接收的非線性聲波，而虛線為對(duì)應(yīng)位置的FDTD計(jì)算信號(hào)。顯然，實(shí)驗(yàn)和理論計(jì)算的仿真信號(hào)正如式(8)預(yù)期，畸變程度逐漸變大。非線性聲波波前將最終演化為沖擊波。

圖3 有限差分法得到的距離聲源2、4、6 cm處接收的非線性超聲信號(hào)與實(shí)驗(yàn)結(jié)果比較Fig.3 Comparison of the nonlinear ultrasonic waves measured at 2, 4, 6 cm from source obtained by FDTD and EXP

利用信號(hào)處理可得到不同聲傳播距離下的基頻信號(hào)和二次諧波相對(duì)(位移)幅度，如圖 4所示。EXP表示實(shí)驗(yàn)信號(hào)，P2、P4、P6、P8為相應(yīng)階次的高階攝動(dòng)解(Perturbation, PERT)，F(xiàn)IT (Fitting Solution)代表前期工作中的擬合解[6]，F(xiàn)DM為對(duì)有限差分信號(hào)的處理結(jié)果，PSEU為偽線性解，BLKS為Blackstock橋函數(shù)的解。需要說明的是，圖4中是位移的相對(duì)幅度，對(duì)于二次諧波，考慮線性解和式(8.1)，其值為，而文獻(xiàn)[18, 20]中振速的相對(duì)幅度為，它們之間存在兩倍的差異[6,18-19]。

基頻信號(hào)的相對(duì)位移幅度位于圖4中的上方位置，在接近聲源位置時(shí)，因?yàn)樯形从须S傳播距離積累的諧波出現(xiàn)，即基頻信號(hào)還未向高次諧波轉(zhuǎn)化，同時(shí)也無損耗、擴(kuò)散衰減發(fā)生，相對(duì)幅度從1開始隨傳播距離減小，標(biāo)志著基頻信號(hào)隨傳播距離開始向諧波傳遞能量。在遠(yuǎn)離聲源位置處，考慮到擴(kuò)散衰減和吸收衰減，實(shí)驗(yàn)信號(hào)在所有理論計(jì)算結(jié)果中的相對(duì)幅度最小。

圖4 不同方法計(jì)算和實(shí)驗(yàn)得到的不同傳播位置的基頻(上)和二次諧波(下)相對(duì)(聲源)位移幅度Fig.4 Relative displacement amplitudes of fundamental(upper) and second harmonic (lower) at different propagation distances obtained by different calculation methods and experiment

二次諧波的相對(duì)幅度位于圖4下方，同樣由于傳播距離和衰減的原因，其相對(duì)幅度從 0開始增加，增加到一定幅度后逐漸減小，實(shí)驗(yàn)的二次諧波相對(duì)幅度在遠(yuǎn)距離處也同樣比幾種理論計(jì)算值小。

4 討 論

對(duì)于攝動(dòng)法，其一階攝動(dòng)展開方程的解為線性方程的解，信號(hào)的振幅為常數(shù)，相對(duì)幅度一直為1。而在3階、5階、7階等奇數(shù)階次攝動(dòng)展開方程的解中均有一次波出現(xiàn)，比如三次諧波攝動(dòng)展開方程式(7.3)中，二次諧波和一次耦合成三次的同時(shí)，還能耦合出一次諧波(見式(8.2))。文中稱高次攝動(dòng)展開方程得到的一次波為一次諧波。偶數(shù)次諧波解中不包含一次波。將所有攝動(dòng)展開方程的一次波疊加，才能得到基頻信號(hào)隨傳播距離的變化規(guī)律。

對(duì)于二次諧波，P2為只考慮二次諧波的攝動(dòng)解，它與傳播距離成線性關(guān)系，僅能在較接近聲源時(shí)才與實(shí)驗(yàn)結(jié)果吻合。而P4、P6、P8為分別考慮了對(duì)應(yīng)的高次諧波中的二次諧波成分(所有偶數(shù)階的攝動(dòng)展開方程中均能耦合二次諧波)，能在更寬泛的測(cè)量條件下比低階近似解更貼近實(shí)驗(yàn)結(jié)果。但圖中P8相對(duì)P6似乎并未提高精度，這在數(shù)學(xué)上也能夠理解，因?yàn)閷⑺须A次的攝動(dòng)解累加起來是非線性聲波方程的解，而這里討論的只是二次諧波。但總體上考慮的階次越高，這些解越能接近實(shí)驗(yàn)結(jié)果。

攝動(dòng)解(也包括其他解)在傳播距離較遠(yuǎn)時(shí)與實(shí)驗(yàn)結(jié)果的差異仍相對(duì)比較明顯，一方面是由于在遠(yuǎn)距離處應(yīng)考慮更多高階的諧波解，另一方面，遠(yuǎn)距離處的吸收衰減和擴(kuò)散也不能忽略。另外，源信號(hào)的形狀也需要考慮，攝動(dòng)解(和偽線性解)是連續(xù)波，而實(shí)驗(yàn)與數(shù)值計(jì)算中采用的是脈沖串。同時(shí)，文中未涉及更高階次的材料非線性系數(shù)，這些也會(huì)對(duì)結(jié)果造成一定影響，其他相關(guān)文獻(xiàn)中[9]也有類似的實(shí)驗(yàn)結(jié)果，呈現(xiàn)了非線性問題的復(fù)雜性。

就二次諧波的相對(duì)幅度而言，幾種理論解的精度在測(cè)量的距離范圍內(nèi)均相對(duì)較為理想，與實(shí)驗(yàn)結(jié)果的吻合程度也較高，這是因?yàn)槲恼驴紤]的是介質(zhì)非線性，諧波隨傳播距離產(chǎn)生并與其他階次的諧波(也包括一次波)相互耦合，最終形成與聲波幅度、傳播距離等許多物理參量相關(guān)的極為復(fù)雜的數(shù)學(xué)關(guān)系。在傳播距離不是非常遠(yuǎn)的情況下，這些解具有一定的精度，但在傳播距離較遠(yuǎn)時(shí)，這些解之間以及它們與實(shí)驗(yàn)值之間的誤差仍需要進(jìn)一步的研究和分析。

文中介紹的五種非線性聲波方程解法之中，有限元和有限差分都是數(shù)值解，攝動(dòng)方法是一種近似的解析解，雖然偽線性解是形式上的解析解，但只能寫出一種隱函數(shù)式，波形計(jì)算需要一定程度的數(shù)值計(jì)算。Blackstock橋函數(shù)能描述不同傳播距離下所有高次諧波的質(zhì)點(diǎn)相對(duì)振速幅度，然而并不能得到傳播過程中波形的變化趨勢(shì)。

攝動(dòng)解限制于能求解出的階次，盡管已經(jīng)得到了高達(dá)八階的諧波解，但合成的非線性時(shí)域波形的畸變程度仍遠(yuǎn)遠(yuǎn)不夠。但它的優(yōu)點(diǎn)也是其他方法不具有的，即定量分析。許多實(shí)驗(yàn)仍運(yùn)用二階攝動(dòng)近似解(式(8.1))表征介質(zhì)的非線性系數(shù)[5,8-9]，但這不能在傳播距離較遠(yuǎn)(或介質(zhì)較長(zhǎng))和激勵(lì)信號(hào)幅度較高時(shí)的情況下與實(shí)驗(yàn)結(jié)果取得較好的一致性[9]。非線性聲波的傳播距離越遠(yuǎn)，諧波就越復(fù)雜，更高階諧波解中的二次諧波的成分不能被二階攝動(dòng)近似解囊括，它的局限性就越能暴露出來。顯然，文中的8階非攝動(dòng)解已能大大擴(kuò)展非線性系數(shù)實(shí)驗(yàn)的測(cè)量范圍，當(dāng)然也能提高其測(cè)量精度。

有限差分和有限元方法在意義上略有不同，然而它們的結(jié)果幾乎一致，直觀觀察時(shí)域聲場(chǎng)的演化過程比較容易，傳播性質(zhì)也可通過信號(hào)處理較容易地得到，但卻不利于量化分析。偽線性解實(shí)際是一種半解析方法，在它的非線性解中，聲傳播速度隨時(shí)間和傳播位置(或質(zhì)點(diǎn)振速的空間梯度)的變化而變化。

5 結(jié) 論

文中介紹了求解一維非線性聲波方程的五種方法，比較分析了它們的優(yōu)劣和適用范圍。有限差分法和有限元方法均能直觀地觀察時(shí)域波形隨傳播距離的演化過程。而偽線性解和攝動(dòng)法得到的是無限長(zhǎng)的穩(wěn)態(tài)解，Blackstock橋函數(shù)不能得到時(shí)域波形。

在分析高次諧波相對(duì)幅度隨傳播距離(或激勵(lì)信號(hào)幅度)的傳播性質(zhì)上，除攝動(dòng)解外，其他四種方法的計(jì)算結(jié)果都與實(shí)驗(yàn)結(jié)果有較好的一致性。然而，由于攝動(dòng)解被廣泛應(yīng)用在非線性系數(shù)的定量計(jì)算上，二階攝動(dòng)近似解僅在接近聲源的區(qū)域內(nèi)適用性高，更高階次諧波攝動(dòng)解中的二次諧波可擴(kuò)展非線性系數(shù)表征的測(cè)量范圍，同時(shí)也能提高它的測(cè)量精度。

與實(shí)驗(yàn)結(jié)果的對(duì)比驗(yàn)證了幾種非線性聲波方程解的有效性，預(yù)期可為測(cè)量流體和固體介質(zhì)非線性系數(shù)提供更精確的方法。

附錄