譜元方法求解波動方程及影響其數值精度的相關因素

摘 要：為探討波動方程的高精度數值模擬，采用Chebyshev譜元方法結合隱式Newmark時間積分方法求解波動方程，求解一個具體算例驗證了數值方法的可行性，討論了時間步長、Newmark因子以及計算區域的網格剖分方式對數值精度的影響，結果表明：和差分法相比，譜元方法求解波動方程具有所用網格節點少，數值精度高的特點；數值誤差隨時間步長減小而減小；在滿足穩定性要求的前提下，數值誤差隨著Newmark因子的減小而減小；當總網格節點數相同時，不同的網格剖分方式所得數值誤差不同，所述方法和結論可用于模擬聲波在空氣中的傳播。

關鍵詞：波動方程；譜元方法；時間積分方法；氣動聲學

中圖分類號：O353 文獻標志碼：A 文章編號：0253-987X(2008)01-0056-04

譜元方法是結合譜方法和有限元法各自特點的一種數值方法，它既具有譜方法的高精度和收斂特性，又具有較好的幾何區域適應性，20世紀80年代初由Patera提出，其后主要應用于計算流體力學和計算傳熱學，隨著對其研究的深入和工程實際問題的需要，該方法近來被應用于更多領域，如氣動聲學、地球物理學等領域。

計算氣動聲學是一門新興的交叉學科，主要通過數值計算研究氣動噪聲的產生機理和傳播特性、流體和聲的相互作用等，目前面，臨的主要問題之一是傳統數值方法的精度不能完全滿足計算氣動聲學的需要，如差分方法和有限元方法，聲學量是流動中脈動量相互抵消的結果，其值較平均流場量小3～5個量級，故需較高的數值精度，針對這一要求，譜元方法在這一領域的應用受到關注，譜元方法求解時間相關問題時，常牽涉求解形如下式的常微分方程組

Mu_t+Gu_t+Ku_t=Q_t (1)

4．1 時間步長h對計算精度的影響

為探討步長A對數值精度的影響，對式(18)分別采用人h=0.2，0.1，0.01，0.001s進行求解，網格剖分均為N_m=4，N_n=4，N_x=4，N_y=4，響應時間均為t=2s，圖1和圖2分別給出不同步長對應的各網格點上的絕對誤差，對比可得，誤差隨步長減小

4．2 δ因子對計算精度的影響

為探討δ對數值精度的影響，δ分別取0.6、0.7、0.8、0.9，對式(18)進行求解，網格剖分方式同4.1節，h=0.01s，響應時間均為t=2s，圖3和圖4分別給出不同δ值對應的各網格點上的絕對誤差，對比可得誤差隨s的減小而減小。

4．3 區域剖分方式對計算精度的影響

為探討當總節點數相同，但單元數和網格數均不同時，區域剖分方式對于計算精度的影響，分別采用不同剖分方式求解式(21)，剖分方式分別為：N_m=N_n=4，N_x=N_y=6(稱4×6剖分格式，反之亦然)；N_m=N_n=3，N_y=N_y=8(3×8剖分格式)；N_m=N_n=2，N_x=N_y=12(2×12剖分格式)，網格總數N_G=(N_mN_x+1)(N_nN_y+1)=625，h=0.01s，δ=0.5，t=2s。

定義ε＝(N_xN_y)/(N_mN_n)，當ε＞1時，表明單元內網格較密而單元較疏，當ε＜1時，表明單元數目較多而單元內網格較疏。

圖5～圖7給出ε＞1時的絕對誤差分布，比較可得，盡管以上幾種剖分方式對應的最大誤差幾乎相等，且產生在相同的總體節點處，但在其他節點上，ε越小則誤差越小。圖8～圖10給出ε＜1時的絕對誤差分布，比較可得當ε＜1時，計算區域采用不同的剖分格式時各節點的誤差均相差不大。

5 結 論

1)采用Chebyshev譜元方法結合Newmark方法求解波動方程是可行的，盡管Newmark方法牽涉到矩陣求逆，但計算結果仍然具有較高的數值精度。

(2)數值誤差隨著時間步長的減小而減小，受時間步長的影響很大，當h=0.2s時，最大絕對誤差e_a=0.272190257025；當h=0.001s時，最大絕對誤差e_a=0.000986182775，兩者相差約276倍。

(3)在滿足穩定性要求的前提下，數值誤差隨著δ的減小而減小，可以推斷，當δ＝0.5時，即中心加速度法，可得最小數值誤差。

(4)當總網格數相同而單元數不同且ε＞1時，最大數值誤差幾乎相同，但在其他網格節點上，ε越小則計算精度越高；當ε＜1時，計算精度幾乎不受ε影響，可以通過譜元方法的誤差分析得出其中原因。

(編輯 王煥雪)