多聯通封閉空間聲場響應的基于核重構的最小二乘無網格解法

李鴻秋，陳國平，史寶軍

（1.金陵科技學院 機電工程學院，南京 211169；2.南京航空航天大學 航空宇航學院，南京 210016；3.山東建筑大學 機電學院，濟南 250101）

研究封閉聲場響應問題的一個重要任務就是精確求解域內的亥姆霍茲方程。目前，此類問題的分析方法多為有限元法以及邊界元法，這兩種方法都是基于單元的方法，隨著頻率的增加，將聲域離散成大量單元，因此基于單元的方法主要集中在低頻段。近年來，Leuven研究小組［1］基于間接Trefftz方法推導了波疊加法（WB法），相比較于單元方法，WB法的精度和收斂速度都顯著提高，但是，WB法的收斂的充分條件是所求聲域為凸形［2］。對于復雜封閉聲腔特別是多聯通封閉聲腔響應的研究一直是難點問題。

無網格方法不需要網格劃分，只需要知道節點信息，因此在涉及網格畸變、網格移動等問題時具有明顯的優勢［3］?，F有的無網格方法主要有以下幾類：光滑粒子法［4］，核重構法［5-6］，散射元法［7］，無網格 Galerkin法［8］，Hp-Clouds 無網格法［9］，單位分解有限元法［10］，無網格局部 Petrov-Galerkin 法［11］以及有限點法［12］等。

本文將核重構無網格法應用于多聯通封閉空間聲場響應的研究中，應用最小二乘配點法離散控制微分方程，建立了亥姆霍茲方程問題的最小二乘配點格式。此方法的系數矩陣是對稱正定的，因此結果具有較好的穩定性。而本文采取的最小二乘配點格式對于邊界的處理方便直接。對于多聯通域問題，如果邊界形狀改變則只需改變邊界參數方程?；谂潼c分布的隨機性，可以根據邊界情況對節點數目進行合理有效的分配。提高計算的精確性。

1 多通域封閉空間聲場

假設有一封閉空間，結構分為彈性和剛性表面兩部分如圖1所示，考慮簡諧激勵情況忽略聲腔內的流體阻尼，則腔內聲壓分布p應滿足亥姆霍茲方程及封閉空間邊界條件：

圖1 亥姆霍茲方程問題的幾何表示Fig.1 The geometric representation of Helmholtz equation

2 核重構方法

核重構方法的關鍵就是通過建立一個修正核函數來構造函數u［13］，對于一維函數來說，帶有修正核函數的近似函數可表示為：

對于二維函數，帶有修正核函數的近似函數可表示為：

式中，C（x，s），C（x，y，s，t）稱為修正函數，通過構造適當的修正函數可以得到滿足不同精度要求的近似函數，本文中修正函數采用多項式基函數的線性組合表示，其最高階次數取決于控制微分方程的最高階導數，當C=1時，就是經典的光滑粒子法（SPH法），SPH法用配點格式建立離散方程組，比較容易實現，但是精度較低，穩定性也較差。

對于一維問題，修正函數的形式如下所示，：

二維問題：

式（5），式（6）中的c0（x）、c1（x）、c2（x）及c0（x，y）、c1（x，y）、…c5（x，y）表示未知的修正系數，如果控制微分方程只包含一階導數項則只需前三項即可滿足精度要求。修正系數可以通過設定近似函數的重構條件來確定。對于一維問題，將u（s）做泰勒展開：

將式（7），式（4）代入式（2），有：

其中，

并且：

于是，式（9）可以寫成如下的矩陣形式：

同理，對于二維問題，將u（s，t）做泰勒展開得：

將式（6）式（4）代入式（3）中得：

其中：且：

同樣的，可以得到重構條件如下：

因此，式（15）可以寫成如下的矩陣形式，且重構函數的系數可以由此式求出：

式（12），式（18）可縮寫為：

3 核函數系數導數求解

3.1 核函數系數一階導數的求解

近似函數ua（x，y）對x的一階導數為：

將式（4）代入上式，則：

上式可寫成如下的矩陣形式：

近似函數對y的一階導數為：

對y的一階導數的精確表達的重構條件為：

上式可寫成如下的矩陣形式：

3.2 核函數系數二階導數求解

