李鴻秋，姜金輝，陳國平，智淑亞

基于拓展單位分解有限元法分析聲波在多域場內的響應

李鴻秋1，姜金輝2，陳國平2，智淑亞1

(1. 金陵科技學院機電工程學院，江蘇南京 211169；2. 南京航空航天大學結構力學和控制國家重點實驗室，江蘇南京 210016)

基于拓展單位分解有限元方法，將平面波函數(shù)和貝塞爾函數(shù)作為基函數(shù)進行拓展。將亥姆霍茲方程離散，求解時不變情況下多域場內聲波的響應，并分析基函數(shù)對求解精度的影響。將波動方程的時間導數(shù)利用二階中心差分方法離散，得到方程的隱式表達式，劃分時間步迭代求解時變情況下聲波在多域場內的響應，分析迭代時間間隔對計算精度的影響，與典型算例的精確解進行比較，驗證精確性。結果表明，平面波函數(shù)作為拓展基函數(shù)，利用二階中心差分法離散時間導數(shù)，分析時不變以及時變情況下多域場內高波數(shù)聲波的響應問題，具有較高的計算精度和計

波動方程；亥姆霍茲方程；單位分解有限元；拓展方法；時變和時不變

0 引言

研究聲波的問題歸根結底是求解不同邊界條件下的波動方程，數(shù)值解法大多基于有限元方法[1-2]、有限差分法[3-4]和邊界元方法[5-6]。當聲波的響應問題涉及到復雜區(qū)域，多域場以及高波數(shù)等情況時，上述方法就會面臨計算量過大或者求解效率較低的問題[7]。尤其對于高波數(shù)問題，BUBASKA等[8]的研究表明，即使布置足夠多的節(jié)點，由于存在污染誤差(pollution error)，精度也很難得到保證。近年來，研究者們又發(fā)展了譜有限元、邊界元方法以及無網格方法計算聲波的傳播和響應問題[9-11]，但是對于高波數(shù)和多域或者復雜域等情況，計算效率和計算精度依然不甚理想。

單位分解有限元方法(Partition of Unity Finite Element Method, PUFEM)能夠在現(xiàn)有網格基礎上通過附加函數(shù)提高計算效率。LAGHROUCHE等[12]的研究表明，PUFEM求解高波數(shù)問題時不必細化有限元網格，SHAID等[13]利用高斯函數(shù)作為拓展基函數(shù)探討了熱場中熱的傳播問題，并將此方法應用到瞬態(tài)熱的求解即時域問題。LUOSTARI等[14]應用極弱變分方法，TEZAUR等[15]應用非連續(xù)拓展方法以及MONK等[16]利用最小二乘法將單位分解方法應用于邊界元方法求解聲波和彈性波問題。MOHAMED等[17]分析了上述幾種方法的求解局限性和計算效率問題，同時利用平面波作為拓展基函數(shù)用于求解不同邊界條件下的波動方程，包括多個散射源的問題，取得了不錯的計算精度和計算效率。同時，也有其他拓展有限單元方法[18]用于求解波的傳播問題，但這些研究往往局限于時間諧波，分析時不變情況下的聲波的傳播和響應問題。

本文利用四邊形等參單元的形函數(shù)，采用平面波函數(shù)或者貝塞爾函數(shù)作為基函數(shù)，經拓展建立單位分解有限元的整體形函數(shù)。構造求解亥姆霍茲方程整體剛度矩陣。將波動方程的時間導數(shù)利用二階中心差分方法離散，得到波動方程的隱式表達式，求解時變情況下聲波在多域場內的響應。

1 聲波的響應

1.1 波動方程

聲波的傳播過程可以描述成波動方程：

1.2 亥姆霍茲方程

波的傳播問題都可以歸結為在給定的邊界條件和初始條件下求波動方程的解。用分離變量的方法求解波動方程，則可以得到不考慮時間問題的亥姆霍茲方程(2)，以及邊界條件式(3)。

2 拓展單位分解有限元法

