參數不確定性內聲場分析的封閉區間有限元方法

向育佳， 史治宇， 馮雪磊

(南京航空航天大學 航空學院 機械結構力學及控制國家重點實驗室，南京 210016)

內聲場分析是傳統的確定性分析方法[1]的重要研究內容，但是聲場響應會受到環境及溫度的影響。由于流體的聲學參數(如：密度ρ0、聲速c)以及聲學邊界的參數(如：邊界速度v0，邊界阻抗Z0)在其取值區間內波動，經典的聲學有限元等方法無法對這類參數不確定的聲場建模分析。隨著區間分析[2-3]理論的提出，通過建立“不確定但有界”輸入參數和輸出結果的上、下邊界的映射關系，區間有限元方法(interval finite-element method, IFEM)成為一種高效的分析不確定問題數值方法。

針對不確定性聲場分析問題，只能依靠第二類方法來獲取容納更多真實解的保守結果。為了提高求解速度，本文基于“混合節點-單元”(mixed-nodal-element, MNE)組裝方式，提出了一種新的封閉區間有限元方法(enclosing-interval finite element method,Enclosing-IFEM)去預測不確定性聲學問題的響應，在結合Enclosing-IFEM的前處理方式，又將基于S-M-W(Sherrman-Morrison-Woodbury)級數的區間攝動有限方法[14](Sherrman-Morrison-Woodbury series based interval perturbation finite element method,SMW-IPFEM)拓展到了動力學、聲學分析中。最后，結合兩個算例，通過與其他區間有限元方法，以及蒙特卡羅法(TN-MIPFEM(Taylor-Neumann series based modified interval perturbation finite element method)[15-16], SMW-IPFEM, 有理級數區間攝動有限元法RSE-IPFEM(rational series expansion based interval perturbation finite element method)[17])結果的對比驗證說明了提出的Enclosing-IFEM的計算優勢，最后得出Enclosing-IFEM在獲得和EBE-IFEM一致的計算結果時計算效率提高了一倍。

1 參數不確定性內聲場的控制方程

1.1 含有區間變量的聲場動力學平衡方程

在無內聲源穩態內聲場的控制方程為“亥姆霍茲方程”即

?2p(x,y,z)+k2·p(x,y,z)=0

(1)

式中：k為波數,k=ω/c；ω為角頻率；c為聲速。如圖1所示，一般聲學域V的邊界條件分為3種：Dirichlet邊界(聲壓邊界)Ωp，Neumann邊界(速度邊界)Ωv，Robin邊界(阻抗邊界)，即

圖1 內聲場的聲學邊界Fig.1 Boundary of interior acoustic domain

p=p0onΩp

(2)

vn=j/(ρ0ω)·?p/?n=v0onΩv

(3)

p=vn/A0=j/(ρ0ωA0)·?p/?nonΩZ

(4)

式中：p0,v0分別為邊界上的輸入聲壓、速度激勵；A0為邊界導納,A0=1/Z0。

將聲學域有限元離散，并引入形函數N，式(1)經過加權余數法(伽遼金方法)的等價積分變換轉化為“穩態”動力學方程

(K-ω2M+jωC)p=F

(5)

在僅考慮速度激勵和邊界阻抗情況下，剛度陣K、質量陣M、阻尼陣C以及載荷向量F可以表示為

(6)

(7)

(8)

(9)

如果假設密度ρ、聲速c，邊界速度v0，邊界導納Z0為區間參數，并用上、下邊界來描述，那么系統的不確定性即可以表示為用一個r維封閉的區間向量

(10)

式中： 上標I為區間量； Δ為區間半徑；uc為區間向量中點值；上、下劃線分別為區間量的上、下邊界。因此，考慮了區間參數的動力學方程式(5)可以改寫成

DI(uI)pI(uI)=FI(uI)DI=(KI-ω2MI+jωCI)

(11)

1.2 區間單元矩陣單峰量分解方法

不同于傳統IPFEM采用泰勒級數展開法[18]來分解不確定性，改進的區間數學方法引入結構靜力分析中的“單峰量分解”方法可以消除有限階泰勒級數展開式的“截斷誤差”。根據EBE-IFEM的聲學單元矩陣分解為“并矢積”形式，對于二維問題，單元剛度陣和質量陣的并矢積表達式為

(12)

(13)

對于阻尼邊界上的單元i，如圖2所示，阻尼陣的并矢積形式為

圖2 邊界上單元的單峰量分解Fig.2 Unimodal components’ decomposition of elements on the acoustic boundary

(14)

(15)

(16)

同理，速度邊界上的單元j，由式(9)可以分解為

(17)

1.3 混合節點-單元組裝方法

