基于直接微分法的水下結(jié)構(gòu)聲學(xué)快速敏感度分析

陳磊磊， 趙文暢， 陳海波， 肖琪聃

(1.信陽(yáng)師范學(xué)院 土木工程學(xué)院，河南 信陽(yáng) 464000； 2.中國(guó)科學(xué)技術(shù)大學(xué) 近代力學(xué)系，合肥 230026)

基于直接微分法的水下結(jié)構(gòu)聲學(xué)快速敏感度分析

陳磊磊1， 趙文暢2， 陳海波2， 肖琪聃1

(1.信陽(yáng)師范學(xué)院 土木工程學(xué)院，河南 信陽(yáng) 464000； 2.中國(guó)科學(xué)技術(shù)大學(xué) 近代力學(xué)系，合肥 230026)

在降低結(jié)構(gòu)振動(dòng)噪聲方面，結(jié)構(gòu)聲學(xué)優(yōu)化方法展現(xiàn)了很高的應(yīng)用潛力，并越來(lái)越受到研究者們的重視。通過(guò)結(jié)構(gòu)聲學(xué)敏感度分析建立輻射聲場(chǎng)與結(jié)構(gòu)材料屬性和結(jié)構(gòu)形狀參數(shù)的關(guān)系，結(jié)合優(yōu)化算法可設(shè)計(jì)出低輻射噪聲和高聲隱身性能的水下航行結(jié)構(gòu)。發(fā)展耦合有限元法與快速多極高階非連續(xù)邊界元法進(jìn)行流固耦合問(wèn)題的敏感度分析，克服傳統(tǒng)耦合算法的計(jì)算精度低，內(nèi)存占有量高，難以進(jìn)行大規(guī)模實(shí)際問(wèn)題的計(jì)算等問(wèn)題。設(shè)計(jì)變量可選為流體密度、結(jié)構(gòu)密度、泊松比、楊氏模量、殼厚度和形狀參數(shù)等。針對(duì)不同設(shè)計(jì)變量，采用直接微分法分別推導(dǎo)出結(jié)構(gòu)聲學(xué)系統(tǒng)和輻射聲功率敏感度表達(dá)式。通過(guò)帶有解析解的水下點(diǎn)激勵(lì)球殼算例驗(yàn)證算法的正確性與有效性，并通過(guò)簡(jiǎn)化潛艇算例展示該算法在大規(guī)模實(shí)際問(wèn)題上的應(yīng)用潛力。

輻射聲功率；敏感度分析；耦合有限元與邊界元；非連續(xù)邊界元

水中航行結(jié)構(gòu)振動(dòng)輻射噪聲不僅影響結(jié)構(gòu)內(nèi)部人員的身心健康，也對(duì)海底生物造成干擾。降低振動(dòng)噪聲和提高水下結(jié)構(gòu)的聲隱身性能是一項(xiàng)非常重要的研究課題。進(jìn)行水中結(jié)構(gòu)(尤其是薄殼結(jié)構(gòu))的振動(dòng)響應(yīng)分析，流體附加作用不可忽略，同時(shí)也增加了數(shù)值計(jì)算的難度?？紤]到有限元在結(jié)構(gòu)振動(dòng)響應(yīng)分析和邊界元在無(wú)限域外聲場(chǎng)分析時(shí)的各自優(yōu)點(diǎn)，國(guó)內(nèi)外許多學(xué)者采用耦合有限元與邊界元方法(FE/BE)進(jìn)行流固耦合分析[1-4]。然而現(xiàn)有研究大多采用傳統(tǒng)邊界元進(jìn)行聲場(chǎng)分析，計(jì)算量大，內(nèi)存占有量高，難以進(jìn)行大規(guī)模實(shí)際問(wèn)題的數(shù)值計(jì)算。文獻(xiàn)[5-6]使用快速多極算法(Fast Multipole Method,FMM)加速流體邊界元計(jì)算，提高了流固耦合問(wèn)題的計(jì)算效率，降低了內(nèi)存占有量。然而這些研究大多采用常量邊界元進(jìn)行聲場(chǎng)計(jì)算，為了提高計(jì)算精度，加大結(jié)構(gòu)離散，導(dǎo)致新的計(jì)算困難產(chǎn)生。相對(duì)于已經(jīng)被廣泛使用的連續(xù)線性和二次邊界元，高階非連續(xù)邊界元不僅計(jì)算精度高而且避免處理角點(diǎn)問(wèn)題，尤其在耦合系數(shù)矩陣、濕表面面積陣及其逆矩陣的計(jì)算上，更加簡(jiǎn)單和有效[7]。本文結(jié)合高階非連續(xù)邊界元與快速多極算法形成快速多極高階非連續(xù)邊界元進(jìn)行結(jié)構(gòu)聲學(xué)計(jì)算，提高傳統(tǒng)耦合算法的計(jì)算精度和計(jì)算效率。另外使用單Helmholtz邊界積分方程進(jìn)行外聲場(chǎng)計(jì)算時(shí)，會(huì)產(chǎn)生解的非唯一性問(wèn)題，本文使用Burton-Miller方法克服這一問(wèn)題[8]。