拓展單位有限元法的核心思想是由單位分解函數(shù)和基函數(shù)共同逼近求解空間，基函數(shù)不依賴有限元網格，不增加網格，可以通過拓展基函數(shù)的階次提高計算效率和計算精度[19]?；瘮?shù)的拓展方式可以是均勻的，也可以是不均勻的，不同的基函數(shù)有不同的計算精度[20]。

2.1 拓展基函數(shù)

從而得到：

將基函數(shù)和有限元形函數(shù)的乘積看作一個整體，形成新的形函數(shù)即：

2.2 誤差分析

3 時不變情況下聲波在多域場的響應分析

3.1 構造給定邊界條件亥姆霍茲方程的積分弱形式

(a) 平面波入射內含圓形邊界(Robin邊界或Neummn邊界)的求解區(qū)域

(b) 平面波入射經圓形邊界(Neummn邊界)反射后的求解區(qū)域 圖1 求解區(qū)域及求解邊界示意圖 Fig.1 Schematic diagram of the solving models

平面聲波入射，，用權重系數(shù)分別去乘亥姆霍茲方程式(2)和邊界條件式(3)的左右兩端，得到：

(10)

(11)

對式(10)應用格林公式得到：

(12)

將式(11)代入式(12)，得到積分弱形式(13)：

(13)

3.2 構造拓展單位分解有限元的整體形函數(shù)

采用四邊形等參元作為單位分解單元，劃分單元(表示單元序號)，求解域為，用平面波函數(shù)作為基函數(shù)在單位分解單元的節(jié)點處進行拓展，拓展個數(shù)為個，均勻分布，因此，，其中，。則拓展單位分解有限元的整體形函數(shù)可以寫成：

(14)

局部節(jié)點響應值為： ，則單元內其它點響應值可用式(15)求得：

(15)

3.3 構造時不變情況下拓展單位分解有限元法的剛度矩陣和求解式

設有轉換矩陣，使，將各單元自由度按對號入座的方式擴展到整體自由度矩陣，可得下式：

(16)

將式(15)和(16)代入亥姆霍茲方程積分弱形式(13)，得到矩陣形式(17)：

(17)

顯然，剛度矩陣為。利用式(17)求。

3.4 數(shù)值算例

算例1

如圖1(a)所示，入射平面聲波為，入射角度為，振幅為單位幅值，求解域為半徑為1的圓形區(qū)域，其中心位于(1, 1)，域內包含兩個半徑分別為0.5、圓心位置(0.5, 1.2)和半徑為0.25、圓心位置(1.5, 0.5)的圓。入射角度分別為α=0°、30°、60°、90°，波數(shù)分別為k=5π、8π、15π、20π。模型的單元個數(shù)為60。利用高階Gauss-Legendre積分方法對式(15)積分。圖2~5給出了其中4種情況下的聲響應值，表1給出了部分計算結果的誤差。

算例2

如圖1(a)所示，入射平面波經由滿足Neumann邊界條件(8b)的圓形表面反射后的響應問題，外截斷邊界滿足Robin邊界條件式(8a)。本例中[21]，，其中r、θ分別表示所在位置的極坐標，、分別表示n階漢克爾函數(shù)和貝塞爾函數(shù)，、是其對應的微分形式。此算例精確解包含貝塞爾函數(shù)，可以采用貝塞爾函數(shù)作為基函數(shù)，即令。算例2 分別采用平面波函數(shù)和貝塞爾函數(shù)作為基函數(shù)，拓展個數(shù)均為8個，計算當入射角度分別為α=0°、60°、90°，波數(shù)分別為k=5π、15π、20π，振幅為單位幅值時的聲響應值。結果表明，使用平面波函數(shù)作為基函數(shù)拓展，計算結果優(yōu)于使用貝塞爾函數(shù)作為基函數(shù)拓展的情況(見表2)，(a)為基函數(shù)為平面波函數(shù)時的誤差值，(b)為基函數(shù)為貝塞爾函數(shù)時的誤差值。圖6~9所示為采用平面波函數(shù)為基函數(shù)求解得到的聲波響應數(shù)值。

