水下聲波傳播過程的時域間斷Galerkin有限元方法

李瑞敏, 郭 攀*, 張景飛, 武文華, 王 飛

(1.鄭州大學 力學與安全工程學院，鄭州 450000；2.大連理工大學 運載工程與力學學部 工程力學系，工業裝備結構分析國家重點實驗室，大連 116024)

1 引 言

由于聲壓力波水下良好的傳播性能，聲波技術廣泛應用于海洋環境監測、海洋軍事、海洋航運、海洋科學等領域[1]。尤其是在水下發射、水下爆炸、水下爆破和水下氣泡射流等問題中，近場強間斷、高梯度復雜多變的沖擊聲壓力波往往引起結構和艦船的毀傷[2]，使得沖擊條件下聲壓力場結構的準確描述問題成為海洋工程和軍事等領域一直以來的研究熱點[3,4]；另外，梯度變化聲速的海洋聲波傳播問題一直以來是海洋聲學研究的重點[5-7]；復雜的海底地形邊界條件下，海底地質勘探、障礙物的探測識別和目標定位等技術領域的研究近年來從未間斷[8-12]。不難發現，這些問題往往伴隨著空間或時域高梯度、強間斷變化的特征。

目前，隨著計算機軟硬件的發展，在聲場的研究領域出現了諸多數值計算方法，如有限體積方法[13]、有限元法、有限差分法[14]、邊界元法[15]以及相關的混合方法等。其中有限元方法在工程領域中應用最多。以往研究表明，在模擬沖擊作用下間斷壓力波的傳播過程時，通常的時域連續Galerkin有限元方法如Newmark法會在壓力波波后出現虛假的數值振蕩，誤差的累積使得問題求解的精度大幅降低[16]。時域間斷Galerkin有限元方法 DGFEM(Discontinuous Galerkin Finite Element Method)方法在沖擊問題中體現了區別于傳統有限元法的優越性，其可以自動引入數值耗散，濾去虛假的高階模式和數值振蕩。Li等[17,18]在原有DGFEM基礎上，分別采用三次Hermit插值函數和線性插值函數對基本未知量及其導數進行插值，將基本未知量及其導數解耦，大幅提高了計算效率。但是在高頻振動的波傳播問題仿真中，這類方法在波前仍有較小的數值振蕩現象，Guo等[19]針對這種現象進行了進一步的研究，引入高頻波動比較敏感的比例剛度阻尼，并將改進的DGFEM運用到了諸多領域，在層合壓電材料沖擊問題、非傅里葉熱彈性問題及彈性固體沖擊問題等算例的求解中均表現出了該方法的優越性[20-22]。

本文在以往研究基礎上，基于聲壓波動傳播方程，構建出水下聲波傳播問題的改進DGFEM(以下簡稱DGFEM)，針對幾類典型的水下高梯度和強間斷問題進行了探討。

2 聲傳播控制方程及時域間Galerkin有限元方法

聲傳播方程為

(1)

式中p為流體壓力，c0為流體中的聲速。

對于聲波傳播方程，流體內單元分布可表示為

(2)

式中mf為流體單元的節點數，pe為單元的結構點壓力向量，Ni為對應節點i的插值函數。

(3)

(4)

(5)

考慮到聲源和水中阻尼的作用，式(5)可以略去下標f，進而推廣為

(6)

式中C為阻尼矩陣，Q(t)為各個節點上的聲源向量。聲波傳播方程的時域可以離散為

0<…

對任一時刻tn，其時刻處的階躍函數可以表示為

(7)

(8)

(9)

(10)

基本未知函數的導數v在時間步內采用線性插值：

(11)

省略上標，式(11)可簡化為

v=vnλ1+vn + 1λ2

(12)

(13)

式(13)分別對基本未知量及其時域一階導數進行變分，可得

(14)

式(14)解耦后，可得

(15)

改進的時域間斷Galerkin有限元方法是為了消除DGFEM方法波前虛假數值振蕩。其基于高頻振蕩敏感于比例剛度阻尼和低頻振蕩敏感于比例質量阻尼的原理，在式(15)引入比例剛度阻尼Ck=βK，對原有DGFEM方法進行進一步發展。

其中β為阻尼系數:

β=2ξ/ωn

(16)

ξ為結構的阻尼比，ωn為結構的固有頻率。