通過(guò)敏感度分析研究水中振動(dòng)結(jié)構(gòu)輻射聲場(chǎng)與結(jié)構(gòu)本身的材料屬性和形狀參數(shù)的關(guān)系，可指導(dǎo)結(jié)構(gòu)材料的選取和形狀設(shè)計(jì)，以降低輻射噪聲和提高水下結(jié)構(gòu)的聲隱身性能。文獻(xiàn)[9-10]使用快速多極邊界元法進(jìn)行聲學(xué)敏感度分析，沒(méi)有考慮聲振耦合，文獻(xiàn)[11]采用耦合FEM/BEM(Boundary Element Method)方法進(jìn)行聲振耦合計(jì)算，然而在進(jìn)行敏感度分析時(shí)忽略了耦合作用。實(shí)際上對(duì)聲振耦合方程進(jìn)行敏感度分析是非常困難的，既要進(jìn)行結(jié)構(gòu)有限元敏感度計(jì)算，又要進(jìn)行聲學(xué)邊界元敏感度計(jì)算，同時(shí)還需考慮耦合矩陣的敏感度結(jié)果，因此現(xiàn)有研究很少進(jìn)行考慮聲振耦合的敏感度分析。本文克服這些困難，對(duì)此進(jìn)行深入研究，推導(dǎo)出結(jié)構(gòu)聲學(xué)耦合敏感度表達(dá)式，設(shè)計(jì)變量可選取為流體密度、結(jié)構(gòu)密度、泊松比、楊氏模量和形狀尺寸參數(shù)等。針對(duì)不同設(shè)計(jì)變量分別推導(dǎo)出結(jié)構(gòu)振動(dòng)輻射聲功率敏感度表達(dá)式，進(jìn)而找到目標(biāo)函數(shù)與設(shè)計(jì)變量的對(duì)應(yīng)關(guān)系，為基于敏感度分析的結(jié)構(gòu)聲學(xué)優(yōu)化設(shè)計(jì)奠定基礎(chǔ)。最后使用帶有解析解的簡(jiǎn)單算法驗(yàn)證本文算法的正確性和有效性，使用簡(jiǎn)化潛艇模型展現(xiàn)本文算法在大規(guī)模實(shí)際問(wèn)題中的應(yīng)用潛力。

1 結(jié)構(gòu)聲學(xué)耦合系統(tǒng)

1.1 水下結(jié)構(gòu)聲學(xué)系統(tǒng)

組合水下結(jié)構(gòu)振動(dòng)有限元方程與流體邊界元方程得到表達(dá)式

(1)

對(duì)于實(shí)際工程問(wèn)題，直接求解式(1)，可得到位移響應(yīng)u和聲壓響應(yīng)p，然而計(jì)算量和內(nèi)存占有量非常大，難以收斂，并且計(jì)算精度不高。為了提高耦合算法的計(jì)算效率，將式(1)中的結(jié)構(gòu)振動(dòng)有限元方程代入到流體邊界元方程中，形成耦合邊界元表達(dá)式[7]

Hp-GWCsfp=pi+GWfs

(2)

其中，

W=ω2ρfS-1CfsD-1

D=K-ω2M+iωC

(3)