盡管EBE-IFEM在內聲場分析問題上已經取得了最佳的保守結果，由于逐單元(element-by-element,EBE)的組裝方式，不僅擴增了系統的自由度(不同單元的公共節點被分裂為多個節點，并從屬于不同單元)，而且為了強制約束節點連續性，還引入了拉格朗日乘子。因此，EBE-IFEM 在迭代求解時，中間變量階數過高，占用了過多的內存，計算時間相較于其他IPFEMs過長。

受到結構靜力學領域的區間分析方法SMW-IPFEM的啟發，本文提出一種新的“混合節點-單元”組裝方法，不需要額外的拉格朗日乘子約束，對于執行“單峰量分解”的單元矩陣在組裝成系統矩陣后仍滿足“并矢積”形式，例如對剛度陣組裝，即

(18)

其中，Φk的組裝如圖3所示。

圖3 矩陣Φk混合節點單元組裝方法Fig.3 The mixed-nodal-element assembly for Φk

2 封閉區間有限元法

2.1 區間動力學方程的“迭代封閉式”

通過第1章的“單峰量分解”和“混合節點-單元”組裝，內聲場的區間動力學方程(式(11))的系數矩陣DI可以改寫成并矢積的形式，按照系統有無阻尼，可以分為兩種情況討論。

(1) 若考慮阻尼則系數矩陣需要實部、虛部分離[19]

(19)

(20)

Λmne=blkdiag(Λk,Λm,Λc,Λk,Λm,Λc)

(21)

(22)

(23)

(2) 若系統無阻尼，則系數矩陣為

(24)

Amne={Φk,-ω2Φm}

(25)

Λmne=blkdiag(Λk,Λm)

(26)

(27)

Bmne={Φk,Φm}T

(28)

不同于傳統的直接法[20](如：高斯消去法，組合法，不等式法等)對于區間線性方程組的形式沒有要求，迭代求解方法只能適用于“迭代封閉式”的區間線性方程組。通過并矢積的變換，聲學系統的區間動力學方程(式(11))可以轉化為“迭代封閉式”[21]。

[Dc+Adiag(Λ·ΔαI)B]pI=FδI

(29)

(30)

式中： mid(·)為取區間變量中間值； rad(·)為區間變量減去中間值后的正負波動區間。

2.2 Neumaier-Pownuk迭代計算方法

“封閉迭代式”的求解方法一般分為兩種，一種是基于“不動點理論”[22]的“Rump’s方法”，另一種是Neumaier在嚴格數學證明基礎上建立的“Neumaier-Pownuk 方法”(N-P法)。Neumaier等的研究已經證明，N-P法在區間分析的工程實踐中殘差收斂性更好，結果的“過估”問題更小。本文采用N-P法求解式(29)，其流程圖如圖4所示。

圖4 Neumaier-Pownuk方法的流程圖Fig.4 The flowchart of the Neumaier-Pownuk method

為說明該方法的步驟，分別定義區間量uI的零長度和幅值以及區間線性方程組的殘差，如下

(31)

(32)

res(p*)=FI-DIp*

(33)

式中，p*為pI的近似解。

步驟1對于聲學系統輸入的A,B,Λ, ΔαI,Dc,F,δI, 有限元法產生的動力學系數陣Dc一定是可逆的，簡化表示對角陣{λ}

{λ}=diag(Λ·Δα)

(34)

步驟2以單位列向量w為初始向量，試函數為

w′=w-rad({λ})mag[B(Dc)-1A]

(35)

步驟3判斷步驟2的向量w′中是否存在正數的元素，若存在，則

w″={λ}mag[B(Dc)-1FδI]

(36)

否則，計算矩陣M的最大特征值對應的特征向量，并以此作為初始向量w，重新計算步驟2～步驟3，即

M=rad({λ})mag[B(Dc)-1A]

(37)

w=eigvect(M)i,ρi=max[eigval(M)]

(38)

w′=w-rad({λ})mag[B(Dc)-1A]

(39)

w″={λ}mag[B(Dc)-1FδI]

(40)

步驟4計算循環初值H0,Y0

H0=w″/w′

(41)

Y0=B(Dc)-1(FδI+AH0)

(42)

步驟5代入步驟4的初值，并執行循環體：式(43)～式(46)，n=1,2，…

(43)

(44)

(45)

(46)

步驟6通過循環中零長度改變量和迭代預設步數，設置第一個“終止準則”，即

max[wid(Hn)-wid(Hn-1)]>ε或n

(47)

步驟7計算每步的區間結果和殘差的幅值，并給出第二個“終止準則”來加快循環速度

(pI)n=(Dc)-1(FδI+AHn) Δdn=mag{res[(pI)n]}

(48)

‖Δdn-Δdn-1‖∞≤‖η·Δdn‖∞

