結構聲學耦合隨機性分析的等幾何有限元-邊界元法研究

胡昊文， 王中王， 徐延明， 陳磊磊

(1.信陽師范學院 建筑與土木工程學院，河南 信陽 4640002.黃淮學院 建筑工程學院，河南 駐馬店 463000)

在實際工程問題中存在著諸多不確定因素，如制造誤差、材料參數不確定、隨機載荷等，這些不確定因素會對結構響應產生影響。近幾十年來，隨著概率論的發(fā)展，很多基于概率論的不確定性分析方法已被提出和應用到工程分析中。蒙特卡羅模擬(Monte Carlo simulation,MCs)[1-4]是求解不確定性問題中的最具代表性的方法之一，其主要思想是通過大量采樣確定系統(tǒng)輸入和輸出之間的統(tǒng)計學關系。該方法操作簡單，因而得到了廣泛應用。但MCs需要大量的樣本來保證其準確性，這導致其計算成本高，難以進行大規(guī)模實際問題分析。Ding等[5]采用本征正交分解(proper orthogonal decomposition，POD)對基于蒙特卡羅模擬的多維不確定性等幾何隨機方法進行加速計算；Doucet等[6]對蒙特卡羅模擬的采樣方法進行改進，提出了自適應抽樣方法，在一定程度上減少了取樣個數，降低了計算量。本文通過POD構造代理模型，利用徑向基函數(radial basis function，RBF)近似求解物理場，從而高效地獲得大量樣本點，并用于基于蒙特卡羅模擬的隨機性分析。

水下結構的輻射與散射聲場分析是一項非常重要的研究課題，在多個工程應用領域中都很受重視。在船舶或潛艇這樣大型復雜結構的振動聲學預報中，由于結構的機械阻抗只比水的聲阻抗大數倍至數十倍，所以水對結構振動的影響不可忽略。目前在結構振動分析領域中，有限元方法是使用最廣泛的，而針對聲場問題，特別是外聲場問題，邊界元法[7-8]因其獨特的優(yōu)勢倍受學者關注。因此，將有限元與邊界元相結合構成耦合有限元與邊界元法(finite element method and boundary element method,FEM-BEM)進行聲振耦合分析兼具了兩種方法的優(yōu)點[9]。Zheng等[10-11]將FEM-BEM和圍道積分方法相結合發(fā)展了一套適用于耦合系統(tǒng)特征值分析的計算方法；Chen等[12]將低頻快速多極算法和對角式快速多極算法聯合起來，形成寬頻快速多極算法，并將其應用到FEM-BEM中，進行水下結構振動快速響應分析。這些研究大多采用傳統(tǒng)拉格朗日函數進行幾何形狀近似與物理場插值計算，雖然操作簡單，但模型誤差會導致計算精度的降低；此外，由于需重復進行CAD(Computer-Aided Design)與CAE(Computer-Aided Engineering)模型數據交互，因此計算效率低下。

本文將細分曲面法應用于耦合有限元-邊界元法中進行水下復雜結構振動聲學響應分析。采用MCs進行參數不確定性分析，并采用POD與RBF加速基于MCs的隨機性分析計算。

1 蒙特卡羅模擬

蒙特卡羅模擬的本質是通過抽取一定量的輸入參數樣本，并計算在這些樣本點處的響應值，最終獲得響應函數的統(tǒng)計特征(例如數學期望、標準差)。這些統(tǒng)計特征，本質上是一個積分問題。我們先給出蒙特卡羅積分

(1)

(2)

(3)

針對N個輸入變量的隨機樣本xi(i=1,2,…,N)，將狀態(tài)函數在這些樣本點處的值p(xi)的算術平均值作為Q的估計值，即

(4)

根據式(4)可知，蒙特卡羅模擬的預測精度取決于樣本數，樣本點數量越大則預測精度越高，但計算量也會越大。這一特點限制了蒙特卡羅方法的適用性。為了克服這個困難，本文組合本征正交分解和徑向基函數，對基于MCs的隨機性分析進行加速求解。