觀察式(2)發(fā)現(xiàn)存在D-1項(xiàng)，實(shí)際上直接求解結(jié)構(gòu)矩陣的逆是非常耗時(shí)的，本文在計(jì)算式(2)右邊D-1fs時(shí)，先通過(guò)求解線性系統(tǒng)方程Dx=fs，得到x向量即為D-1fs的結(jié)果。另外，本文使用GMRES迭代方法求解式(2)，假設(shè)第k次迭代步對(duì)應(yīng)的節(jié)點(diǎn)聲壓向量為pk,在計(jì)算式(2)左邊D-1Csfpk時(shí)，先計(jì)算Csfpk=b，然后通過(guò)求解線性系統(tǒng)方程Dx=b獲得x向量以得到D-1Csfpk的結(jié)果。在每次迭代過(guò)程中該系統(tǒng)方程都需要被重新求解，由于D矩陣是稀疏帶狀的，使用稀疏直接求解器(Sparse Direct Solver)能十分有效且精確地求解這一系統(tǒng)方程。實(shí)際上在每次迭代計(jì)算中最耗時(shí)的是式(2)的求解，如果使用傳統(tǒng)邊界元算法進(jìn)行該方程的求解，對(duì)于N自由度問(wèn)題(流體邊界元自由度為N)，計(jì)算復(fù)雜度為O(N2)，內(nèi)存占有量也為O(N2)，這無(wú)疑是非常耗時(shí)和耗內(nèi)存的，限制了該方法在大規(guī)模實(shí)際問(wèn)題中的應(yīng)用。本文采用寬頻快速多極算法加速傳統(tǒng)邊界元算法中的矩陣向量相乘，同時(shí)避免計(jì)算流體邊界元系數(shù)矩陣G與H，對(duì)于N自由度問(wèn)題，計(jì)算復(fù)雜度為O(N)，內(nèi)存占有量也為O(N)，降低了計(jì)算量和內(nèi)存占有量，推進(jìn)了耦合有限元與邊界元法在大規(guī)模實(shí)際問(wèn)題中的應(yīng)用。

1.2 水下結(jié)構(gòu)聲學(xué)系統(tǒng)敏感度分析

在進(jìn)行敏感度分析時(shí)，首先需要確定設(shè)計(jì)變量，本文針對(duì)不同設(shè)計(jì)變量，例如流體密度、結(jié)構(gòu)密度、泊松比、楊氏模量、結(jié)構(gòu)形狀尺寸等，分別給出相應(yīng)的敏感度計(jì)算表達(dá)式。在下面的計(jì)算過(guò)程中，v代表設(shè)計(jì)變量。

(1) 當(dāng)選取流體密度ρf為設(shè)計(jì)變量，直接將式(2)對(duì)設(shè)計(jì)變量ρf進(jìn)行微分，可得到如下敏感度表達(dá)式

(4)

(3) 當(dāng)選取結(jié)構(gòu)形狀參數(shù)為設(shè)計(jì)變量，直接將式(2)對(duì)設(shè)計(jì)變量進(jìn)行微分，可得

(5)

a=WCsfp+Wfs

(6)

使用式(2)得到節(jié)點(diǎn)聲壓向量值p，然后利用式(4)或式(5)可以得到節(jié)點(diǎn)聲壓向量敏感度值。觀察式(6)可發(fā)現(xiàn)存在D矩陣求逆后再求導(dǎo)的計(jì)算，實(shí)際上直接計(jì)算D-1已經(jīng)是十分耗時(shí)的，再求導(dǎo)更是非常困難和耗時(shí)。為了避免直接計(jì)算D-1的導(dǎo)數(shù)，本文推導(dǎo)出如下表達(dá)式來(lái)代替D-1的導(dǎo)數(shù)

(7)

本文采用有限差分法來(lái)計(jì)算式(7)中的D的導(dǎo)數(shù)項(xiàng)，設(shè)定的有限差分步長(zhǎng)為10-4。如果使用傳統(tǒng)邊界元方法進(jìn)行數(shù)值計(jì)算，會(huì)形成系數(shù)矩陣G與H，采用有限差分法可以得到其敏感度結(jié)果，然而正如前文所說(shuō)，該方法計(jì)算效率低，實(shí)際上對(duì)于G或H的敏感度矩陣與向量的相乘，類似于G或H與向量的相乘，可以使用快速多極算法進(jìn)行加速計(jì)算，無(wú)需直接計(jì)算G、H及其敏感度矩陣，提高了計(jì)算效率，降低了內(nèi)存占有量，為大規(guī)模實(shí)際耦合問(wèn)題的敏感度計(jì)算及其優(yōu)化分析提供了有效的數(shù)值方法。

使用邊界元法可得到流體域內(nèi)任一點(diǎn)處的聲壓p

(8)

式中：G(x,y)=iρfωeikr/(4πr);F(x,y)為G(x,y)在點(diǎn)x處的法向?qū)?shù)值;y為域內(nèi)計(jì)算點(diǎn);S為流固耦合面;r為點(diǎn)x與y的距離。使用式(8)可以得到流體域內(nèi)任意點(diǎn)處的聲壓值，將式(8)對(duì)設(shè)計(jì)變量進(jìn)行微分，可得到流體域內(nèi)任意點(diǎn)處的聲壓敏感度值。