對二維問題，式（23）進一步對x和y分別求導，式（26）進一步對y求導，同樣可以得到如下的矩陣形式：

3.3 基于核重構的配點法

核重構配點法是指將方程（1）的離散形式用于微分方程的求解，設在區域Ω中有N個離散點，則方程（1）的離散形式可寫做：

式中，xi表示節點I的坐標，于是，修正核函數的離散形式可表示為：

式（30）可寫成如下的形函數形式：

對于二維問題，方程（2）的離散形式為：

式（33）的形函數的形式寫成：

NI（x，y）的形式如下：

對于多維問題，其核函數通過一維問題核函數的乘積來求得，如對二維問題，有：

式中，αx，αy分別為沿x，y方向的精化參數，本文根據文獻［14］，取 αx=1.17Δx，αy=1.17Δy，Δx、Δy分別為x、y方向的節點間距。核函數取三次樣條函數，即：

于是形函數的導數可根據下式計算：

4 亥姆霍茲方程最小二乘配點格式

以方程（1）為例，說明最小二乘配點無網格方法。將近似函數ua（x，y）代入控制微分方程。最小二乘配點法的原理就是以控制微分方程在各離散點處誤差平方總和為最小的條件來建立求解試函數的方程。

設在作用域Ω內部有Ni個點，在邊界Γd上有Nd個點，在邊界Γn上有Nn個點，則共計布有Ni+Nd+Nn個點，以近似值代入式（1），則在Ni，Nd，Nn各離散點上的誤差分別為：

上式可寫做如下矩陣形式的殘量方程：

各離散點方差總和為最小的條件為：

于是可得到：

根據上式即可計算ui（i=1，2，…，NP），代入式（34）就可得到封閉空間聲場的響應值。

極坐標下Helmholtz方程的解為：

式中，Jn（kr）、Yn（kr）分別表示宗量為（kr）的n階柱貝塞爾函數和諾埃曼函數。如果多聯通域的內部空穴不是圓形時則一般不能求得解析解，為了與解析解進行對比判斷本文方法計算結果的精確性，特選取中心為圓孔的圓環形多聯通域，對此典型算例進行計算檢驗該方法的有效性。

對于多連通圓環域，內孔邊界上滿足：

環保高壓、農藥減量增效、行業競爭加劇，給農藥企業發展帶來機遇與挑戰。立足新形勢，企業如何實現持續發展，成為行業熱議的重要話題之一。

外圓邊界上滿足：

對于圓環狀多聯通域這樣的典型問題。式（44）對于所有的θ都成立，從而得到多聯通圓環區域的Helmholtz方程的通解為：

將式（45），式（46）代入式（44），根據邊界條件確定A0，B0即可。

為了對誤差進行分析和比較，本文采用如下的誤差范數：

根據上述原理編制了亥姆霍茲方程二維問題的最小二乘配點法的計算程序，并給出節點均勻分布以及隨機分布時的典型算例，與多聯通圓環區域的Helmholtz方程解析解進行比較，檢驗了該方法的精確性以及穩定性。

5 數值算例

設求解區域為一同心圓環，中心孔洞的半徑為0.2；外圓的半徑為1.0，示意圖如圖2所示。

例1 微分方程及邊界條件如下：

例2 微分方程及邊界條件如下：

對此問題采取兩種離散方案：一是離散點沿半徑均勻分布，二是離散點沿半徑隨機分布。

根據邊界條件求得例1的解析解為：

圖2 圓環作用域Ω示意圖Fig.2 A scheme ring of the scope Ω

例2的解析解為：

圖3和圖6分別給出了節點均勻分布與隨機分布時的分布情況。圖4和圖5分別是例1通過本文無網格數值方法以及解析方法得到的u值。圖7和圖8分別是例2通過本文無網格數值方法以及解析方法得到的u值。

圖3 24×24節點均勻分布圖Fig.3 A uniform point distributions of 24 ×24 points

圖4 節點均勻分布無網格方法的計算結果Fig.4 Numerical solutions with points uniform distribution for Example 1

圖5 例1解析解的計算結果Fig.5 Analytic solutions for Example 1

誤差情況：

表1 算例誤差Tab.1 Errors of the two examples

考慮到外周與內孔尺寸上的差別，因此在隨機分布節點時，讓外邊界上的點是內孔邊界節點數目的兩倍，其它各點在內域隨機分布?？梢钥闯?，當節點隨機分布以及均勻分布時所得出的結果都是穩定且精確的。