2 降階計算

使用部分樣本點的響應，構造一個解空間Λ，然后將解空間Λ分解成一組正交基底和增補基底，為構建代理模型提供子空間。引入m個輸入變量αs及其響應λ(αs)，構造一個解空間Λ，如下

(5)

其中：Λ∈Rn×m；n為單個輸入變量下的響應個數；m為變量個數。使用本征正交分解中的奇異值分解法(singular value decomposition,SVD)可將矩陣Λ分解成如下式子

(6)

式中：r=min(m,n)；U∈Rn×n，V∈Rm×m為正交矩陣；uij，vij分別為這兩個矩陣里對應的元素；uj為ΛΛT的特征向量，也稱作Λ的左奇異特征向量；vj為ΛTΛ的特征向量，也叫做Λ的右奇異特征向量；Σ∈Rn×m是一個對角矩陣，對角元素σj按從大到小排列。通過定義γj=uj，Aj(αi)=σjvij，可得到如下表達式

(7)

式中：γj被稱為正交基；Aj(αi)為相應的基底。使用式(7)，系統(tǒng)響應可以通過簡化表達的γj和Aj(αi)的線性組合來表示，且其規(guī)模比初始模型要小得多。遵循Li等[26]和Wang等[27]的思想，我們使用徑向基函數對原始輸入變量響應函數進行插值，可獲得任意參數下的系統(tǒng)響應，如下

(8)

式中：φ為基函數；η為其系數。本文選取常用的高斯核函數作為基函數，如下

(9)

式中：記號‖·‖為歐幾里得范數；系數γi由基函數的寬度決定[28]。

使用式(8)我們可以得到下列線性方程組

(10)

通過求解方程組，我們可以得到權系數η。因此，(7)式可改寫為

(11)

基于上述分析，任何隨機變量下的系統(tǒng)響應都可以通過這種線性組合來近似而不必使用完整的數值計算。根據這個特點，我們可以使用很小一部分取樣數據構成的矩陣來近似全階系統(tǒng)響應，從而避免蒙特卡羅模擬對全階系統(tǒng)的使用，也就避免了大量樣本點響應的重復計算。

3 細分曲面法

細分規(guī)則分為拓撲規(guī)則和幾何規(guī)則。拓撲規(guī)則描述新頂點與其單元和相鄰頂點之間的關系，幾何規(guī)則描述新頂點位置的計算。

3.1 Loop細分模式的細分規(guī)則分析

(12)

(13)

圖1 頂點V和邊緣節(jié)點E分布Fig.1 Distribution of edge point E and vertex point V

圖2 Loop細分幾何規(guī)則Fig.2 Loop subdivision refinement

3.2 曲面擬合分析

經過Loop細分之后，須對細分網格進行擬合操作。以規(guī)則單元為例，介紹曲面擬合操作過程。首先需要對該單元(標記為e)的鄰近單元進行重新編號，如圖3所示?？梢园l(fā)現，單元e的鄰近單元有12個，因此具有12個頂點，與這些頂點相對應的12個基函數用于進行曲面擬合插值計算，曲面擬合技術如下式

(14)

式中：xe為局部參數坐標系中點(θ1，θ2)的曲面坐標；Bi為箱樣條基函數[29]；Ci為12個控制頂點的坐標。圖4顯示參數空間和幾何空間的轉換關系。針對不規(guī)則單元的擬合方法見陳磊磊等[30]的工作。最終由Loop細分方法構造出的曲面在規(guī)則頂點處具有C2連續(xù)性，在不規(guī)則頂點上具有C1連續(xù)性。

圖3 有12個控制頂點的規(guī)則三角形單元Fig.3 A regular triangle with 12 control vertices

圖4 用3個規(guī)則頂點擬合三角形單元Fig.4 Fitting of the triangle element with three regular vertices

4 基于細分曲面的結構聲學耦合分析