(49)

步驟8滿足步驟6～步驟7之一則輸出近似的區間解

p*=(pI)n=(Dc)-1(FδI+AHn)

(50)

3 SMW-IPFEM的動力學拓展

自Impollonia和Muscolino基于S-M-W級數提出SMW-IPFEM以來，因為動力學控制方程的系數矩陣難以轉化為“并矢積”形式，SMW-IPFEM只能被用于求解結構靜力學問題。Xia等因此尋求替代方案，根據SMW-IPFEM的核心思想，即在Neumann級數展開式中盡可能的保留更多的高階項，發展出改進方法TN-MIPFEM，從而改善備受詬病的Neumann級數求區間逆陣的非保守近似的缺陷。然而，Taylor級數分解不確定性的非保守近似問題在TN-MIPFEM中同樣存在，這使得基于“單峰量分解”的SMW-IPFEM成為區間攝動有限元方法(IPFEMs)中控制“級數展開式非保守近似問題”最優途徑。

根據第2章Enclosing-IFEM的區間動力學方程的構造方法，即式(19)～式(30)，顯然“并矢積”形式的動力學系數陣可以成為SMW-IPFEM的動力學推廣的重要前提，那么系數可以改寫為“秩一矩陣”和的形式

(51)

(52)

A=[a1,…,ar],B=[b1,…,br]

(53)

利用SMW級數近似式，區間系數矩陣的近似逆陣為

(54)

其中

(55)

忽略式(54)中k階展開式的交叉項，式(54)可以進一步表示為

(56)

這樣根據攝動法的“區間單調性”假設，可以算出聲場的區間響應

(57)

(58)

4 二維不確定性內聲場的數值算例

為了比較不同方法的計算精度以及效率，本章的所有算法程序都運行在MATLAB R2017a的軟件平臺上，計算機配置：主頻為2.90 GHz的Intel酷睿CPU i5-9400F，對于EBE-IFEM和Enclosing-IFEM需要調用區間運算工具箱INTLAB V12的區間數學函數庫[23]。

4.1 二維管道聲場

算例一為如圖5所示的二維管道聲場，尺寸為1.0 m×0.1 m，除了最左端為速度邊界并受到1 m/s的法向速度激勵，其余邊界為剛性壁面。管道聲場內充滿空氣，聲速和空氣的密度為區間量，即：cI=[334.3,343.2]m/s,ρI0=[1.204,1.269]kg/m3。采用兩類共5種不同的IFEMs(TN-MIPFEM, RSE-IPFEM, SMW-IPFEM, EBE-IFEM 以及Enclosing-IFEM)來預測聲壓的響應邊界，模型被離散成640個單元，729個節點。為了比較不同IFEMs結果的精度以及計算效率，結合確定性聲學有限元法進行10 000次蒙特卡羅樣本試驗，文獻[24]表明，當隨機樣本足夠多時，蒙特卡羅結果的集合可以被視為真實解的參照解。

圖5 二維管道聲場的有限元網格Fig.5 The finite element mesh for 2D acoustic tube

圖6和圖7中，兩種改進的區間數學方法Enclosing-IFEM和EBE-IFEM的響應在全頻率完全吻合，表明本文提出的“混合節點單元”(MNE)組裝方法和逐單元(EBE)的組裝方法是等價的，不影響計算結果。并且，可以看出沒有一種區間攝動有限元法可以給出保守性的解的邊界，總會出現蒙特卡羅解落在邊界外的情況。這是由于區間攝動有限元法在分解不確定性和求區間系數矩陣的近似逆階段，分別在Taylor級數展開式和Neumann級數展開式中用有限階項去近似真實值，這樣不可避免的產生“非保守近似”的問題。特別是對于動力學拓展的SMW-IPFEM，相較于傳統的兩種Taylor級數分解不確定性的區間有限元法TN-MIPFEM和RSE-IPFEM，解邊界的更寬一些，在全頻段的表現也更穩定。但是，RSE-IPFEM的效果要差很多，額外的誤差來自構造過程中利用“模態疊加法”去計算響應。

圖6 中心線150 Hz的聲壓響應的虛部Fig.6 The imaginary part of sound pressure along the central-axis at f=150 Hz

圖7 中心線300 Hz的聲壓響應Fig.7 The imaginary part of sound pressure along the central-axis at f=300 Hz

計算時間的比較，如表1所示。TN-MIPFEM顯然是最快的，RSE-IPFEM次之，Enclosing-IFEM僅慢于前兩種IPFEMs。而且Enclosing-IFEM保持和EBE-IFEM 同樣的精度，卻只有EBE-IFEM不到一半的計算時間。因為MNE組裝自由度為729，EBE組裝自由度為2 560，MNE組裝可以顯著的降低EBE組裝帶來的龐大的迭代計算自由度，Enclosing-IFEM可以被認為是EBE-IFEM的降階技術。SMW-IPFEM的計算時間稍長于Enclosing-IFEM，因為大量的“秩一矩陣”的分解運算和累加求和運算，相比于矩陣運算顯得十分的低效。

