非線性波動(dòng)方程的新數(shù)值迭代方法*

曹娜 陳時(shí) 曹輝 王成會(huì) 劉航

(陜西師范大學(xué)物理學(xué)與信息技術(shù)學(xué)院, 陜西省超聲學(xué)重點(diǎn)實(shí)驗(yàn)室, 西安 710119)

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

1 引 言

非線性聲學(xué)是聲學(xué)的一個(gè)重要分支, 當(dāng)系統(tǒng)的聲強(qiáng)比較強(qiáng)時(shí)就會(huì)產(chǎn)生各種非線性效應(yīng), 如諧波產(chǎn)生、沖擊波形成、聲輻射力的出現(xiàn)等.非線性聲學(xué)在某些領(lǐng)域有著重要的運(yùn)用, 如高聲強(qiáng)聚焦超聲[1?4]、超聲懸浮[5?7]、聲空化[8]、聲諧波成像[9,10]、參量發(fā)射陣[11,12]等.在這些領(lǐng)域中聲波的非線性方程的求解是非常重要的.

關(guān)于非線性聲學(xué)系統(tǒng)波動(dòng)方程的求解, 現(xiàn)階段用到的方法一般包括：1)完全的數(shù)值計(jì)算方法, 如有限元和有限差分法[13?17].這類(lèi)方法得到解的物理意義并不明確, 很難揭示非線性作用的物理本質(zhì), 而且在很多情況下還會(huì)引起數(shù)值發(fā)散問(wèn)題, 并非適用于全部的非線性問(wèn)題.2)嚴(yán)格的解析方法[18?21].這種方法只能處理極少數(shù)系統(tǒng)的非線性聲學(xué)問(wèn)題,如理想流體中非線性聲波的傳播.3)微擾法[22].它的優(yōu)點(diǎn)是方法簡(jiǎn)單和解的物理意義清晰, 但是只適合處理低聲強(qiáng)時(shí)的非線性效應(yīng).且它只考慮低階諧波對(duì)高階諧波的作用, 而忽略其反作用, 因此并不滿(mǎn)足能量守恒定律.

對(duì)于聲學(xué)非線性方程的求解問(wèn)題, 本文提出了一種新的、半解析的數(shù)值迭代方法.它是在頻域內(nèi)把聲場(chǎng)展開(kāi)為傅里葉級(jí)數(shù)的形式, 實(shí)現(xiàn)時(shí)間變量和空間坐標(biāo)的分離.然后根據(jù)計(jì)算精度的具體需求,截?cái)喔哳l諧波而實(shí)現(xiàn)方程的求解.它的解具有非常清晰的物理意義, 即是各階諧波的組合.經(jīng)過(guò)研究發(fā)現(xiàn), 在微擾法適用的聲強(qiáng)范圍內(nèi), 本文提出的方法也是適用, 且滿(mǎn)足能量守恒定律(無(wú)耗散的系統(tǒng)).在微擾法不適用的一個(gè)較寬的聲強(qiáng)范圍內(nèi), 迭代法依然適用且滿(mǎn)足能量守恒定律(無(wú)耗散的系統(tǒng)).只是在極高聲強(qiáng)的情況下, 本文提出的方法才不適用.

2 理論方法

2.1 非線性波動(dòng)方程的迭代數(shù)值方法

在拉格朗日坐標(biāo)系下, 黏性液體中一維非線性聲波的位移滿(mǎn)足下式[23]：

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

在許多情況下(如求解穩(wěn)態(tài)問(wèn)題), 把u中的變量t和x分離開(kāi)來(lái)是有利的, 一般情況下u可以表示為(可以稱(chēng)為頻域內(nèi)的傅里葉級(jí)數(shù)展開(kāi))：

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

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

將方程(2)代入方程(1)中, 因?yàn)橛邢嗤瑫r(shí)間因子( e xp(inωt) , n =0,±1,±2,··· )的項(xiàng)之和必須為零, 所以可以得到下面的方程：

注意此處只給出了場(chǎng)變量的方程, 共軛場(chǎng)的方程并沒(méi)有列出來(lái), 只要對(duì)方程(3)—(9)取復(fù)數(shù)共軛就可以得到共軛場(chǎng)的方程, 因此方程(3)—(9)是完備的.

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

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

2.2 本文迭代方法的具體運(yùn)用過(guò)程

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

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

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

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

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

3 數(shù)值計(jì)算和討論

通過(guò)數(shù)值計(jì)算分析了非線性聲波在液體中的傳播性質(zhì), 得到了本文提出的新數(shù)值方法的適用范圍, 并證明了其有效性.在所有的計(jì)算中, 如果沒(méi)有特別說(shuō)明, 那么用到的參數(shù)是：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.如果文中或圖中對(duì)某個(gè)參數(shù)有特別說(shuō)明,那么該參數(shù)就替換為特殊說(shuō)明處的數(shù)據(jù).

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

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

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

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

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

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

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

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

圖3 入射波能流不同時(shí)二階諧波的相對(duì)能流隨迭代次數(shù)的變化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 )中傳播時(shí), 黏度不同的情況下各階諧波的相對(duì)能流 Pn/Pi(n =1, 2)隨角頻率的變化.從圖5(a)和圖5(b)可見(jiàn)：1)基頻波的相對(duì)能流隨角頻率的增加而減少, 二次諧波的相對(duì)能流隨角頻率的增加而增加; 2)隨黏度的增加基頻波和二次諧波的相對(duì)能流均減小,但基頻波的影響程度比二次諧波的影響程度更大一些.

圖5 黏度不同的情況下, 各階諧波的相對(duì)能流隨角頻率的變化(圖中曲線的黏度分別是 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 結(jié) 論

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

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

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

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