4.1 結構有限元分析

水下薄殼結構振動響應模型如圖5所示，考慮外部激勵為簡諧波荷載的情況，采用薄殼平衡方程的弱形式建立有限元方程，如下表示

(K-iωC-ω2M)u=Kdu=f

(15)

式中：K，C和M分別為系統(tǒng)剛度矩陣、阻尼矩陣、質量矩陣，均由單元矩陣集合而成；u為節(jié)點位移向量；f為節(jié)點荷載向量。

圖5 水下結構聲學耦合模型Fig.5 Structural-acoustic interaction model

4.2 聲學邊界元分析

對于外聲場問題，控制方程為Helmholtz方程，對應的邊界積分方程表達式如下

(16)

式中：x為源點；y為場點；p(x)為聲壓；pi(x)為入射波聲壓；q=?p/?n，n為邊界外法線方向；G(x,y)是格林函數。系數c(x)由點x處的幾何特征決定，若點x處光滑，c(x)=1/2(下述推導均假設x處光滑)。式(16)被稱為常規(guī)邊界元方程(CBIE)，將其對外法向n(x)求導，可得法向導數邊界積分方程(NDBIE)，表達式如下

(17)

僅使用式(16)或式(17)求解外聲場問題時，在某些特殊頻率(虛假頻率)處會出現計算結果不準確現象，即解的非唯一性問題。通過組合CBIE與NDBIE構成Burton-Miller表達式，可在這些虛假頻率處獲得準確結果。

對式(16)和(17)進行線性組合，并采用細分曲面網格對結構邊界進行離散，可以得到如下線性方程組

Hp=Gq+pin

(18)

式中：H和G為基于細分曲面法的邊界元系數矩陣；p和q為邊界上的聲壓和聲通量向量；pin為入射波對應的節(jié)點向量。

4.3 基于細分曲面法有限元邊界元耦合分析

考慮聲場對結構振動的影響，將式(15)中的薄殼結構外部節(jié)點荷載向量f分解為激勵荷載fs和聲載荷fp。聲荷載表達式如下

(19)

(20)

式中，e1，e2，e3為3個正交的笛卡爾基向量，那么全局荷載向量表示為

f=Csfp+fs

(21)

流體和結構的法向速度分別表示為vn,f和vn,s。如果結構域和流體域之間沒有能量損失，則流體速度vn,f和結構速度vn,s相等

vn,f-vn,s=0

(22)

并且有

vn,s=iωCfsu

(23)

(24)

式中，ρf為流體密度。將式(24)代入到邊界元系統(tǒng)方程式(18)中，可以得到

Hp=Gω2ρfCfsu+pin

(25)

將式(21)代入式(15)，可以得到結構動力學問題的耦合方程

Kdu=Csfp+fs

(26)

將式(26)代入方程式(25)以消除未知位移u，可得

[H-GY]p=Gqs+pin

(27)

式中：Y為稠密矩陣，可看作為結構部分等效全局阻抗矩陣；qs為不考慮聲荷載時由結構載荷作用產生的輻射聲通量。Y和qs的表達式分別為

(28)

奇異值分解或模態(tài)分析法都可用來近似求解方程式(28)中矩陣Kd的逆。獲得Kd的逆之后，使用一些常用的迭代求解器(例如GMRES)來求解式(27)可獲得耦合面上的聲壓值p。

5 數值算例

5.1 球模型

本節(jié)考慮的是水下球殼結構在點激勵載荷作用下振動輻射聲場分析，載荷作用在(r,0,0)點處，r為球殼半徑，如圖6所示。結構形狀尺寸參數與材料屬性參數如表一所示。該模型具有如下解析表達式[31]

(29)

式中：(R,θ)為測試點的極坐標；Pn為n階勒讓德多項式；hn為第n階球漢克爾函數；zn為輻射阻抗；Zn為模態(tài)阻抗。

圖6 水下球體模型示意圖Fig.6 Schematic diagram of underwater sphere model