圖2 聲波響應勢能值的實部(k=8π, α=0°) Fig.2 The real part of the potential energy of acoustic wave response(k=8π, α=0°)

圖3 聲波響應勢能值的實部(k=8π, α=60°) Fig.3 The real part of the potential energy of acoustic wave response(k=8π, α=60°)

圖4 聲波響應勢能值的實部(k=15π, α=90°) Fig.4 The real part of the potential energy of acoustic wave response(k=15π, α=90°)

圖5 聲波響應勢能值的實部(k=20π, α=0°) Fig.5 The real part of the potential energy of acoustic wave response(k=20π, α=0°)

表1 不同波數(shù)和入射角情況下的相對誤差(算例1) Table 1 The relative errors for different wave numbers and different incident angles (Example 1) kα/(°)εL2 /% 5π02.2570×10-4 5π604.1592×10-4 8π60 2.40×10-2 10π90 1.89×10-2 15π01.5956×10-5 15π901.1530×10-5

表1和表2的誤差結果表明，雖然單元網格劃分較少，即使單元尺寸已經遠遠超過了波長，依然能獲得較高的計算精度。且當內邊界為Neumann邊界條件時，精確解包含貝塞爾函數(shù)情況下，采用平面波函數(shù)作為基函數(shù)拓展，同樣具有較高的計算精度。

圖6 聲波響應勢能值的實部(k=5π, α=0°) Fig.6 The real part of the potential energy of acoustic wave response(k=5π, α=0°)

圖7 聲波響應勢能值的實部(k=5π, α=60°) Fig.7 The real part of the potential energy of acoustic wave response(k=5π, α=60°)

圖8 聲波響應勢能值的實部(k=10π, α=60°) Fig.8 The real part of the potential energy of acoustic wave response(k=10π, α=60°)

圖9 聲波響應勢能值的實部(k=10π, α=90°) Fig.9 The real part of the potential energy of acoustic wave response(k=10π, α=90°)

表2 不同波數(shù)和入射角情況下的相對誤差(算例2) Table 2 The relative errors for different wave numbers and different incident angles (Example 2) kα/(°)εL2(a)/%εL2(b)/% 5π02.2570×10-22.67 5π604.1592×10-24.34 10π901.896.54 10π602.715.74

4 時變情況下聲波在多域場的響應分析

4.1 構造給定邊界條件亥姆霍茲方程的積分弱形式

設聲波在空間傳播的時間為，域內的波動方程和邊界條件可以寫成式(1)的形式。則其初始狀態(tài)為

(18)

時間步長為，表示在時間步時的聲響應值。使用隱式方法可以避免顯式方法中需要的小時間步長，本文采用二階中心差分法離散時間導數(shù)，迭代公式寫成式(19)的形式：

(19)

同樣，是在時間步時的響應值。顯然，受時間步長和網格劃分的影響。同樣用權函數(shù)w去乘式(19)的兩端，并在空間和邊界上積分，利用散度定理和邊界條件(8a)，得到此類邊界條件下波動方程的積分弱形式：

貴州是我國紅粘土的主要分布區(qū)，紅粘土厚度與中國北部的黃土相比厚度較薄，一般厚度為8-10米，局部地區(qū)厚度可以達到20多米。厚度的變化主要與下部基巖起伏有關，在巖溶緩坡、盆地等地厚度較大，峰林高地、陡坡等地厚度較薄甚至沒有。在紅粘土堆積厚的地方有著明顯的垂直分帶結構，從上到下可分為表土層、全風化層、半風化層和基巖四個結構，而在紅粘土堆積較薄的地方垂直分帶結構不明顯。為了更加詳細準確的反映紅粘土剖面特征，本次實驗選取的8個剖面都發(fā)育良好，其中包括貴陽白云巖紅粘土剖面兩組、石灰?guī)r紅粘土剖面一組，安順白云巖紅粘土剖面一組、石灰?guī)r紅粘土剖面兩組，遵義白云巖紅粘土剖面、石灰?guī)r紅粘土剖面各一組。