1.3 水下結(jié)構(gòu)振動(dòng)輻射聲功率的敏感度計(jì)算

在流體中，圍繞振動(dòng)結(jié)構(gòu)的任意面上的輻射聲功率表達(dá)式為

(9)

(10)

式中，v為流固耦合面上流體法向速度向量，與結(jié)構(gòu)節(jié)點(diǎn)位移向量u具有如下轉(zhuǎn)換關(guān)系

v=-iωS-1Cfsu

(11)

使用式(2)得到節(jié)點(diǎn)聲壓向量p，然后利用式(1)中的第一式(結(jié)構(gòu)有限元方程)，得到結(jié)構(gòu)節(jié)點(diǎn)位移響應(yīng)u，最后將式(11)代入到式(10)中，即可得到流固耦合面上的結(jié)構(gòu)振動(dòng)輻射聲功率。

將式(10)對(duì)設(shè)計(jì)變量進(jìn)行微分得到如下表達(dá)式

(12)

(13)

式中:()H為向量取共軛后再轉(zhuǎn)置;w1與w2表達(dá)式為

w1=iωCfsu*

(14)

耦合系數(shù)矩陣Cfs由結(jié)構(gòu)形狀參數(shù)決定，與流體密度、結(jié)構(gòu)材料屬性等參數(shù)沒(méi)有關(guān)系。因此，當(dāng)結(jié)構(gòu)形狀參數(shù)被選為設(shè)計(jì)變量時(shí)，Cfs的導(dǎo)數(shù)值存在，否則為0。利用式(4)或式(5)可以得到節(jié)點(diǎn)聲壓向量敏感度值，為了得到輻射聲功率的敏感度值，還需計(jì)算得到節(jié)點(diǎn)位移的敏感度值，下面針對(duì)不同設(shè)計(jì)變量，分別推導(dǎo)出u的敏感度表達(dá)式。

(1) 當(dāng)選取流體密度ρf為設(shè)計(jì)變量，直接將式(1)中的第一式對(duì)設(shè)計(jì)變量ρf進(jìn)行微分，可得到如下敏感度表達(dá)式

(15)

(2) 當(dāng)選取結(jié)構(gòu)密度ρs、泊松比v、楊氏模量E、殼厚度等參數(shù)為設(shè)計(jì)變量時(shí)，得到

(16)

(3) 當(dāng)選取結(jié)構(gòu)形狀參數(shù)為設(shè)計(jì)變量時(shí)(例如球殼半徑)，得到

(17)

使用式(15)～式(17)可以得到節(jié)點(diǎn)位移敏感度值，然后利用式(12)即可得到輻射聲功率敏感度值。

2 數(shù)值算例

2.1 數(shù)值計(jì)算離散單元

在數(shù)值計(jì)算中，先使用有限元軟件ANSYS建立和離散結(jié)構(gòu)模型，導(dǎo)出每個(gè)單元的剛度、質(zhì)量矩陣，然后使用一維變帶寬存儲(chǔ)方法進(jìn)行所有單元矩陣的整合，最后形成整體剛度和質(zhì)量矩陣，由于結(jié)構(gòu)矩陣的稀疏與對(duì)稱性，只需存儲(chǔ)整體矩陣的上三角非零元素，使用3個(gè)一維數(shù)組可有效存儲(chǔ)整體矩陣，第一個(gè)數(shù)組存儲(chǔ)整體矩陣上三角每行非零元素的數(shù)量，第二個(gè)數(shù)組存儲(chǔ)上三角每行非零元素的列編號(hào)，第三個(gè)數(shù)組存儲(chǔ)上三角每行非零元素的值。同時(shí)，采用相同網(wǎng)格進(jìn)行結(jié)構(gòu)有限元與流體邊界元的計(jì)算。使用Shell281殼單元進(jìn)行水下薄殼結(jié)構(gòu)有限元計(jì)算，在每個(gè)節(jié)點(diǎn)處有3個(gè)位移自由度和3個(gè)旋轉(zhuǎn)自由度；為了克服傳統(tǒng)連續(xù)邊界單元的計(jì)算精度低，角點(diǎn)處處理困難等問(wèn)題，本文使用帶有二次幾何形狀近似的非連續(xù)線性邊界單元BE94進(jìn)行流體邊界元的計(jì)算，BE94單元表示使用9節(jié)點(diǎn)進(jìn)行幾何插值計(jì)算，使用內(nèi)部4節(jié)點(diǎn)進(jìn)行物理插值計(jì)算，詳細(xì)的單元描述和插值過(guò)程請(qǐng)參考文獻(xiàn)[7]。

