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

楊慶年，鄭俊杰，苗 雨

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

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

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

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

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

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

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

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

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

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

式中：h為檢驗函數(shù).

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

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

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

域內(nèi)變量使用基本解近似：

基本解的形式為

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

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

對所有節(jié)點列出上面方程，則有

式中：

2 雙互易雜交邊界點法

2.1 雙互易法

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

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

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

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

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

相應的面力為

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

2.2 雙互易雜交邊界點方法

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

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

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

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

2.3 附加方程

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

域內(nèi)點未知變量用下式表達：

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

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

由式（32）求得

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

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

聯(lián)立方程（36）和（37）得

其中

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

3 數(shù)值算例

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

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

3.1 懸臂梁的振動

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

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

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

對于一維問題，變形模態(tài)為

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

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

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

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

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

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

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

圖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］

對于一維問題，變形模態(tài)為

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

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

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

4 結 論

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

通過數(shù)值算例可以得出以下結論：

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

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

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

［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.