4.2 構造時變情況下拓展單位分解有限元法的剛度矩陣和求解式

在求解域劃分單元，求解域表示為，采用平面波函數(shù)作為基函數(shù)，拓展個數(shù)為。同樣將式(16)、(17)代入波動方程的積分弱形式(20)，可以得到矩陣求解式：

(21)

式中，。顯然，如果能夠得到在時間步、時的聲響應值和，則可以獲得在時間步時的聲響應值。

我做夢也想不到李小樹會千里迢迢去尋找許春花，可是他真的走了，走的時候，除了給大黑貓帶來一些他家里剩余下的貓食外，他還向我討要許春花的那幅肖像?？吹剿麩崆衅笈蔚难凵瘢冶惆旬嫺逅徒o了他。李小樹小心翼翼地把畫稿收藏在他事先準備好的一個畫筒里，就火急火燎地走了，像陣風似的消失在我的視線里。

4.3 數(shù)值算例

算例3

如圖1(a)所示，平面波入射，本算例的求解區(qū)域為凹字形狀的不規(guī)則區(qū)域，初始條件為：。時變問題解包含時間，其精確解為，分別計算k=5π、10π，α=0°、60°情況下在時間的聲響應值。圖10、11分別給出了在T時刻的部分結果，相對誤差見表3。

圖10 聲波響應勢能值的實部(k=5π, α=0°,T=0.01 s) Fig.10 The real part of the potential energy of acoustic wave response(k=5π, α=0°,T=0.01s)

圖11 聲波響應勢能值的實部(k=5π, α=60°, T=0.005 s) Fig.11 The real part of the potential energy of acoustic wave response(k=5π, α=60°, T=0.005 s)

表3 不同波數(shù),入射角和不同時間步情況下的相對誤差(算例3) Table 3 The relative errors for different wavenumbers, different incident angles and different time steps (Example 3) kα/(°)T/sεL2 /% 5π00.0010.001 5 5π00.0020.000 2 10π600.0010.005 6 10π600.0050.004 7

算例4

平面波入射，被圓心為(0, 0)、半徑為1的圓形剛性壁反射后的響應，令求解區(qū)域是圓心為(4,4)、半徑分別為的多域空間。初始條件為：。精確解包含時間項：。

精準醫(yī)療是針對于患者醫(yī)療保健和健康的個性化醫(yī)學模式，它通過醫(yī)生的醫(yī)療決策和實踐制定出適合不同疾病人群的治療方案。隨著對CRSwNP的發(fā)病機制的不斷深入了解，精準醫(yī)療分析整合疾病的診斷和治療并能制定出最優(yōu)化的治療方案[28]。而實現(xiàn)精準醫(yī)療的基礎必須具備的要素有：患者參與治療方案的決定；預判初始治療的成功率；防治疾病進展的有效策略和疾病內在型為驅動的個性化治療[29]。為了實現(xiàn)疾病內在型為驅動的治療目的，必須對疾病的內在型有著充分且標準化的認識，而且能夠洞察用于評估或預測療效、指導完善臨床策略的生物標記物[10]。

分別計算k=5π、6π，α=0°、60°情況下的時變響應值，部分結果圖見圖12和圖13所示，相對誤差見表4。算例3和算例4的結果顯示，利用二階差分離散時間導數(shù)構造整體計算矩陣，基于拓展單位分解有限元法計算時變情況下聲波在多域場的響應問題，計算結果精度較高，且當時間步長滿足小于1.0×10-3的情況下，如果繼續(xù)降低時間步長，計算精度變化不大。

圖12 聲波響應勢能值的實部(k=5π, α=0°, T=0.01 s) Fig.12 The real part of the potential energy of acoustic wave response(k=5π, α=0°, T=0.01 s)