2.2 水下點(diǎn)激勵(lì)球殼算例

本節(jié)進(jìn)行水下薄球殼在單點(diǎn)激勵(lì)時(shí)的振動(dòng)響應(yīng)分析，球殼半徑為r，坐標(biāo)原點(diǎn)在球心處，單位幅值的點(diǎn)激勵(lì)作用在點(diǎn)(r,0,0)上，詳細(xì)模型圖見(jiàn)文獻(xiàn)[7]，流體與結(jié)構(gòu)參數(shù)如表1所示。文獻(xiàn)[7]給出了該模型振動(dòng)位移響應(yīng)與輻射聲壓的解析表達(dá)式。由于該模型沿X軸軸對(duì)稱分布，流固耦合面上的振動(dòng)輻射聲功率可表達(dá)成

(18)

式中:θm=mπ/n;v*f(θm)=iωu*(θm)，足夠大的n能產(chǎn)生高精度的計(jì)算結(jié)果。將上式對(duì)設(shè)計(jì)變量進(jìn)行微分可得到輻射聲功率敏感度的解析表達(dá)。

表1 水下球殼的材料與幾何參數(shù)

在圖1中，“FE88/FMBE94”表示使用Shell281單元(8節(jié)點(diǎn))進(jìn)行結(jié)構(gòu)有限元分析，使用帶有二次幾何形狀近似的非連續(xù)線性單元BE94進(jìn)行流體邊界元計(jì)算，同時(shí)使用FMM算法進(jìn)行耦合邊界元方程的快速求解?！癋E44/FMBE41”表示使用Shell63單元(4節(jié)點(diǎn))進(jìn)行結(jié)構(gòu)有限元分析，使用帶有線性幾何形狀近似的非連續(xù)常量單元BE41進(jìn)行流體邊界元計(jì)算，同時(shí)使用FMM算法進(jìn)行耦合邊界元方程的快速求解。設(shè)計(jì)變量分別為流體密度和結(jié)構(gòu)密度，計(jì)算頻率為50 Hz,計(jì)算點(diǎn)分布在距離球心R=2r的圓環(huán)上，離散單元數(shù)為384。觀察這兩個(gè)圖發(fā)現(xiàn)，高階耦合單元的計(jì)算結(jié)果與解析解符合的比較好，然而低階單元的誤差比較大，驗(yàn)證了本文算法的高精度特點(diǎn)。

(a) 流體密度

(b) 結(jié)構(gòu)密度

Fig.1 Radiated sound pressure sensitivity at distanceR=2rat 50 Hz

在圖2中，后綴“BM”表示使用Burton-Miller方法進(jìn)行流體邊界元計(jì)算，“CBIE”表示使用傳統(tǒng)單Helmholtz邊界積分方程進(jìn)行流體聲場(chǎng)計(jì)算。解析解曲線的計(jì)算頻率步長(zhǎng)為0.01 Hz，數(shù)值解步長(zhǎng)為2 Hz，離散單元數(shù)為864。觀察該圖可以發(fā)現(xiàn)在100 Hz以內(nèi)的水下結(jié)構(gòu)本征頻率依次為55.9 Hz、70.5 Hz、80.6 Hz、88.4 Hz和94.8 Hz。另外，使用“CBIE”方法在148 Hz左右得到的數(shù)值結(jié)果偏離解析解，實(shí)際上是由于傳統(tǒng)單Helmholtz邊界積分方程在計(jì)算外聲場(chǎng)問(wèn)題時(shí)會(huì)產(chǎn)生非唯一解的問(wèn)題，然而使用Burton-Miller法可以有效的解決這個(gè)問(wèn)題。

圖2 在點(diǎn)(2r,0,0)處的聲壓隨計(jì)算頻率變化圖

