參數(shù)不確定性內(nèi)聲場(chǎng)分析的封閉區(qū)間有限元方法

向育佳， 史治宇， 馮雪磊

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

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

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

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

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

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

?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)聲場(chǎng)的聲學(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分別為邊界上的輸入聲壓、速度激勵(lì)；A0為邊界導(dǎo)納,A0=1/Z0。

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

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

(5)

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

(6)

(7)

(8)

(9)

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

(10)

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

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

(11)

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

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

(12)

(13)

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

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

(18)

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

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

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

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

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

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

(19)

(20)

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

(21)

(22)

(23)

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

(24)

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

(25)

Λmne=blkdiag(Λk,Λm)

(26)

(27)

Bmne={Φk,Φm}T

(28)

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

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

(29)

(30)

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

2.2 Neumaier-Pownuk迭代計(jì)算方法

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

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

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

(31)

(32)

res(p*)=FI-DIp*

(33)

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

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

{λ}=diag(Λ·Δα)

(34)

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

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

(35)

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

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

(36)

否則，計(jì)算矩陣M的最大特征值對(duì)應(yīng)的特征向量，并以此作為初始向量w，重新計(jì)算步驟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計(jì)算循環(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通過(guò)循環(huán)中零長(zhǎng)度改變量和迭代預(yù)設(shè)步數(shù)，設(shè)置第一個(gè)“終止準(zhǔn)則”，即

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

(47)

步驟7計(jì)算每步的區(qū)間結(jié)果和殘差的幅值，并給出第二個(gè)“終止準(zhǔn)則”來(lái)加快循環(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的動(dòng)力學(xué)拓展

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

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

(51)

(52)

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

(53)

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

(54)

其中

(55)

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

(56)

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

(57)

(58)

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

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

4.1 二維管道聲場(chǎng)

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

圖5 二維管道聲場(chǎng)的有限元網(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à)的，不影響計(jì)算結(jié)果。并且，可以看出沒(méi)有一種區(qū)間攝動(dòng)有限元法可以給出保守性的解的邊界，總會(huì)出現(xiàn)蒙特卡羅解落在邊界外的情況。這是由于區(qū)間攝動(dòng)有限元法在分解不確定性和求區(qū)間系數(shù)矩陣的近似逆階段，分別在Taylor級(jí)數(shù)展開(kāi)式和Neumann級(jí)數(shù)展開(kāi)式中用有限階項(xiàng)去近似真實(shí)值，這樣不可避免的產(chǎn)生“非保守近似”的問(wèn)題。特別是對(duì)于動(dòng)力學(xué)拓展的SMW-IPFEM，相較于傳統(tǒng)的兩種Taylor級(jí)數(shù)分解不確定性的區(qū)間有限元法TN-MIPFEM和RSE-IPFEM，解邊界的更寬一些，在全頻段的表現(xiàn)也更穩(wěn)定。但是，RSE-IPFEM的效果要差很多，額外的誤差來(lái)自構(gòu)造過(guò)程中利用“模態(tài)疊加法”去計(jì)算響應(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

計(jì)算時(shí)間的比較，如表1所示。TN-MIPFEM顯然是最快的，RSE-IPFEM次之，Enclosing-IFEM僅慢于前兩種IPFEMs。而且Enclosing-IFEM保持和EBE-IFEM 同樣的精度，卻只有EBE-IFEM不到一半的計(jì)算時(shí)間。因?yàn)镸NE組裝自由度為729，EBE組裝自由度為2 560，MNE組裝可以顯著的降低EBE組裝帶來(lái)的龐大的迭代計(jì)算自由度，Enclosing-IFEM可以被認(rèn)為是EBE-IFEM的降階技術(shù)。SMW-IPFEM的計(jì)算時(shí)間稍長(zhǎng)于Enclosing-IFEM，因?yàn)榇罅康摹爸纫痪仃嚒钡姆纸膺\(yùn)算和累加求和運(yùn)算，相比于矩陣運(yùn)算顯得十分的低效。

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

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

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

圖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

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

表2 對(duì)于二維車內(nèi)聲場(chǎng)10 000次蒙特卡羅取樣，基于攝動(dòng) 和改進(jìn)的區(qū)間算術(shù)的區(qū)間有限元法的單頻求解時(shí)間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)的定義來(lái)評(píng)價(jià)區(qū)間解的精度，如下

(59)

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

圖11 不同區(qū)間有限元法對(duì)于二維管道聲場(chǎng)中心線 位置的相對(duì)性誤差的比較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)。通過(guò)將Enclosing-IFEM的區(qū)間動(dòng)力學(xué)矩陣等價(jià)變形為“并矢積”式，本文還解決了基于S-M-W級(jí)數(shù)的區(qū)間攝動(dòng)有限元方法(SMW-IPFEM)的動(dòng)力學(xué)分析的拓展問(wèn)題，將SMW-IPFEM由結(jié)構(gòu)靜力學(xué)分析領(lǐng)域推廣到一般動(dòng)力學(xué)研究領(lǐng)域?；诓淮_定性聲場(chǎng)的算例分析，本文得到以下結(jié)論：

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

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

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

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