圖13 聲波響應勢能值的實部(k=6π, α=60°,T=0.03 s) Fig.13 The real part of the potential energy of acoustic wave response(k=6π, α=60°,T=0.03 s)

表4 不同波數(shù),入射角和不同時間步情況下的相對誤差(算例4) Table 4 The relative errors for different wavenumbers, different incident angles and different time steps (Example 4) kα/(°)T/sεL2 /% 5π00.000 10.050 3 5π00.000 20.048 2 6π600.000 10.136 2 6π600.000 30.0125 1

5 結論

基于拓展單位分解有限元方法求解高波數(shù)聲波在多域場內的響應，分別利用平面波函數(shù)，貝塞爾函數(shù)作為基函數(shù)進行拓展，分析時不變情況下聲波在多域場內的響應問題。利用二階中心差分法離散時間導數(shù)，同樣基于單位分解有限元法構造波動方程的積分弱形式，用于求解時變情況下聲波在多域場內的響應問題，計算結果表明：本文構造的拓展基函數(shù)用于波動方程和亥姆霍茲方程的求解計算效率和計算精度較高，并且適用于不規(guī)則、多域、復雜區(qū)域以及高波數(shù)情況下聲波響應的求解，此方法對網格的劃分數(shù)量要求較低，即使網格數(shù)量較少，甚至波長大于單元長度依然可以取得較高的計算精度，不增加單位個數(shù)通過增加拓展基個數(shù)也可以提高計算精度。

參考文獻

[1] KOHNO H, BATHE K J, WRIGHT J C. A finite element procedure for multiscale wave equations with application to plasma waves[J]. Computers & Structures, 2010, 88(1-2): 87-94.