圖3呈現(xiàn)了流固耦合面上的輻射聲功率敏感度幅值隨計(jì)算頻率變化關(guān)系，設(shè)計(jì)變量依次選取為流體密度、結(jié)構(gòu)密度、泊松比、楊氏模量、球殼厚度和半徑。觀察這些圖可以發(fā)現(xiàn)在低頻處敏感度幅值比較小，在結(jié)構(gòu)本征頻率處迅速變化，在一些頻率處急速降低(表現(xiàn)為近似向下垂直線)是由于敏感度值在該頻率附近迅速變化并且發(fā)生算法符號(hào)的改變，例如由正變負(fù)，或由負(fù)變正。另外，在這些圖中，數(shù)值解與解析解都符合的比較好，驗(yàn)證了本文算法的正確性和有效性。

(a) 流體密度

(b) 結(jié)構(gòu)密度

(c) 泊松比

(d) 楊氏模量

(e) 殼厚度

(f) 球半徑

Fig.3 Amplitude of radiated sound power sensitivity on the fluid-structure interaction surface in terms of frequency

圖4給出了在計(jì)算頻率為50 Hz時(shí)使用傳統(tǒng)算法與快速算法在計(jì)算流固耦合面上的輻射聲功率值時(shí)所耗費(fèi)時(shí)間的對(duì)比。其中“剛性解”和“加速剛性解”分別表示使用傳統(tǒng)邊界元和寬頻快速多極算法加速計(jì)算剛性球模型時(shí)所需的時(shí)間(不考慮流固耦合)；“彈性解”和“加速?gòu)椥越狻狈謩e表示使用傳統(tǒng)耦合有限元與邊界元法和耦合有限元與寬頻快速多極邊界元算法計(jì)算彈性球模型時(shí)所需的時(shí)間。由圖4可發(fā)現(xiàn)，網(wǎng)格數(shù)越多，寬頻快速多極算法的加速效果越明顯。另外，耦合分析比非耦合分析耗費(fèi)更多時(shí)間，這是因?yàn)樵隈詈戏治龅牡?jì)算中，每一步都需要求解方程Dx=b，并且彈性耦合的迭代收斂次數(shù)遠(yuǎn)大于剛性分析次數(shù)。

圖4 計(jì)算時(shí)間隨單元數(shù)變化關(guān)系

表2給出了采用三種不同類型耦合單元在相近流體邊界節(jié)點(diǎn)自由度下(約1 536)計(jì)算輻射聲功率時(shí)的精度對(duì)比，表中的誤差數(shù)據(jù)表示百分誤差，選取5個(gè)頻率進(jìn)行數(shù)值計(jì)算。表中“FE44/FMCBE44”表示采用線性有限殼單元(Shell63)進(jìn)行結(jié)構(gòu)分析，線性連續(xù)邊界單元進(jìn)行聲場(chǎng)分析；“FE88/FMCBE88”表示采用二次有限殼單元(Shell281)進(jìn)行結(jié)構(gòu)分析，二次連續(xù)邊界單元進(jìn)行聲場(chǎng)分析。觀察表2可以發(fā)現(xiàn)FE44/FMCBE44單元表現(xiàn)最差，F(xiàn)E88/FMBE94單元表現(xiàn)最好。另外在80 Hz本征頻率處計(jì)算誤差非常大，在200 Hz以外計(jì)算誤差隨頻率迅速增加，總之耦合非連續(xù)單元表現(xiàn)最好。

表2 三種不同離散耦合單元計(jì)算精度對(duì)比

Tab.2 Comparison of calculation accuracy for three different types of coupled elements %

2.3 簡(jiǎn)化潛艇結(jié)構(gòu)聲學(xué)響應(yīng)分析

簡(jiǎn)化潛艇的形狀與尺寸如圖5所示，分別進(jìn)行在平面波作用下(沿著正X軸方向傳播，幅值為1 Pa)和在模型兩端的A點(diǎn)與B點(diǎn)處施加集中力F作用下(幅值為1 N)的結(jié)構(gòu)聲學(xué)響應(yīng)分析，流體密度，結(jié)構(gòu)材料屬性參數(shù)和殼厚度見(jiàn)表1。

圖5 兩端分別用錐形殼和半球殼封閉的圓柱殼模型

Fig.5 An elastic cylindrical shell with a half-spherical shell and a conical shell

首先進(jìn)行平面波入射時(shí)的聲散射分析。圖6給出了在計(jì)算點(diǎn)(-50,0,0)處的聲壓隨計(jì)算頻率的變化關(guān)系，“剛性解”表明不考慮流固耦合作用，僅采用FMBE94流體單元進(jìn)行聲學(xué)計(jì)算，得到流體域內(nèi)計(jì)算點(diǎn)處的聲壓；“彈性解”是采用FE88/FMBE94耦合單元進(jìn)行流固耦合計(jì)算。觀察該圖可以發(fā)現(xiàn)，剛性解與數(shù)值解的偏差隨著頻率增加而增大，表明進(jìn)行水下薄殼結(jié)果的振動(dòng)響應(yīng)分析時(shí)必須考慮流體的附加作用。

