非線性波動方程的新數值迭代方法*

曹娜 陳時 曹輝 王成會 劉航

(陜西師范大學物理學與信息技術學院, 陜西省超聲學重點實驗室, 西安 710119)

提出了一種新的求解非線性波動方程的數值迭代法, 它是一種半解析的方法.與完全的數值計算方法(如有限元、有限差分法)相比, 這種迭代法的解具有非常清晰的物理含義, 即它的解是各階諧波的組合.與微擾法相比, 它能夠考慮各階諧波的相互作用, 且能夠滿足能量守恒定律.用它研究了非線性聲波在液體中的傳播性質, 結果表明, 在微擾法適用的聲強范圍內迭代法也適用, 在微擾法不適用的一個較寬的聲強范圍內迭代法依然適用.

1 引 言

非線性聲學是聲學的一個重要分支, 當系統的聲強比較強時就會產生各種非線性效應, 如諧波產生、沖擊波形成、聲輻射力的出現等.非線性聲學在某些領域有著重要的運用, 如高聲強聚焦超聲[1?4]、超聲懸浮[5?7]、聲空化[8]、聲諧波成像[9,10]、參量發射陣[11,12]等.在這些領域中聲波的非線性方程的求解是非常重要的.

關于非線性聲學系統波動方程的求解, 現階段用到的方法一般包括：1)完全的數值計算方法, 如有限元和有限差分法[13?17].這類方法得到解的物理意義并不明確, 很難揭示非線性作用的物理本質, 而且在很多情況下還會引起數值發散問題, 并非適用于全部的非線性問題.2)嚴格的解析方法[18?21].這種方法只能處理極少數系統的非線性聲學問題,如理想流體中非線性聲波的傳播.3)微擾法[22].它的優點是方法簡單和解的物理意義清晰, 但是只適合處理低聲強時的非線性效應.且它只考慮低階諧波對高階諧波的作用, 而忽略其反作用, 因此并不滿足能量守恒定律.

對于聲學非線性方程的求解問題, 本文提出了一種新的、半解析的數值迭代方法.它是在頻域內把聲場展開為傅里葉級數的形式, 實現時間變量和空間坐標的分離.然后根據計算精度的具體需求,截斷高頻諧波而實現方程的求解.它的解具有非常清晰的物理意義, 即是各階諧波的組合.經過研究發現, 在微擾法適用的聲強范圍內, 本文提出的方法也是適用, 且滿足能量守恒定律(無耗散的系統).在微擾法不適用的一個較寬的聲強范圍內, 迭代法依然適用且滿足能量守恒定律(無耗散的系統).只是在極高聲強的情況下, 本文提出的方法才不適用.

2 理論方法

2.1 非線性波動方程的迭代數值方法

在拉格朗日坐標系下, 黏性液體中一維非線性聲波的位移滿足下式[23]：

式中 β 為液體的非線性系數; c0是靜態時(即不存在聲波時)液體的聲速;其中 μ 為液體的體積黏滯系數, ρ0是靜態時(即不存在聲波時)液體的密度.下標中逗號后的坐標x和時間t表示對它們求偏導數.令 μ =0 , 則方程(1)退化為理想液體中一維非線性聲波的波動方程.

在許多情況下(如求解穩態問題), 把u中的變量t和x分離開來是有利的, 一般情況下u可以表示為(可以稱為頻域內的傅里葉級數展開)：

其中i是一個虛單位, ω 表示波的角頻率; A0為一個實數場變量, A0/2 表示聲波的“直流”部分;An(n≥1)為第n階諧波的復數場變量(即復振幅),Anexp(inωt)的實部是第n階諧波真實的位移,是 An的復數共軛場變量.注意 An和中已經不包含時間變量t, 它們只是空間坐標x的函數.

一般情況下高階諧波是比較弱的, 根據計算精度的需要可以忽略掉某些高階諧波.為了簡化理論的敘述, 本文只考慮階數小于或者等于N(N≤6)的各階諧波, 忽略掉其他的高階諧波, 稱之為N階近似.

將方程(2)代入方程(1)中, 因為有相同時間因子( e xp(inωt) , n =0,±1,±2,··· )的項之和必須為零, 所以可以得到下面的方程：