式(15)可以改寫為

(17)

3 算例及結果分析

算例1選取結構長度為50 m，截面積為0.5 m2的一維流體區域，如圖1所示。在區域最左端為完全反射邊界條件，最右端受到脈沖階躍荷載的作用，如圖2所示。流體介質密度為1000 kg/m3，聲波在該流體中的傳播速度為1500 m/s。結構劃分為100個單元，時間步為0.000 1 s，分別采用 Newmark 方法與改進的DGFEM進行模擬。

圖1 流體區域受沖擊作用模型

Fig.1 Model of fluid region by impact action

圖2 外載階躍荷載

圖3給出了傳統的Newmark方法與改進的DGFEM方法在0.005 s，0.01 s和0.02 s時水中壓力波分布的情況。可以看出，壓力波從右向左傳播，傳統的Newmark法計算得到的壓力波在波后和波前均出現了虛假的數值振蕩，而且從波前和波后數值振蕩的區域分布來看，這類數值振蕩隨著時間增長數值振蕩范圍不斷增加，誤差累積不斷增大，這勢必降低問題求解的準確性和穩定性。而改進的DGFEM方法可以很好地濾去不同時刻波前和波后虛假數值振蕩，相對較好地保證了這類間斷聲波傳播問題求解的精度和穩定性。

算例2為了驗證新方法在復雜邊界問題求解中的可靠性，考慮到水下障礙物的聲探測技術，設計了如圖4所示的模型。隨機分布的孔表示水中的障礙物，設定為全反射邊界。在數值模擬中，流體的密度為1000 kg/m3，聲音在流體中的傳播速度為1500 m/s，模擬區域的尺寸為20 m×20 m，將模型劃分為1681個單元，時間步長為1.0×10-4s。聲源位于點O(10.8 m,10.8 m)，聲壓力荷載如圖5所示。

圖3 不同時刻下壓力波傳播情況

圖4 含障礙物二維聲傳播模型網格劃分

圖5 外載階躍荷載

圖6和圖7是在t=0.002 s，0.005 s，0.008 s，0.01 s四個不同時刻，DGFEM和Newmark方法所模擬的流體區域中聲壓分布情況。可以看出，聲源激發聲波后，向周圍傳播，當它們接觸障礙物時，聲波發生反射并繼續傳播。當采用Newmark方法進行計算時，數值振蕩較為明顯，區域內壓力碎片狀分布相對無序。壓力波傳播過程勢必會受到大型虛假數值振蕩的干擾，而壓力信號真實分布特征受到擾亂。使用改進的DGEFM進行模擬時，壓力波整個時程清晰明確。從云圖可以明顯觀察到壓力波擴散傳播過程中反射和疊加等特征。

圖6 不同時刻聲壓分布情況(Newmark)

圖7 不同時刻聲壓分布情況(DGFEM)

圖8 層合液體網格劃分

圖9比較了Newmark方法和改進的DGFEM在t=0.003 s，t=0.006 s，t=0.01 s，t=0.013 s，t=0.016 s，t=0.02 s時的聲波傳播情況。可以看出，Newmark方法得到的結果隨著壓力波的傳播出現了數值振蕩，當壓力波到達層合界面處時，層合界面處表現出更強的數值振蕩，隨著壓力波的傳播，這種數值振蕩不斷增強，誤差累積增大。而改進的DGFEM可以很好濾除具有間斷特征的壓力波傳播過程以及層合介質間斷面兩側出現的數值振蕩，三維壓力波梯度結構明確，壓力波間斷特征較好描述。

4 結 論

本文構建了水下聲波傳播問題的改進時域間斷Galerkin有限元方法，并針對具有強間斷、高梯度變化特征的水下含多障礙物復雜邊界問題以及層合液體問題進行了求解。結果表明，通常的時域連續Galerkin有限元方法如Newmark方法求解這類問題時會在波前波后出現虛假的數值振蕩現象，隨著波的傳播，數值振蕩范圍會不斷擴大，誤差累積也不斷增大，影響問題計算精度。而改進的DGFEM能較好地消除聲波在不同時刻以及不同介質層合間斷面處波前波后的數值振蕩，使整個傳播過程清晰明確，具有較高的精度和可靠性。問題的求解，為進一步流固聲耦合問題研究奠定了基礎。

圖9 不同時刻聲壓分布情況