雙互易雜交邊界點方法求解Helmholtz方程＊

楊慶年，鄭俊杰，苗 雨

（1.華中科技大學土木工程與力學學院，湖北武漢 430074；2.南陽理工學院土木工程系，河南南陽 473004）

亥姆霍茲（Helmholtz）方程▽2u＋μ2u＝f是一個橢圓形偏微分方程［1－2］，在高波數的情況下它的解呈現高振蕩的特性，因此它的求解對數值方法是一個挑戰.傳統數值方法的求解精度有待提高，如有限元法的計算精度隨著方程波數的增大而急劇降低，這是因為有限元法普遍采用了低階的多項式作為位勢函數的插值函數，而低階多項式不可能對高頻震蕩的波傳播問題給出很好的近似，要提高計算精度，必須大幅度增加單元網格的數量，計算工作量大，效率低，不適應實際復雜大規模問題的求解.解決這個問題，單靠計算機硬件的發展是遠遠不夠的，必須尋找適應性較強的高效計算方法，以提高解決實際問題的能力.

雜交邊界點方法［3－6］是一種純邊界類型的無網格方法，該方法無論插值還是積分都不需要網格，計算時僅需要邊界上離散點的信息，前處理簡單，是一種很有發展前景的計算方法.然而，該算法同其他邊界類型方法一樣，在求解Helmholtz方程這樣的非齊次方程時需要域內劃分網格，從而失去了其降維的優勢.本文將雜交邊界點方法同雙互易法［7－8］結合，提出了一種求解Helmholtz方程的新方法，該方法將問題的解分為通解和特解兩部分，通解使用雜交邊界點法求解，特解則使用徑向基函數插值得到，數值算例表明該方法計算精度高，穩定性好，計算時需要域內布置節點，但僅用來徑向基函數插值，因此不影響該方法降維的優勢.

1 Helmholtz方程求解的雜交邊界點法

為方便起見，可假定f＝0，Helmholtz方程可表示為

式（1）的左端項利用雜交邊界點法求解拉普拉斯方程，右端項使用徑向基函數插值.根據雙互易定理，式（1）的解u可分為通解uh和特解up兩部分，即

特解up僅僅需要滿足非齊次方程

通解uh必須滿足拉普拉斯方程和修正邊界條件：

通解利用雜交邊界點法求解.雜交邊界點法基于修正變分原理，在修正變分原理中，假定域內位移場u，邊界位移場和邊界面力是相互獨立的，域Ω的邊界為Γ＝Γu＋Γt，和分別是Γu和Γt上的邊界值，相應的修正變分泛函ΠAB定義為

由δΠAB＝0，可得下面弱積分方程：式（8）和式（9）在域Ω的任何一部分都成立.例如，對于圖1所示的以sJ為圓心的子域Ωs［9］，邊界為Γs和Ls.對于子域Ωs，我們采用下面的弱積分形式代替式（8）和式（9）.

式中：h為檢驗函數.

圖1 子域Ωs以及相應基本解的源點SJFig.1 Local domainΩsand source point of fundamental solution corresponding to SJ

為了使Ls上積分為零，檢驗函數h取類似于移動最小二乘權函數的形式：

式中：dJ是域內任意點Q到節點sJ的距離.和用移動最小二乘近似［10］插值：

域內變量使用基本解近似：

基本解的形式為

式中：r是場點Q到源點PI之間的距離.

當u采用式（15）時，式（10）和（11）可寫成：

對所有節點列出上面方程，則有

式中：

2 雙互易雜交邊界點法

2.1 雙互易法

利用雙互易法求解Helmholtz問題是將非齊次應用出現的域內積分轉化為等價的邊界積分，對非齊次項進行插值，－μ2u用式（22）近似.

式中：αj是一系列初始未知參數；fj是近似函數；N和J分別是邊界點總數和域內點總數.

作為方程的相同形式，特解根據特解的基本形式進行插值：

如果up滿足方程（3），可得如下方程

選取fj＝1＋r＋r2，顯而易見，滿足式（24）的特解=u為

相應的面力為

通過式（22），（23）和（24），特解寫成矩陣形式：

2.2 雙互易雜交邊界點方法

雜交邊界點方法采用移動最小二乘插值形成形函數，不滿足δ特性，故邊界條件需作處理.本文采用實值虛值交換法處理修正的邊界條件［11］，有：

式中：RI，J＝［Φ，J（SI）］－1，N是該邊界節點總數；和是節點通解.

將式（27）～（30）代入式（2），然后將結果代入式（20）和（21），得

式（31）和（32）是雙互易雜交邊界點法求解亥姆霍茲問題的系統方程，假設N個點位于邊界上，通過式（31）和（32）可得到N個方程.但是，上面的方程包含L個內部點的位移，所以附加方程是必須的.

2.3 附加方程

利用式（31）和（32）不能求出域內點的變量，因此，本節研究求解Helmholtz問題的附加方程.

域內點未知變量用下式表達：

通解uc用基本解插值，特解up用式（27）表示，因此式（33）可表示為

式中：u＊是域內點的位移；uS是域內點基本解矩陣；是特解基本形式的值的矩陣.

由式（32）求得

將式（35）代入式（34）得

將式（35）代入式（31）得

聯立方程（36）和（37）得

其中

以上就是雙互易雜交邊界點法求解Helmholtz問題的理論推導過程.像邊界元法，本文提出的方法也存在“邊界層效應”，本文使用一種自適應積分方案來解決這一問題.雙互易雜交邊界點方法是一種僅在邊界上布點的無網格方法，無論插值還是積分都不需要邊界單元，域內的點僅僅是為了特解插值的需要，不影響該方法是一種邊界類型的方法.