表1 球模型材料參數表

圖7給出40 Hz時位于xoy平面上半徑為10 m的圓環(huán)上的若干測試點處的聲壓幅值數值解與解析解的對比。圖7中符號為“Sub-IGABEM”表示基于細分曲面的等幾何分析算法得到的數值解，“Analytical”代表解析解。觀察圖7可以發(fā)現數值解與解析解吻合，驗證了本文算法的正確性。圖8給出位于(2r,0,0)處的測試點的聲壓值隨頻率變化圖，觀察該圖發(fā)現數值解與解析解吻合，再次驗證本文算法的正確性。

圖7 40 Hz時圓環(huán)上計算點處的聲壓值Fig.7 Sound pressure at calculation points at 40 Hz

圖8 計算點處聲壓隨頻率變化Fig.8 Sound pressure at a computing point with frequency

接下來測試POD與RBF結合近似物理場響應的準確度。以激勵球模型為例，保持其他參數不變，選取球殼半徑r∈(0.5 m,3.5 m)為輸入變量, 位于xoy平面上半徑為10m的圓環(huán)上均勻分布的20個測試點處的聲壓值(單位：Pa)為響應值，考察該輸入變量對響應結果的影響。在輸入變量的取值范圍上均勻取樣(輸入變量密度函數為均勻密度函數)，取樣步長分別為0.4，0.2，0.1和0.05(單位：m)，產生的樣本點數分別為9,16,31和60。因此，構造出的Λ矩陣的空間維數分別為20×9，20×16，20×31和20×60。采用POD對前3個小規(guī)模Λ矩陣進行分解，并采用RBF近似獲得步長為0.05時的樣本點處的響應值，并分別與解析解對比，如圖9所示。圖9顯示測試結果與解析解符合，表明POD與RBF的組合可以準確預測任意輸入變量的響應結果。

表2給出7個輸入變量的響應結果的預測值與解析解的相對誤差。(選擇的7個輸入半徑均為預測點而非樣本點)。觀察該表2可發(fā)現：相對誤差整體較小，驗證了POD與RBF組合進行響應結果預測的可靠性；插值步長越小(原始樣本點數越多)，預測結果越精確。需要注意的是，雖然增加原始樣本點數可以提高預測結果的精度，但是需要對每個原始樣本點進行響應分析，這在實際工程中往往會提高計算成本。因此，選取合適的樣本點數來進行預測是有必要的。

圖9 在極坐標為(10,0.35)的點處的聲壓值隨球殼半徑變化Fig.9 Sound pressure at a point with polar coordinate (10,0.35) in terms of radius

表2 不同插值步長下的近似值與解析值的對比

使用奇異值分解獲得的矩陣Σ是對角矩陣，主對角線上元素從大到小排列且迅速降低。較大的奇異值包含較大的系統(tǒng)特征，根據這個特點，我們可以對矩陣進行降階處理。選取輸入變量為r，步長為0.1時構造的響應矩陣(維數20×31)進行奇異值分解，然后再通過RBF操作插值得到若干樣本點處的響應。對Σ矩陣分別進行3種操作，方案1:不改變Σ矩陣的維數；方案2:Σ矩陣的維數選取為min(20,31)； 方案3:Σ矩陣的維數選取為0.5×min(20,31)。

圖10和11分別給出了坐標為(10,0,0)點處的聲壓值隨輸入變量r變化時的期望與標準差。圖10、圖11中，x軸表示用于MCs計算期望與標準差的樣本點數(預測點數)，方案4表示解析解。觀察這兩個圖可以發(fā)現，不同方案預測效果不同，但都隨著預測點數的增加而逼近解析解。雖然方案1和方案2具有相同的計算表現，但方案2的矩陣維數要比方案1小，因此預測效率更高。方案3的降階度最高，對應的預測精度也最差，但在40個預測點時也收斂到解析解。