圖6 在計(jì)算點(diǎn)(-50,0,0)處的聲壓隨頻率變化圖

Fig.6 Sound pressure at a computing point (-50,0,0) in terms of frequency

圖7呈現(xiàn)了在計(jì)算點(diǎn)處聲壓敏感度幅值隨頻率變化關(guān)系，殼厚度選為設(shè)計(jì)變量。使用快速算法得到的計(jì)算結(jié)果與傳統(tǒng)算法得到的計(jì)算結(jié)果吻合一致，表明快速算法的使用沒(méi)有降低計(jì)算精度，保持了傳統(tǒng)邊界元方法的高計(jì)算精度優(yōu)點(diǎn)。另外，在進(jìn)行平面波作用下的結(jié)構(gòu)散射聲場(chǎng)敏感度計(jì)算時(shí)，若不考慮流固耦合，只須求解聲學(xué)敏感度邊界積分方程，若考慮流固耦合須求解耦合邊界元敏感度方程，增加的計(jì)算困難就是耦合矩陣與結(jié)構(gòu)矩陣敏感度的高效求解，尤其當(dāng)高階耦合單元被使用，更是增加了這一困難。本文使用ANSYS軟件進(jìn)行結(jié)構(gòu)矩陣的計(jì)算，并采用有限差分法計(jì)算結(jié)構(gòu)矩陣的敏感度結(jié)果。

圖7 在計(jì)算點(diǎn)處對(duì)殼厚度的聲壓敏感度幅值隨頻率變化

Fig.7 Sound pressure sensitivity with respect to the thickness of the shell at the computing point

然后考察在集中力作用下，該模型的振動(dòng)響應(yīng)。圖8和圖9分別給出了流固耦合面上的輻射聲功率和其敏感度幅值隨頻率變化關(guān)系，設(shè)計(jì)變量為殼厚度。同樣的，剛性解與彈性解差別很大，實(shí)際上結(jié)構(gòu)越薄，耦合影響就越大。

圖8 流固耦合面上輻射聲功率隨頻率變化

Fig.8 Radiated sound power on the fluid-structure interaction surface in terms of frequency

圖9 流固耦合面上輻射聲功率敏感度幅值隨頻率變化

Fig.9 Amplitude of radiated sound power sensitivity on the fluid-structure interaction surface

3 結(jié) 論

本文首先推導(dǎo)出用于進(jìn)行流固耦合分析的耦合邊界元方程和結(jié)構(gòu)振動(dòng)輻射聲功率表達(dá)式，然后針對(duì)不同設(shè)計(jì)變量分別推導(dǎo)出相應(yīng)的結(jié)構(gòu)聲學(xué)敏感度和結(jié)構(gòu)振動(dòng)輻射聲功率敏感度表達(dá)式，通過(guò)帶有解析解的點(diǎn)激勵(lì)球殼算例驗(yàn)證本文算法的有效性和正確性，通過(guò)簡(jiǎn)化潛艇算例驗(yàn)證本文算法的應(yīng)用潛力，下一步對(duì)整個(gè)算法過(guò)程與程序編寫進(jìn)行優(yōu)化，然后對(duì)大規(guī)模實(shí)際問(wèn)題進(jìn)行快速敏感度分析和結(jié)構(gòu)聲學(xué)優(yōu)化算法研究。

[1] 俞孟薩, 史小軍, 陳克勤. 采用有限元和邊界元方法分析彈性加肋圓柱殼的聲學(xué)相似性[J]. 中國(guó)造船, 1999(3): 65-71.

YU Mengsa, SHI Xiaojun, CHEN Keqin. Study on acoustic similitude of elastic stiffened cylindrical shells by FEM and BEM[J]. Shipbuilding of China, 1999(3): 65-71.

[2] 黎勝, 趙德有.用耦合有限元/邊界元方法研究加筋板的聲傳輸[J].振動(dòng)工程學(xué)報(bào), 2001, 14(3)：364-367.

LI Sheng, ZHAO Deyou. Coupled FEM/BEM approach for analysis of sound transmission through stiffened plates[J].Journal of Vibration Engineering, 2001, 14(3)：364-367.