圖6 24×24節點隨機分布圖Fig.6 A random point distributions of 24 ×24 points

圖7 例2節點隨機分布無網格方法的計算結果Fig.7 Numerical solutions with points random distribution for Example 2

圖8 例2解析解的計算結果Fig.8 Analytic solutions for Example 2

表1是無網格方法計算結果與解析解相比較的誤差值，可以看出，無網格方法的計算精度還是比較高的。而且無論是節點均勻分布還是隨機分布所得到的結果都是比較精確且穩定的。本文給出的方法對節點分布情況并不敏感，且隨著節點數目的增加，計算精度也進一步增加。基于本文方法對節點的不敏感性，在外邊界或內邊界或者需要特別關注的區域設置節點更具靈活性。

6 結論

本文應用核重構法構造近似函數，并應用配點法以及最小二乘原理構造亥姆霍茲方程二維問題的最小二乘格式。相比較于Galerkin方法，配點法在求得區域內點的形函數及其導數后實施很直接。本文無網格方法給出的系數矩陣是對稱正定的，因此最終代數方程組的求解具有良好的穩定性。節點的分配無論是均勻還是隨機都具有較好的收斂性與精確性。

邊界條件問題一直以來都是無網格方法的難題，而本文采取的最小二乘配點格式的無網格方法對于邊界的處理是方便直接的。對于多聯通域問題，如果邊界形狀改變則只需改變邊界參數方程即可。此方法由于其對節點分布的不敏感，且可以在邊界上根據情況而調整節點數目，因此具有很大的優越性。通過數值算例可以發現，無論是節點均勻分布還是隨機分布的情況，計算精度都比較高，且具有較好的穩定性與收斂性。

［1］Pluymers B，Desmet W，Vandepitte D，et al，Application of an efficient wave based prediction technique for the analysis of vibro-acousstic radiation problems［J］.Journal of Computational and Applied Mathematics，2004，168：353-364.

［2］Desmet W.A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems［J］.Advances in Engineering Software，2002，33：527-540.

［3］Belytschko T，Krongauz Y，Organ D，et al.Meshless method：an overview and recent developments［J］.Compute Method in Applied Mechanics and Engineering，1996，139：3-47.

［4］Gingold R A，Moraghan J J.Smoothed particle hydrodynamics：theory and applications to non-spherical stars［J］.Mon Not Roy Astrou Soc，1977，181：375-389.

［5］Liu W K，Jun S，Zhang Y F.Reproducing kernel particle methods［J］.International Journal of Numerical Methods Fluids，1995，20：1081-1106.

［6］ Liu W K，Han W M，Lu H S，et al.Reproducing kernel element method part I：theoretical formulation［J］.Computer Methods in Applied Mechanics and Engineering，2004，193：933-951.

［7］Nayroles B，Touzot G，Villon P.Generalizing the finite element method：diffuse approximation and diffuse elements［J］.Computational Mechanics，1992，10：307-318.

［8］Belytschko T，Lu Y Y，Gu L.Element free galerkin methods［J］.International Journal of Numerical Methods in Engineering，1994，37：229-256.

［9］Duarte C A，Oden J T.Hp clouds：a H-p meshless method［J］.Numerical Methods for Partical Differential Equations，1996，12：673-705.

［10］ Melenk J M，Babuska I.The partition of unity finite element methods：basic theory and application［J］.Compute Method in Applied Mechanics and Engineering，1996，139：263-288.

［11］ Atluri S N，Zhu T.A new meshless local petrov-galerkin（MLPG）approach in computational mechanics［J］.Compute Method in Applied Mechanics and Engineering，1998，22：117-127.

［12］ Onate E，Idelsohn S，Zienkiewicz O C，et al.A finite point method in computational mechanics：applications to convective transport and fluid flow［J］.International Journal of Numerical Methods in Engineering，1996，39：3839-3866.

［13］ Aluru N R.A point collocation method based on reproducing kernek approximations［J］.International Journal for Numerical Methods in Engineering，2000，47：1083-1121.

［14］史寶軍，袁明武，李 君.基于核重構思想的最小二乘配點型無網格方法［J］.力學學報，2003，35（6）：697-706.