表1 對于二維管道聲場10 000次蒙特卡羅取樣，基于攝動 和改進的區間算術的區間有限元法的單頻求解時間Tab.1 Time-consuming of 10 000 MC, perturbation based and modified interval arithmetic based IFEMs for 2D single frequency of 2D acoustic tube

4.2 車內二維聲場

圖8 二維車內聲場的有限元網格Fig.8 The finite element mesh for 2D interior-acoustic of a car

首先，根據圖9的單頻分析結果可得，Enclosing-IFEM和EBE-IFEM計算結果吻合。而且，比其他3種區間攝動有限元法，結果能囊括真實解的替代解(蒙特卡羅解集)，體現出改進的區間數學方法的保守性優勢。與文獻結論一致，“N-P方法”的結果“過估”問題也被控制在較小范圍內。再者，從圖10中也可以得到同樣的驗證，Enclosing-IFEM和EBE-IFEM的保守性在全頻段都保持的很好，其他3種IPFEMs的解的帶寬在共振頻率附近過窄，如f=130 Hz和f=180 Hz附近的頻響，攝動方法的非保守近似很明顯。

圖9 底部輪廓線位置200 Hz的聲壓響應邊界Fig.9 Bounds of the sound pressure along the bottom-line at f=200 Hz

圖10 R1位置10～230 Hz的聲壓響應邊界Fig.10 Bounds of the sound pressure at R1 in the frequency band f=10-230 Hz

計算效率，如表2所示。由于自由度降低了，Enclosing-IFEM相比于EBE-IFEM縮短了一半時間，相較于TN-MIPFEM和RSE-IPFEM，該方法在效率上不如攝動方法。但是，同蒙特卡羅的真實解相比，Enclosing-IFEM在計算效率和保守性的邊界都是無法替代的唯一方案。

表2 對于二維車內聲場10 000次蒙特卡羅取樣，基于攝動 和改進的區間算術的區間有限元法的單頻求解時間Tab.2 Time-consuming of 10 000 MC, perturbation based and modified interval arithmetic based IFEMs for 2D single frequency of 2D acoustic field of a car

4.3 相關性誤差分析

根據其他文獻的誤差分析研究[25]，本文給出上、下邊界相關性誤差(ER)的定義來評價區間解的精度，如下

(59)

式中，pMC為蒙特卡羅解集，在真實解無法給出時可以用來替代真實解。由于本文給出的對比方法比較多，故將不同方法相比于蒙特卡羅解的相對誤差刻畫在同一圖中。選取第一個算例作為分析對象，如圖11所示，改進的區間數學方法Enclosing-IFEM和EBE-IFEM的相關性誤差最小，且小于1，而其他方法的誤差很大。這表明，不同于IPFEMs，基于改進的區間數學方法Enclosing-IFEM和EBE-IFEM確實能將不確定性和“過估”現象控制在逐單元水平。

圖11 不同區間有限元法對于二維管道聲場中心線 位置的相對性誤差的比較Fig.11 Comparison of the relative errors of responses of the 2D acoustic tube along the central-axis by different IFEMs

5 結 論

本文基于提出的“混合節點-單元”(MNE)組裝方式，發展出了一種新的基于改進區間數學策略的區間有限元方法，即封閉區間有限元法(Enclosing-IFEM)。通過將Enclosing-IFEM的區間動力學矩陣等價變形為“并矢積”式，本文還解決了基于S-M-W級數的區間攝動有限元方法(SMW-IPFEM)的動力學分析的拓展問題，將SMW-IPFEM由結構靜力學分析領域推廣到一般動力學研究領域。基于不確定性聲場的算例分析，本文得到以下結論：

(1) Enclosing-IFEM在不損失改進區間數學方法的保守性解的計算優勢情況下，將傳統EBE-IFEM的計算效率提高了一倍。

(2) Enclosing-IFEM的結果與EBE-IFEM一致，“混合節點單元”組裝方法和“逐單元”組裝方法是等價的。

(3) 在不確定性分解方面，Taylor級數展開的方法存在非保守近似，通過SMW-IPFEM與TN-MIPFEM在管道聲場的計算比較，采用“單峰量分解”的SMW-IPFEM精度較高。

(4) 盡管SMW-IPFEM同樣可以應用于聲學分析，由于“秩一矩陣”的分解及求和十分低效，其計算時間是TN-MIPFEM的100～200多倍。