接下來我們考慮輸入變量r(單位：m)服從高斯分布時測試點處聲壓的期望與標準差結果。高斯分布函數的均值為μ=3，方差為σ=0.1，選取25個輸入變量樣本來預測50個預測點處的響應值，然后通過這50個預測點響應值進行基于MCs的期望與標準差計算。預測得到的坐標為(10,0,0)點處的聲壓值的期望為2.41×10-3Pa，標準差為3.09×10-3Pa，與解析解吻合。將殼的厚度t(單位：m)設定為輸入變量，服從高斯分布，均值為μ=0.03，方差為σ=0.01，其他變量與前面一致。選取25個輸入變量樣本來預測50個預測點處的響應值，并與解析解對比，如圖12所示，可以發(fā)現數值解與解析解吻合。然后通過這50個預測點響應值進行基于MCs的期望與標準差計算，得到的期望數值解為2.009×10-3Pa，標準差數值解為3.962×10-4Pa，也與解析解一致。這兩個例子中，樣本點減少50%，從而減少了一半的網格劃分和耦合算法的迭代，且對原始數據進行了數據降維，減少了數據存儲所消耗的內存，直接提高了計算效率。

圖10 期望值隨樣本點數變化Fig.10 Expectations in terms of number of samples

圖11 標準差隨樣本點數變化Fig.11 Standard deviation with number of samples

圖12 變量為t時的預測結果與解析解對比Fig.12 Comparison of prediction results and analytical solutions with variable t

5.2 復雜薄壁艇模型

本節(jié)對復雜結構進行隨機性分析，使用個人PC端電腦處理器為i7-8700CPU。分析結構采用薄壁艇模型，坐標原點位于艇艏與艇體連接處截面圓圓心，x軸正方向由艇艏指向艇艉，初始模型由9 510個三角形單元構成，在艇艏頂端處受到沿x軸負方向，幅值為1 N的簡諧點激勵作用。通過細分曲面技術構造出光滑結構模型，如圖13所示。圖14給出40 Hz時，沿著x軸方向傳播的平面波作用下潛艇表面聲壓值分布，分布圖呈現對稱性，單個響應下的計算時間為30 min。

在該激勵頻率作用下，對結構振動響應隨殼厚度變化的不確定性進行分析研究，考察殼厚度變化對測試點處聲壓值的影響，20個測試點均勻分布于yoz面內半徑10 m的圓環(huán)上，如圖13所示。把艇殼的厚度t(單位：m)，滿足期望μ=0.01，方差σ=0.005的高斯分布)作為設計變量。從該范圍內均勻的選取出140個樣本點，使用耦合算法計算這些變量下的聲壓值響應，總的計算時間為4 200 min。隨后通過降階快速算法，均勻的從該140個樣本點中選取25個作為訓練點，插值計算出140個變量下(隨機分布)的響應?？焖偎惴偟挠嬎銜r間為750 min。

圖15顯示使用25個原始樣本數據預測140個樣本的響應結果。從圖15中可以看出，使用POD和RBF得到的預測數據與使用細分曲面等幾何數值方法獲得的結果吻合，驗證本文算法的可靠性。

圖13 細分曲面潛艇模型Fig.13 Submarine model constructed by subdivision surfaces

圖14 在40 Hz時潛艇表面聲壓分布Fig.14 Sound pressure contour on structural surface at 40 Hz

圖15 聲壓隨殼厚度變化Fig.15 Sound pressure in terms of thickness

6 結 論

本文基于等幾何有限元和邊界元方法，建立了聲振耦合系統(tǒng)的分析模型，使用細分曲面法對復雜模型進行建模，并對水下振動的復雜結構的輻射聲場進行了準確預報。此外，對由模型的輸入不確定造成的輸出不確定，進行基于蒙特卡羅模擬的隨機性分析。數值算例結果表明，本文提出的算法能夠對系統(tǒng)響應較為準確的預測，并可通過加速蒙特卡羅模擬來實現隨機變量下響應結果的統(tǒng)計特征的快速計算。