3 數值算例

本文以兩種不同邊界條件的振動梁作為研究對象，為了驗證該算法的計算精度，相對誤差定義為：

式中：（e）和（n）分別表示解析解和數值解；u表示位移.

3.1 懸臂梁的振動

如圖2所示的懸臂梁，梁長0.9m，用40個邊界點和24個域內點離散，在固定端中部的邊界點（點40）處施加一個任意小的位移u0，其余所有點t＝0.

圖2 懸臂梁的計算模型Fig.2 Discretization of the cantilever

固有頻率的解析解為［8］

對于一維問題，變形模態為

式中：m是初始期望階次的整數值；l是梁的長度；C是任意常數.

根據式（39），梁的固有頻率的表達式為

式中：ρ和E是材料常數；當C被定義時，式（40）求出變形模態，本文取C＝1.

當利用式（38）計算時，取μ＝0，▽μ＝0.1.如果連續兩次迭代所有點的u值平均增長30%，Δμ可減為0.01；當u的平均值開始下降時，Δμ提高到0.5，直到u值再增加時為止.

計算時，把振動梁的邊界分為4片光滑的邊界段.在邊界上均勻布置40個點（在AD和BC上各布置19個點，在AB和CD上各布置1個點），即N＝40，在域內布置24個點，即L＝24.圖3為第1～4級固有頻率時梁的變形圖.從圖中可以看出，數值解和解析解吻合得非常好.

圖3 1～4級固有頻率懸臂梁的變形圖Fig.3 Deformed shape for first four natural frequencies of cantilever

為了分析計算精度對邊界點數量的敏感性，我們分為6種不同情況，每種情況的相對誤差如圖4所示.從圖中可以看出，域內點的數量對這類問題的影響還是比較明顯的，特別是高階振動模態.

圖4 1～4級固有頻率懸臂梁的變形圖Fig.4 Deformed shape for first four natural frequencies of cantilever

3.2 兩端固支梁的振動

如圖5所示的兩端固定梁，梁長0.9m，點的布置與圖2相同，在左端中部的邊界點（點40）處施加一個任意小的位移u0，其余所有點t＝0.

固有頻率的解析解為［8］

對于一維問題，變形模態為

表1為在不同離散方案情況下μ數值解與解析解的比較.

表1 μ數值解與解析解的比較Tab.1 Comparison of numerical and analytical solution forμ

圖5 兩端固定梁計算模型Fig.5 Discretization of the clamped－clamped beam

4 結 論

本文將雙互易法和雜交邊界點法結合求解Helmholtz方程，該算法無論插值還是積分都不需要網格，是一種邊界類型的純無網格法，域內布置少量節點僅僅用來徑向基函數插值.數值算例表明，該方法在求解Helmholtz方程時計算精度高，前處理簡單，適合求解此類問題.

通過數值算例可以得出以下結論：

1）域內點的位置和數量在低階振動模態分析中，對計算結果的影響并不大，但在高階振動模態的分析中是很重要的，由數值算例看，域內節點數量只要保證不低于邊界節點數量的一半就可以滿足計算精度要求.

2）子域半徑的大小對計算結果影響較大，一般取0.8h～0.9h較為合適，h為相鄰節點之間在參數空間的平均距離.

3）利用該算法求解振動梁的變形模態具有高精度和高收斂性，在雜交邊界點法求解中的邊界層效應通過自適應積分方案來解決.該方法還可用于求解更為復雜的非齊次問題.

［1］ CHEN J T，WONG F C.Dual formulation of multiple reciprocity method for the acoustic mode of a cavity with a thin partition［J］.Journal of Sound and Vibration，1998，217：75－95.

［2］ HARARI I，BARBONE P E，SLAVUTIN M，et al.Boundary infinite elements for the Helmholtz equation in exterior domains［J］.International Journal of Numerical Method in Engineering，1998，41：1105－1131.

［3］ ZHANG J M，YAO Z H，MASATAKA T.The meshless regular hybrid boundary node method for 2Dlinear elasticity［J］.Engineering Analysis with Boundary Elements，2003，27：259－268.

［4］ MIAO Y，WANG Y H.Development of hybrid boundary node method in two dimensional elasticity［J］.Engineering Analysis with Boundary Elements，2005，29：703－712.

［5］ ZHANG J M，YAO Z H，LI H.A hybrid boundary node method［J］.International Journal of Numerical Method in Engineering，2002，53：751－763.

［6］ ZHANG J M，TANAKA M，MATSUMOTO T.Meshless analysis of potential problems in three dimensions with the hybrid boundary node method［J］.International Journal of Numerical Method in Engineering，2004，59：1147－1168.

［7］ WROBEL L C，BREBBIA C A.The dual reciprocity boundary element formulation for non－linear diffusion problems［J］.Computer Methods in Applied Mechanics and Engineering，1987，65（2）：147－164.

［8］ PARTRIDGE P W，BREBBIA C A，WROBEL L C.The dual reciprocity boundary element method［M］.Southampton：Computational Mechanics Publications，1992.

［9］ ZHU T，ZHANG J D，ATLURI S N.A local boundary integral equation（LBIE）method in computational mechanics，and a meshless discretization approach［J］.Computational Mechanics，1998，21（3）：223－235.

［10］LANCASTER P，SALKAUSKAS K.Surfaces generated by moving least squares methods［J］.Mathematics of Computation，1981，37：141－158.

［11］ZHU T，ATLURI S N.A new meshless local Petrov－Galerkin（MLPG）approach in computational mechanics［J］.Computational Mechanics，1998，22（2）：117－127.