注意此處只給出了場變量的方程, 共軛場的方程并沒有列出來, 只要對方程(3)—(9)取復數共軛就可以得到共軛場的方程, 因此方程(3)—(9)是完備的.

方程(3)—(9)是一組耦合的非線性方程, 直接求解它們是很困難的.本文提出求解它們的一種新的簡單迭代方法.用 A(m)和 A?(m)( m ≥0 )表示第m次迭代計算得到的場量.在第m次迭代計算中,采用了如下方法：方程(3)—(9)等號左邊的場量取為 A(m), 右邊的場量取為 A(m-1)和 A?(m-1).在第m次迭代計算中, 用到如下的方程：

分別用 A(m-1)和 A?(m-1)替換F中的A和 A ? , 得到的結果就是 F(m-1).方程(17)—(23)是一組非耦合的方程, 因此可以分別獨立地計算出.這意味著當涉及到更多的高階諧波時, 計算量不會急劇地增加.

2.2 本文迭代方法的具體運用過程

用迭代方法研究非線性聲波在黏性液體中的傳播問題.現設在 x =0 處有一列平面波朝x正向傳播, 其為入射聲波, 且聲場可以表示為

其中 Bi是一個已知量.入射波的能流密度 Pi可以表示為

本文需要計算在 x =L 處出射的各階諧波的聲場.它們可以表示為

其中 Bon就是要計算的量.出射波的能流密度Pon可以表示為

為了求解在 [ 0 ,L] 坐標間隔內的非線性聲場,用有限差分法來求解方程(17)—(23).在迭代計算中, 令 A(0)=0 和 A?(0)=0 , 非 零的 A(m)和A?(m)由邊界激勵條件產生.用到的邊界條件是：在x=0和 x =L 兩個端點處, 各階諧波的位移和垂直應力都是連續的.

3 數值計算和討論

通過數值計算分析了非線性聲波在液體中的傳播性質, 得到了本文提出的新數值方法的適用范圍, 并證明了其有效性.在所有的計算中, 如果沒有特別說明, 那么用到的參數是：c0= 1.5 ×103m/s, ρ0= 103kg/m3, P0= 1.01 × 105Pa,μ=1× 10—3Pa/s, L =0.05 m, ω = 5 × 106rad/s,β=3.5.如果文中或圖中對某個參數有特別說明,那么該參數就替換為特殊說明處的數據.

圖1顯示了非線性聲波在理想液體( μ =0 )中傳播時能量守恒的破壞程度 Ed、二階諧波的相對能流 P2/Pi和三次諧波的相對能流 P3/Pi隨入射能流 Pi的變化情況.能量守恒的破壞程度 Ed定義為所有出射能流密度與入射能流密度的相對差值,即越小能量守恒越能保證, 它越大能量守恒定律破壞程度越大.實線對應著本文提出的迭代法的情況, 虛線對應著微擾法的情況.

由圖1(a)可見, 能量守恒的破壞程度 Ed隨入射波能流 Pi的增大而增大.其他參數不變, 入射波能流相同時, 明顯可以看出本文提出的新數值迭代方法得到的能量守恒的破壞程度比微擾法得到的能量守恒的破壞程度小很多.當入射波能流是1.5×107J/(m2·s)時, 迭代法的破壞程度是 7.7 ×10—3, 微擾法的破壞程度是 0 .506 ; 當入射能流是5.5×107J/(m2·s)時, 迭代法的破壞程度 1.92 ×10—2, 微擾法的破壞程度是 2 .268 .

由圖1(b)和圖1(c)可知, 二次諧波的相對能流 P2/Pi和三次諧波的相對能流 P3/Pi均隨入射能流 Pi的增加而增加.當入射波能流小于1.5×107J/(m2·s)時, 迭代法和微擾法得到的P2/Pi和 P3/Pi的值幾乎相同; 當入射波能流大于5.5×107J/(m2·s) 時, 迭代法和微擾法得到的 P2/Pi的值分別等于 0 .6295 和 1 .4741 及 P3/Pi的值分別大于0.0421和 0 .4847 .

從圖1還可以看出, 當入射波能流小于1.5×107J/(m2·s)時, 兩種方法(本文提出的新數值迭代法和微擾法)均適用; 當入射波能流小于5.5×107J/(m2·s)且大于 1 .5× 107J/(m2·s)時, 迭代法適用, 微擾法不適用; 當入射波能流大于5.5×107J/(m2·s)時, 兩種方法均不適用.