[3] 沈順根, 李琪華, 王大云, 等. 加肋旋轉(zhuǎn)殼結(jié)構(gòu)噪聲聲輻射水彈性研究[J].中國(guó)造船, 1992(2): 53-62.

SHEN Shungen, LI Qihua, WANG Dayun, et al. Hydroelastic analysis of sound radiation from a stiffened shell of revolution[J]. Shipbuilding of China, 1992(2): 53-62.

[4] ZHENG H, LIU G R, TAO J S. FEM/BEM analysis of diesel piston-slap induced ship hull vibration and underwater noise[J]. Applied Acoustics, 2001, 62(4): 341-58.

[5] 陳磊磊, 陳海波,鄭昌軍,等.基于有限元與寬頻快速多極邊界元的二維流固耦合聲場(chǎng)分析[J].工程力學(xué), 2014, 31(8): 63-69.

CHEN Leilei, CHEN Haibo, ZHENG Changjun, et al. FEM/wideband FMBEM coupling analysis for two dimensional acoustic fluid-strucutre interaction problems[J]. Engineering Mechanics, 2014, 31(8)：63-69.

[6] JUNGE M, BRUNNER D, GAUL L. Solution of FE-BE coupled eigenvalue problems for the prediction of the vibro-acoustic behavior of shipe-like strucutres[J]. International Journal for Numerical Methods in Engineering, 2011, 87: 664-676.

[7] PETERS H, MARBURG S, KESSISSOGLOU N. Structural-acoustic coupling on non-conforming meshes with quadratic shape functions[J]. International Journal for Numerical Methods in Engineering, 2012, 91(1): 27-38.

[8] BURTON A J, MILLER G F. The application of integral equation methods to the numerical solution of some exterior boundary-value problems[J]. Proceedings of the Royal Society of London, 1971, 323: 201-210.

[9] CHEN L L, ZHENG C J, CHEN H B. A wideband FMBEM for 2D acoustic design sensitivity analysis based on direct differentiation method[J]. Computational Mechanics, 2013, 52(3): 631-648.

[10] ZHENG C J, MATSUMOTO T, TAKAHASHI T, et al. A wideband fast multipole boundary element method for three dimensional acoustic shape sensitivity analysis based on direct differentiation method[J]. Engineering Analysis with Boundary Elements, 2012, 36(3): 361-371.

[11] FRITZE D, MARBURG S, HARDTKE H J. FEM-BEM-coupling and structural-acoustic sensitivity analysis for shell geometries[J]. Computers and Structures, 2005, 82(2): 143-154.

Fast analysis for underwater structures’ acoustic sensitivity based on direct differentiation method

CHEN Leilei1, ZHAO Wenchang2, CHEN Haibo2, XIAO Qidan1

(1. School of Civil Engineering, Xinyang Normal University, Xinyang 464000, China；2. Department of Modern Mechanics, University of Science and Technology of China, Hefei 230026, China)

The structural-acoustic optimization method receives more attention of designers due to its high potential in reducing structures’ vibration and noise. The structural acoustic sensitivity analysis can be used to establish relations among radiated acoustic field, structural material properties and structural shape parameters, then unwater navigating structures with lower radiated noise and higher sound-stealth performance can be designed with the optimization method. Here, the coupled algorithm was developed by combining FE and fast multi-pole higher order discontinuous BE to perform the sensitivity analysis for fluid-structure interaction problems. The developed algorithm improved the calculation accuracy and decreased the memory storage. Fluid density, structural density, Poisson’s ratio, Yong’s modulus, structural shape size and so on were taken as design variables. The sensitivity formulas for structural-acoustic system and radiated sound power with respect to different design variables were derived based on the direct differentiation method, respectively. A numerical example for an underwater spherical shell under a point excitation with an analytical solution was used to verify the validity and correctness of the developed algorithm, and an example for a simplified submarine was used to show the application potential of the developed algorithm in large scale practical problems.

radiated sound power; sensitivity analysis; coupled finite element(FE) and boundary element(BE); discontinuous boundary element(BE)

國(guó)家自然科學(xué)基金(11172291;U1504505)；河南省科技攻關(guān)項(xiàng)目(172102210453)；河南省高等學(xué)校重點(diǎn)科研項(xiàng)目(16B560009 ；17A560009)

2016-04-07 修改稿收到日期：2016-05-24

陳磊磊 男，博士，講師，1986年11月生

O39

A

10.13465/j.cnki.jvs.2017.13.026