非線性聲波方程的幾種解法對比

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

(1. 貴州理工學院，貴州 貴陽 550003；2. 貴州省醫工交叉工程研究中心，貴州 貴陽 550003；3. 近地面探測技術重點實驗室，江蘇 無錫 214035)

0 引 言

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

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

目前，實驗表征材料非線性系數時常用二階攝動解進行計算[5,8-9]，但二階攝動近似解只能在近場或小振幅激勵信號才適用，當樣品不太小或激勵源不太弱時，得到的實驗結果會有較大誤差[9]。

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

1 氣液固介質中的非線性聲波方程

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

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

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

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

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

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

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

2.1 經典解析解[19-20]

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

式中：v和v0分別為諧波和激勵信號的質點振動速度，即v/ v0表示諧波的相對(激勵信號)振動速度。但式(6)無法在沖擊波產生位置附近取得較好的結果，Blackstock利用橋函數：

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

2.2 攝動方法

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

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

利用符號計算工具，并進行簡單的人為干預，可以得到高達8階的非齊次方程特解：

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

2.3 有限差分方法

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

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

2.4 有限元方法[6]

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

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

通過計算可以得到非線性聲波方程的有限元數值聲場。

2.5 偽線性解

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

式(13)即為非線性聲波方程的偽線性解，但它不能直接獲得，仍需要部分數值計算。在計算下一時刻的聲場時，將當前時刻下的每個位置的聲速利用式(12)進行估計，代入(13)即可計算下個時刻的聲波時域波形。

3 理論計算與實驗結果的對比

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

有限差分數值解(Finite Difference Time Domain, FDTD)和偽線性解(Pseudo Linear Solution,PSEU)的時域波形與實驗信號(Experimental Signal,EXP)的對比如圖 1所示，這些信號均是在距離聲源 6 cm處得到的。理論的偽線性解為無限長的穩態解，圖中的波形是用了梯形窗進行時域上截斷的結果。

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

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

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

圖2 有限差分法和偽線性解得到的距聲源6 cm處非線性聲波質點振速信號與實驗結果比較Fig.2 Comparison of the nonlinear particle velocity signals at 6 cm from source obtained by FDTD, PSEU and EXP

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

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

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

基頻信號的相對位移幅度位于圖4中的上方位置，在接近聲源位置時，因為尚未有隨傳播距離積累的諧波出現，即基頻信號還未向高次諧波轉化，同時也無損耗、擴散衰減發生，相對幅度從1開始隨傳播距離減小，標志著基頻信號隨傳播距離開始向諧波傳遞能量。在遠離聲源位置處，考慮到擴散衰減和吸收衰減，實驗信號在所有理論計算結果中的相對幅度最小。

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

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

4 討 論

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

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

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

就二次諧波的相對幅度而言，幾種理論解的精度在測量的距離范圍內均相對較為理想，與實驗結果的吻合程度也較高，這是因為文章考慮的是介質非線性，諧波隨傳播距離產生并與其他階次的諧波(也包括一次波)相互耦合，最終形成與聲波幅度、傳播距離等許多物理參量相關的極為復雜的數學關系。在傳播距離不是非常遠的情況下，這些解具有一定的精度，但在傳播距離較遠時，這些解之間以及它們與實驗值之間的誤差仍需要進一步的研究和分析。

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

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

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

5 結 論

文中介紹了求解一維非線性聲波方程的五種方法，比較分析了它們的優劣和適用范圍。有限差分法和有限元方法均能直觀地觀察時域波形隨傳播距離的演化過程。而偽線性解和攝動法得到的是無限長的穩態解，Blackstock橋函數不能得到時域波形。

在分析高次諧波相對幅度隨傳播距離(或激勵信號幅度)的傳播性質上，除攝動解外，其他四種方法的計算結果都與實驗結果有較好的一致性。然而，由于攝動解被廣泛應用在非線性系數的定量計算上，二階攝動近似解僅在接近聲源的區域內適用性高，更高階次諧波攝動解中的二次諧波可擴展非線性系數表征的測量范圍，同時也能提高它的測量精度。

與實驗結果的對比驗證了幾種非線性聲波方程解的有效性，預期可為測量流體和固體介質非線性系數提供更精確的方法。

附錄