[2] HARARI I. A survey of finite element methods for time-harmonic acoustics[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(13): 1594-1607.

[3] PERREY-DEBAINE E, TREVELYAN J, BETTESS P. Wave boundary elements: a theoretical overview presenting applications in scattering of short waves[J]. Engineering Analysis with Boundary Element, 2004, 28(2): 131-141.

[4] HUTTUNEN T, MONK P. The use of plane waves to approximate wave propagation in anisotropic media[J]. Mathematics of Computation, 2007, 204: 350-367.

[5] BUFFA A, HIPTMAIR R, SCHWAB C. Boundary element methods for Maxwell transmission problems in Lipschitz domains[J]. Numerical Mathematics, 2003, 95(3): 459-485.

[6] PERREY-DEBAIN E, TREVELYAN J, BETTESS P. Wave boundary elements: a theoretical overview presenting applications in scattering of short waves[J]. Engineering Analysis with Boundary Elements, 2004, 28(2): 131-41.

[7] GUDDATI M N, YUE B. Modified integration rules for reducing dispersion error in finite element methods[J]. Computational Methods in Applied Mathematics, 2004, 193(3-5): 275-287.

[8] BABUSKA I M, SAUTER S A. Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers[J] SIAM REVIEW, 2000, 42(3): 451-484.

[9] CHAKRABORTY A, GOPALAKRISHNAN S. A spectral finite element model for wave propagation analysis in laminated composite plate[J]. Journal of Vibration and Acoustics, 2006, 128(4): 477-488.

[10] HAM S, LAI B, BATHE K J. The method of finite spheres for wave propagation problems[J]. Computers and Structures, 2014, 142: 1-14.

[11] GENECHTEN B, VANDEPITTE D DESMET W. A trefftz-bassed numerical modelling framework for Helmholtz problems with complex multiple scatter configurations[C]//Proceedings of the ICCES Special Symposium on Meshless & Other Novel Computational Methods; 2008, Suzhou, China.

[12] LAGHROUCHE O, BETTESS P, PERREY-DEBAIN E, et al. Wave interpolation finite elements for Helmholtz problems with jumps in the wave speed[J]. Computer Methods in Applied Mechanics and Engineering, 2005, 194(2): 367-381.

[13] MOHAMED S, MOHAMMED S, TREVELYAN J. A partition of unity FEM for time-dependent diffusion problems using multiple enrichment functions[J]. International Journal for Numberical Mehtods in Engineering, 2013, 93(3): 245-265.

[14] LUOSTARI T, HUTTUNEN T, MONK P. The ultra weak variational formulation using Bessel basis functions[J]. Communications in Computational Physics, 2012, 11(2): 400-414.

[15] TEZAUR R, KALASHNIKOVA I, FARHAT C. The discontinuous enrichment method for medium-frequency Helmholtz problems with a spatially variable wavenumber[J]. Computer Methods in Applied Mechanics and Engineering, 2014, 268(1): 126-140.

[16] MONK P, WANG D Q. A least-squares method for the Helmholtz equation[J]. Communications in Numerical Methods in Engineering, 1999, 175(1): 121-136.

[17] MOHAMED M S, MOHAMMED S, LAGHROUCHE O. An enriched finite element model with q-refinement for radiative boundary layers in glass cooling[J]. Journal of Computational Physics, 2014, 258: 718-737.

[18] VERGOT K, VANMAELE C, VANDEPITTE D. An efficient wave based approach for the tiem-harmonic vibration analysis of 3d plate assembly[J]. Journal of Sound and Virbration, 2013, 332(8): 1930-1946.

[19] CHRISTODOULOU K, LAGHROUCHE O, MOHAMED M S. High-order finite element for the solution of Helmholtz problems[J]. Computers and Structures, 2017, 191: 129-139.

[20] MOHAMED M S, LAGHROUCHE O, EL-KACIMI A. Some numerical aspects of the PUFEM for efficient solution of 2D Helmholtz problems[J]. Computers and Structures, 2010, 88(2): 1484-1491.

[21] LINTON C M, EVANS D V. The interaction of waves with arrays of vertical circular cylinders[J]. Journal of Fluid Mechanics, 1990, 215(1): 549-569.

Analysis of acoustic wave response in complex multiple domains based on the enriched partition of unit finite element method

LI Hong-qiu1, JIANG Jin-hui2, CHEN Guo-ping2, ZHI Shu-ya1

(1. College of Mechanical and Electrical Engineering, Jinling Institute of Technology, Nanjing 211169, Jiangsu, China; 2. State Key Laboratory ofMechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, Jiangsu, China)

Abstract: Based on the enriched partition of unity finite element method, the plane wave functions or Bessel functions are enriched as basis functions to discretize Helmholtz equation, the response of acoustic wave in complex multiple domains under time-independent condition is solved and the effect of basis functions on solution accuracy is analyzed. The two-order central difference method is used to discrete time derivative of the wave equation for implicit solution and the time step iterative method is used to solve the response of acoustic wave in complex multiple domains under time-dependent condition. By taking the plane wave function enriched as basis function and using the two-order central difference method to discrete time derivative, the responses of high wavenumber acoustic waves in complex multiple domains under time-dependent and time-independent conditions are calculated. Also, the effect of iterative time interval on calculation accuracy is analyzed. Typical examples are used to compare the calculation results with the exact solutions and to verify the method in this paper having high accuracy and efficiency.

Key words: wave equation; Helmholtz equation; partition of unity finite element;enriched method; time-dependent and time-independent

算效率。

中圖分類號：O429

文獻標識碼：A

文章編號：1000-3630(2019)-05-0481-07

DOI編碼：10.16300/j.cnki.1000-3630.2019.05.001

收稿日期: 2018-05-03;

修回日期: 2018-06-18

基金項目: 江蘇省自然科學基金資助(BK20151099)；金陵科技學院科研基金資助(JIT-B-201219)；江蘇省高校優(yōu)秀中青年教師和校長境外研修計劃資助

作者簡介:李鴻秋(1973－), 女, 河北邯鄲人, 博士研究生, 副教授, 研究方向為振動與噪聲控制。

通訊作者: 李鴻秋,E-mail: lihongqiu@jit.edu.cn