圖2顯示了聲波在理想液體( μ =0 )中傳播時各階諧波的相對能流 Pn/Pi(n=1,2,···,6) 隨入射波能流 Pi的變化情況.從圖2可見, 各階諧波的能流隨入射波能流的變化趨勢相同, 均隨入射波能流的增加而增加, 但基頻波的能流明顯大于其他高階諧波的能流.從圖2還可以看出, 基頻波的相對能流隨入射波能流的增加而減小, 高階諧波的相對能流隨入射波能流的增加而增加.

圖2 各階諧波的相對能流隨入射波能流的變化Fig.2.Relation between relative energy of each order of harmonics and incident wave energy.

圖2顯示的結果是與文獻[20, 24]中基頻波能流與各高階諧波能流之間的關系相一致.從圖2可以明顯得到, 聲波在介質中傳播時, 基頻波能流向各高階諧波傳遞, 這同時也解釋了, 基頻波的相對能流隨入射波能流的增加而減少, 但各高階諧波的相對能流卻隨入射波能流的增加而增加.

圖3顯示了聲波在理想液體( μ =0 )中傳播時入射波能流不同的情況下二階諧波的相對能流 P2/Pi隨迭代次數m的變化關系.m表示迭代次數.由圖3可知, 當入射波能流是0.4687×107J/(m2·s), 迭代次數大于等于3時, 二次諧波的相對能流收斂于 0 .145 ; 當入射波能流是2.7×107J/(m2·s), 迭代次數大于等于4時, 二次諧波的相對能流收斂于 0 .515 ; 當入射波能流是5.418×107J/(m2·s), 迭代次數大于等于 7時, 二次諧波的相對能流收斂于 0 .639 .由圖3可見, 本文提出的新數值迭代方法具有收斂性, 且入射聲強越大, 達到收斂的迭代次數需要越多.

圖3 入射波能流不同時二階諧波的相對能流隨迭代次數的變化Fig.3.Relative energy flow of the second harmonic varies with the number of iterations under different incident wave energy flow.

圖4 黏度不同的情況下, 各階諧波的能流隨入射波能流的變化(圖中曲線的黏度分別是 1 ×10-6 , 6 ×10-1 和 1 0×10-1 Pa/s,箭頭表示黏度減小的方向)Fig.4.Relation of relative energy flow of each order of harmonics with incident wave energy flow under different viscosity.Viscos?ites for different curves are 1 ×10-6 , 6 ×10-1 , 1 0×10-1 Pa/s, respectively.Arrows indicate the direction of decreasing viscosity.

圖5為聲波在黏性液體( μ =0 )中傳播時, 黏度不同的情況下各階諧波的相對能流 Pn/Pi(n =1, 2)隨角頻率的變化.從圖5(a)和圖5(b)可見：1)基頻波的相對能流隨角頻率的增加而減少, 二次諧波的相對能流隨角頻率的增加而增加; 2)隨黏度的增加基頻波和二次諧波的相對能流均減小,但基頻波的影響程度比二次諧波的影響程度更大一些.

圖5 黏度不同的情況下, 各階諧波的相對能流隨角頻率的變化(圖中曲線的黏度分別是 1 ×10-6 和 1 ×10-1 Pa/s,箭頭表示黏度減小的方向)Fig.5.Relative energy flow varies with angular frequency under different visco?sity.Viscosites for different curves are 1×10-6and 1 ×10-1 Pa/s, respectively.Arrows indic?ate the direction of decreasing viscosity.

4 結 論

用一種新的數值計算方法研究了聲波在液體中的傳播特性.得到的主要結論如下：

2)聲波在介質中傳播時, 基頻波的相對能流隨入射波能流(或角頻率)的增加而減少, 但各高階諧波的相對能流卻隨入射波能流(或角頻率)的增加而增加;

3)聲波的各階諧波能流(或相對能流)均隨黏度的增加而減小;

4)在相同的參數下, 本文提出的研究計算方法比微擾法能更好地保證在研究聲傳播的過程中能量守恒, 且當涉及到更多的高次諧波時, 用本文提出的數值方法計算時間不會急劇增加.