向育佳， 史治宇， 馮雪磊

(南京航空航天大學(xué) 航空學(xué)院 機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗室，南京 210016)

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

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

1 參數(shù)不確定性內(nèi)聲場的控制方程

1.1 含有區(qū)間變量的聲場動力學(xué)平衡方程

在無內(nèi)聲源穩(wěn)態(tài)內(nèi)聲場的控制方程為“亥姆霍茲方程”即

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

(1)

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

圖1 內(nèi)聲場的聲學(xué)邊界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為邊界導(dǎo)納,A0=1/Z0。

將聲學(xué)域有限元離散，并引入形函數(shù)N，式(1)經(jīng)過加權(quán)余數(shù)法(伽遼金方法)的等價積分變換轉(zhuǎn)化為“穩(wěn)態(tài)”動力學(xué)方程

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

(5)

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

(6)

(7)

(8)

(9)

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

(10)

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

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

(11)

1.2 區(qū)間單元矩陣單峰量分解方法

不同于傳統(tǒng)IPFEM采用泰勒級數(shù)展開法[18]來分解不確定性，改進(jìn)的區(qū)間數(shù)學(xué)方法引入結(jié)構(gòu)靜力分析中的“單峰量分解”方法可以消除有限階泰勒級數(shù)展開式的“截斷誤差”。根據(jù)EBE-IFEM的聲學(xué)單元矩陣分解為“并矢積”形式，對于二維問題，單元剛度陣和質(zhì)量陣的并矢積表達(dá)式為

(12)

(13)

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

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

(14)

(15)

(16)

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

(17)

1.3 混合節(jié)點(diǎn)-單元組裝方法

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

受到結(jié)構(gòu)靜力學(xué)領(lǐng)域的區(qū)間分析方法SMW-IPFEM的啟發(fā)，本文提出一種新的“混合節(jié)點(diǎn)-單元”組裝方法，不需要額外的拉格朗日乘子約束，對于執(zhí)行“單峰量分解”的單元矩陣在組裝成系統(tǒng)矩陣后仍滿足“并矢積”形式，例如對剛度陣組裝，即

(18)

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

圖3 矩陣Φk混合節(jié)點(diǎn)單元組裝方法Fig.3 The mixed-nodal-element assembly for Φk

2 封閉區(qū)間有限元法

2.1 區(qū)間動力學(xué)方程的“迭代封閉式”

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

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

(19)

(20)

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

(21)

(22)

(23)

(2) 若系統(tǒng)無阻尼，則系數(shù)矩陣為

(24)

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

(25)

Λmne=blkdiag(Λk,Λm)

(26)

(27)

Bmne={Φk,Φm}T

(28)

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

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

(29)

(30)

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

2.2 Neumaier-Pownuk迭代計算方法

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

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

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

(31)

(32)

res(p*)=FI-DIp*

(33)

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

步驟1對于聲學(xué)系統(tǒng)輸入的A,B,Λ, ΔαI,Dc,F,δI, 有限元法產(chǎn)生的動力學(xué)系數(shù)陣Dc一定是可逆的，簡化表示對角陣{λ}

{λ}=diag(Λ·Δα)

(34)

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

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

(35)

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

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

(36)

否則，計算矩陣M的最大特征值對應(yīng)的特征向量，并以此作為初始向量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計算循環(huán)初值H0,Y0

H0=w″/w′

(41)

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

(42)

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

(43)

(44)

(45)

(46)

步驟6通過循環(huán)中零長度改變量和迭代預(yù)設(shè)步數(shù)，設(shè)置第一個“終止準(zhǔn)則”，即

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

(47)

步驟7計算每步的區(qū)間結(jié)果和殘差的幅值，并給出第二個“終止準(zhǔn)則”來加快循環(huán)速度

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

(48)

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

(49)

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

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

(50)

3 SMW-IPFEM的動力學(xué)拓展

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

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

(51)

(52)

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

(53)

利用SMW級數(shù)近似式，區(qū)間系數(shù)矩陣的近似逆陣為

(54)

其中

(55)

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

(56)

這樣根據(jù)攝動法的“區(qū)間單調(diào)性”假設(shè)，可以算出聲場的區(qū)間響應(yīng)

(57)

(58)

4 二維不確定性內(nèi)聲場的數(shù)值算例

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

4.1 二維管道聲場

算例一為如圖5所示的二維管道聲場，尺寸為1.0 m×0.1 m，除了最左端為速度邊界并受到1 m/s的法向速度激勵，其余邊界為剛性壁面。管道聲場內(nèi)充滿空氣，聲速和空氣的密度為區(qū)間量，即：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)來預(yù)測聲壓的響應(yīng)邊界，模型被離散成640個單元，729個節(jié)點(diǎn)。為了比較不同IFEMs結(jié)果的精度以及計算效率，結(jié)合確定性聲學(xué)有限元法進(jìn)行10 000次蒙特卡羅樣本試驗，文獻(xiàn)[24]表明，當(dāng)隨機(jī)樣本足夠多時，蒙特卡羅結(jié)果的集合可以被視為真實(shí)解的參照解。

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

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

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

圖7 中心線300 Hz的聲壓響應(yīng)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可以被認(rèn)為是EBE-IFEM的降階技術(shù)。SMW-IPFEM的計算時間稍長于Enclosing-IFEM，因為大量的“秩一矩陣”的分解運(yùn)算和累加求和運(yùn)算，相比于矩陣運(yùn)算顯得十分的低效。

表1 對于二維管道聲場10 000次蒙特卡羅取樣，基于攝動 和改進(jìn)的區(qū)間算術(shù)的區(qū)間有限元法的單頻求解時間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 車內(nèi)二維聲場

圖8 二維車內(nèi)聲場的有限元網(wǎng)格Fig.8 The finite element mesh for 2D interior-acoustic of a car

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

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

圖10 R1位置10～230 Hz的聲壓響應(yīng)邊界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，該方法在效率上不如攝動方法。但是，同蒙特卡羅的真實(shí)解相比，Enclosing-IFEM在計算效率和保守性的邊界都是無法替代的唯一方案。

表2 對于二維車內(nèi)聲場10 000次蒙特卡羅取樣，基于攝動 和改進(jìn)的區(qū)間算術(shù)的區(qū)間有限元法的單頻求解時間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 相關(guān)性誤差分析

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

(59)

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

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

5 結(jié) 論

本文基于提出的“混合節(jié)點(diǎn)-單元”(MNE)組裝方式，發(fā)展出了一種新的基于改進(jìn)區(qū)間數(shù)學(xué)策略的區(qū)間有限元方法，即封閉區(qū)間有限元法(Enclosing-IFEM)。通過將Enclosing-IFEM的區(qū)間動力學(xué)矩陣等價變形為“并矢積”式，本文還解決了基于S-M-W級數(shù)的區(qū)間攝動有限元方法(SMW-IPFEM)的動力學(xué)分析的拓展問題，將SMW-IPFEM由結(jié)構(gòu)靜力學(xué)分析領(lǐng)域推廣到一般動力學(xué)研究領(lǐng)域?；诓淮_定性聲場的算例分析，本文得到以下結(jié)論：

(1) Enclosing-IFEM在不損失改進(jìn)區(qū)間數(shù)學(xué)方法的保守性解的計算優(yōu)勢情況下，將傳統(tǒng)EBE-IFEM的計算效率提高了一倍。

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

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

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