黃振衛(wèi)，周其斗，紀(jì) 剛，孟慶昌，段嘉希

(海軍工程大學(xué) 艦船工程系，武漢 430033)

加肋圓柱殼結(jié)構(gòu)的FE-BE算法網(wǎng)格尺度劃分原則研究

黃振衛(wèi)，周其斗，紀(jì) 剛，孟慶昌，段嘉希

(海軍工程大學(xué) 艦船工程系，武漢 430033)

研究結(jié)構(gòu)有限元耦合流體邊界元算法(FE-BE算法)中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分原則具有十分重要的意義。以有限長(zhǎng)加肋圓柱殼為研究對(duì)象，采用波數(shù)譜分析方法分離出結(jié)構(gòu)主振型分量(攜帶能量最多的振動(dòng)模式)。根據(jù)分析結(jié)果提出FE-BE算法中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分原則，并采用FE-BE算法驗(yàn)證了該原則的有效性。

加肋圓柱殼； 結(jié)構(gòu)有限元耦合流體邊界元； 振動(dòng)與聲輻射； 網(wǎng)格尺度

結(jié)構(gòu)有限元耦合流體邊界元(FE-BE)算法一直被認(rèn)為是解決水下大型結(jié)構(gòu)低頻振動(dòng)與聲輻射的有力工具[1-3]，因?yàn)椴捎眠吔缭椒ㄇ蠼鉄o限域的聲場(chǎng)問題時(shí)，邊界元只需要在結(jié)構(gòu)濕表面(結(jié)構(gòu)物與水接觸的表面)上劃分，當(dāng)結(jié)構(gòu)有限元的劃分完成后，邊界元?jiǎng)澐謽I(yè)已完成。為了節(jié)約FE-BE算法的計(jì)算成本，一般應(yīng)選取單元尺度較大的有限元單元?jiǎng)澐纸Y(jié)構(gòu)濕表面，但是這必然帶來計(jì)算精度的降低。另一方面，為了保證FE-BE算法的計(jì)算精度，結(jié)構(gòu)濕表面的單元尺度應(yīng)該盡量小以便于描述最短的結(jié)構(gòu)波長(zhǎng)，但是這將導(dǎo)致大型復(fù)雜結(jié)構(gòu)的邊界元求解規(guī)模急劇地增加，使得FE-BE算法的計(jì)算時(shí)間非常長(zhǎng)。因此，對(duì)于FE-BE算法，結(jié)構(gòu)濕表面單元尺度的選取與計(jì)算精度和計(jì)算成本之間有個(gè)權(quán)衡的問題。目前很少有文獻(xiàn)研究FE-BE算法中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分原則的問題。

結(jié)構(gòu)的振動(dòng)可以分解為一系列不同波長(zhǎng)的振動(dòng)分量的疊加，其中攜帶振動(dòng)能量最大的分量即為結(jié)構(gòu)的主振型分量。一般而言，采用有限元單元?jiǎng)澐纸Y(jié)構(gòu)濕表面時(shí)，單元尺度必須要能夠描述結(jié)構(gòu)濕表面上最短的彎曲波長(zhǎng)[4]，但是實(shí)際振動(dòng)結(jié)構(gòu)的最短彎曲波長(zhǎng)可以無限小，如果采用非常小的單元尺度劃分結(jié)構(gòu)濕表面，將會(huì)使得系統(tǒng)的邊界元求解規(guī)模急劇地增加。顯然，不能簡(jiǎn)單地采用結(jié)構(gòu)最短彎曲波長(zhǎng)作為FE-BE算法中濕表面網(wǎng)格劃分的參考標(biāo)準(zhǔn)。由于結(jié)構(gòu)主振型分量的彎曲波長(zhǎng)對(duì)于流固耦合計(jì)算精度有重要的作用，因此，本文重點(diǎn)采用波數(shù)譜方法分離有限長(zhǎng)圓柱殼的主振型分量，根據(jù)分析結(jié)果提出FE-BE算法中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分原則，并采用FE-BE算法驗(yàn)證了該原則的有效性。本文首次提出了以結(jié)構(gòu)主振型分量的彎曲波長(zhǎng)作為FE-BE算法中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分的參考標(biāo)準(zhǔn)，結(jié)論對(duì)于FE-BE算法的有限元模型的網(wǎng)格劃分具有十分重要的參考意義。

1 模型及有限元?jiǎng)澐?/h2>

以有限長(zhǎng)加肋圓柱殼為研究對(duì)象，結(jié)構(gòu)幾何尺寸及材料和流體屬性如表1所示。殼體內(nèi)表面上均勻布置23根肋骨，激勵(lì)力作用于中間肋骨上，沿著徑向指向殼體外部的方向。采用有限元軟件PATRAN進(jìn)行建模，結(jié)構(gòu)濕表面采用三角形單元進(jìn)行有限元?jiǎng)澐?，肋骨采用四邊形單元進(jìn)行有限元?jiǎng)澐郑永邎A柱殼的有限元模型示意圖如圖1所示。72個(gè)計(jì)算聲場(chǎng)點(diǎn)均布于與柱殼中間肋骨同心，半徑為9 m的圓周上。計(jì)算頻率為50 Hz～500 Hz，頻率間隔為5 Hz，模型為自由邊界條件，不考慮自由水面對(duì)結(jié)構(gòu)振動(dòng)與聲輻射的影響。

表1 模型幾何尺寸及材料和流體屬性Tab.1 Geometric size and parameter values for calculation

圖1 加肋圓柱殼有限元模型Fig.1 Finite element model of the ring stiffened cylindrical shell

2 主振型的分離

本節(jié)采用編制的圓柱殼結(jié)構(gòu)振動(dòng)波數(shù)譜程序，分離柱殼結(jié)構(gòu)主振型分量，該程序的正確性已在文獻(xiàn)[5]中得到很好的驗(yàn)證。為了描述有限長(zhǎng)圓柱殼結(jié)構(gòu)做剛體運(yùn)動(dòng)的波數(shù)譜，將柱殼長(zhǎng)度延長(zhǎng)為原柱殼長(zhǎng)度的12倍，假定延長(zhǎng)的部分的外殼體法向速度為0[5]，如圖2示，L為圓柱殼的長(zhǎng)度，a為圓柱殼的半徑，圖中省略了肋骨示意。

圖2 無限長(zhǎng)柱殼有限長(zhǎng)模型Fig.2 An infinite length cylindrical shell with finite vibration distribution

采用編制的FE-BE算法程序計(jì)算圓柱殼濕表面的法向位移復(fù)數(shù)幅值w(z,θ)，即

(1)

式中：wR和wI分別為是結(jié)構(gòu)濕表面法向位移的實(shí)部和虛部。值得注意的是，本節(jié)為了保證FE-BE算法程序計(jì)算結(jié)構(gòu)濕表面振動(dòng)響應(yīng)的精度，加肋圓柱殼的每檔肋骨間距采用6個(gè)單元進(jìn)行劃分，周向采用48個(gè)單元進(jìn)行劃分，最終得到15 552個(gè)結(jié)構(gòu)單元，8 330個(gè)節(jié)點(diǎn)，14 448個(gè)濕表面單元，7 226個(gè)濕表面節(jié)點(diǎn)。

將法向位移的實(shí)部wR沿著柱殼的軸向和周向兩個(gè)方向進(jìn)行傅里葉變換，得到實(shí)部所代表的駐波場(chǎng)

(2)

式中，kz為軸向波數(shù)，n為周向完整波的數(shù)量，且有

(3)

(4)

同理，可以得到虛部所代表的駐波場(chǎng)

(5)

總的法向速度的振動(dòng)功率可表示為

(6)

式中:Ev(n,kz)為法向速度振動(dòng)功率譜，表達(dá)式為：

Ev(n,kz)=

(7)

根據(jù)式(7)，即可繪制出不同的周向振動(dòng)模式n時(shí)，結(jié)構(gòu)法向速度振動(dòng)功率隨著結(jié)構(gòu)軸向波數(shù)的變化曲線，進(jìn)而分離出加肋圓柱殼的主振型分量以及主振型分量的彎曲波長(zhǎng)λb。

圖3為模型的法向速度振動(dòng)功率譜(n=0,1,…,9)，從圖3(a)中可以看出，當(dāng)n=2，無因次軸向波數(shù)大約為0時(shí)，結(jié)構(gòu)波的振動(dòng)能量最大(主振型分量)，對(duì)應(yīng)的結(jié)構(gòu)軸向波長(zhǎng)為∞。分析圓柱殼周向方向上結(jié)構(gòu)波的傳播特性，可以得到n=2時(shí)，周向上的結(jié)構(gòu)波長(zhǎng)為2πa/n=1.99 m。值得一提的是結(jié)構(gòu)中其它非主振型分量攜帶的能量小于主振型分量所攜帶的能量，例如圖3中振動(dòng)能量級(jí)大3 dB，意味著攜帶的能量將增大一倍左右。

采用相同的方法即可分離出其他頻率下結(jié)構(gòu)主振型分量，計(jì)算結(jié)果見表2。

(a) 50 Hz

(b) 300 Hz

(c) 500 Hz

表2 不同頻率下結(jié)構(gòu)主振型分量Tab.2 Main vibration mode at different frequency

3 網(wǎng)格尺度對(duì)FE-BE算法的影響

由于加肋圓柱殼的振動(dòng)與聲輻射缺乏理論解，本文將逐步加密結(jié)構(gòu)濕表面上的網(wǎng)格，觀察FE-BE算法的計(jì)算結(jié)果是否趨于穩(wěn)定，并與NASTRAN軟件中的FE-IE算法[6]計(jì)算結(jié)果進(jìn)行對(duì)比。為了保證無限元法的計(jì)算精度，在結(jié)構(gòu)外域構(gòu)造一個(gè)長(zhǎng)10 m，半徑5 m的圓柱面人工邊界，人工邊界的單元尺度為0.483 m，結(jié)構(gòu)濕表面軸向單元尺度為0.059 m，周向單元尺度為0.083 m，無限元的計(jì)算階數(shù)為10階。根據(jù)第2節(jié)波數(shù)譜的計(jì)算結(jié)果，設(shè)計(jì)工況如表3所示。為了便于研究，所有工況下柱殼圓周方向上的單元尺度選取為L(zhǎng)c=0.083 m，定義Nz=λb/Lz為一個(gè)主振型波長(zhǎng)內(nèi)的單元數(shù)量，其中，Lz為結(jié)構(gòu)濕表面軸向方向上的單元尺度，λb為結(jié)構(gòu)主振型分量軸向彎曲波波長(zhǎng)。

表3 計(jì)算工況Tab.3 Computational cases

圖4～圖6為各工況下輻射聲壓級(jí)的FE-BE算法計(jì)算值對(duì)比，從圖中可以看出：

(1) 激勵(lì)頻率為50 Hz時(shí)，不同結(jié)構(gòu)濕表面網(wǎng)格尺度下FE-BE算法的計(jì)算結(jié)果均與FE-IE算法的計(jì)算結(jié)果吻合良好。主要原因是，此時(shí)結(jié)構(gòu)主振型的軸向彎曲波波長(zhǎng)為無窮大，即使軸向單元尺度選取為一個(gè)肋骨間距一個(gè)單元(Lz=0.238 m)時(shí)，也能準(zhǔn)確描述主振型分量的彎曲波波長(zhǎng)。但是結(jié)構(gòu)濕表面的單元尺度(即邊界元尺度)不能無限增大，根據(jù)文獻(xiàn)[7]，邊界元的單元尺度一般應(yīng)保證一個(gè)聲波長(zhǎng)至少有6個(gè)單元，因此，結(jié)構(gòu)濕表面的單元尺度最大不能超過1/6倍聲波長(zhǎng)；

(2) 激勵(lì)頻率為300 Hz時(shí)，隨著結(jié)構(gòu)濕表面的單元尺度減小，F(xiàn)E-BE算法的計(jì)算結(jié)果趨于穩(wěn)定，且與FE-IE算法的計(jì)算結(jié)果吻合良好。主要原因是300 Hz時(shí)FE-BE算法的結(jié)構(gòu)濕表面網(wǎng)格尺度保證了一個(gè)主振型波長(zhǎng)內(nèi)至少有7個(gè)有限元單元，足以描述結(jié)構(gòu)的主振型分量的波動(dòng)特性；

(3) 激勵(lì)頻率為500 Hz時(shí)，隨著結(jié)構(gòu)濕表面的單元尺度減小，F(xiàn)E-BE算法的計(jì)算結(jié)果趨于穩(wěn)定，但與FE-IE算法的計(jì)算結(jié)果有一定的差別。當(dāng)結(jié)構(gòu)濕表面的軸向單元尺度為L(zhǎng)z=0.238 m時(shí)，一個(gè)主振型彎曲波長(zhǎng)內(nèi)只有3.2個(gè)單元，無法準(zhǔn)確地描述結(jié)構(gòu)主振型分量的波動(dòng)特性，因此，單元尺度為L(zhǎng)z=0.238 m的FE-BE算法計(jì)算結(jié)果與濕表面網(wǎng)格加密后FE-BE算法計(jì)算結(jié)果差別較大。當(dāng)結(jié)構(gòu)濕表面的軸向單元尺度為L(zhǎng)z=0.119 m時(shí)，保證了一個(gè)主振型彎曲波長(zhǎng)內(nèi)有6.5個(gè)有限元單元，此時(shí)FE-BE算法的計(jì)算結(jié)果與濕表面網(wǎng)格加密后FE-BE算法計(jì)算結(jié)果吻合良好。此外，500 Hz時(shí)網(wǎng)格加密后FE-BE算法計(jì)算結(jié)果與FE-IE算法計(jì)算結(jié)果產(chǎn)生一定差別的原因可能是近場(chǎng)效應(yīng)對(duì)FE-IE算法的影響更為明顯，為了改善FE-IE算法的計(jì)算精度，人工邊界應(yīng)該選取在距離結(jié)構(gòu)濕表面更遠(yuǎn)的位置或者采用更小的單元尺度劃分人工邊界。

圖4 加肋圓柱殼在50 Hz時(shí)的聲壓級(jí)Fig.4 Numerical SPL of the ring stiffened cylindrical shell at 50 Hz

圖5 加肋圓柱殼在300 Hz時(shí)的聲壓級(jí)Fig.5 Numerical SPL of the ring stiffened cylindrical shell at 300 Hz

圖6 加肋圓柱殼在500 Hz時(shí)的聲壓級(jí)Fig.6 Numerical SPL of the ring stiffened cylindrical shell at 500 Hz

總的來說，對(duì)加肋圓柱殼結(jié)構(gòu)而言，當(dāng)結(jié)構(gòu)主振型分量的彎曲波波長(zhǎng)小于對(duì)應(yīng)激勵(lì)頻率下的聲波波長(zhǎng)時(shí)，以結(jié)構(gòu)主振型分量的彎曲波長(zhǎng)作為FE-BE算法中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分標(biāo)準(zhǔn)是可行的，即結(jié)構(gòu)濕表面的網(wǎng)格尺度應(yīng)保證一個(gè)主振型分量的彎曲波長(zhǎng)內(nèi)至少有6個(gè)單元。值得注意的是，此時(shí)結(jié)構(gòu)濕表面的單元尺度能夠描述結(jié)構(gòu)中攜帶能量最大的振動(dòng)分量，進(jìn)而保證了FE-BE算法的計(jì)算精度。若將結(jié)構(gòu)濕表面的單元尺度進(jìn)一步減小或者增大，都有可能帶來對(duì)應(yīng)的計(jì)算成本的增加或計(jì)算精度的損失。

為了進(jìn)一步應(yīng)用FE-BE算法的網(wǎng)格劃分原則，本節(jié)將計(jì)算50～500 Hz時(shí)外域聲場(chǎng)中A(0 m,100 m,0 m)點(diǎn)的輻射聲壓級(jí)。根據(jù)以上分析，計(jì)算頻率在500 Hz以內(nèi)時(shí)，加肋圓柱殼每個(gè)肋骨間距采用2個(gè)單元進(jìn)行劃分(Lz=0.119 m)即可保證FE-BE算法的計(jì)算精度。同時(shí)，F(xiàn)E-IE算法中加肋圓柱殼每個(gè)肋骨間距4個(gè)有限元單元進(jìn)行劃分，人工邊界為一個(gè)長(zhǎng)12 m，半徑為8 m的圓柱面，人工邊界的單元尺度為0.48 m，無限元的計(jì)算階次為10階。A點(diǎn)的輻射聲壓級(jí)傳遞函數(shù)頻響曲線如圖7所示，從圖中可以看出，激勵(lì)頻率較低時(shí)，F(xiàn)E-BE算法的輻射聲壓級(jí)傳遞函數(shù)頻響曲線與FE-IE算法的頻響曲線吻合良好；激勵(lì)頻率較高時(shí)，F(xiàn)E-BE算法的輻射聲壓級(jí)傳遞函數(shù)頻響曲線與FE-IE算法的頻響曲線變化趨勢(shì)一致。

圖7 結(jié)構(gòu)聲學(xué)傳遞函數(shù)的FE-BE計(jì)算值Fig.7 Sound transfer function obtained by the FE-BE method

4 結(jié) 論

本文重點(diǎn)采用波數(shù)譜方法分離圓柱殼結(jié)構(gòu)的主振型分量，根據(jù)分析結(jié)果提出FE-BE算法中結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分原則，并采用FE-BE算法驗(yàn)證了該原則的有效性，計(jì)算結(jié)果表明：對(duì)于文中算例而言，當(dāng)結(jié)構(gòu)主振型分量的彎曲波長(zhǎng)小于對(duì)應(yīng)激勵(lì)頻率下的聲波長(zhǎng)時(shí)，以結(jié)構(gòu)主振型分量的彎曲波長(zhǎng)作為結(jié)構(gòu)濕表面的網(wǎng)格尺度劃分標(biāo)準(zhǔn)是可行的，即結(jié)構(gòu)濕表面的網(wǎng)格尺度應(yīng)保證一個(gè)主振型分量的彎曲波長(zhǎng)內(nèi)至少有6個(gè)單元。

[1] ZHOU Q, JOSEPH P F.A numerical method for the calculation of dynamic response and acoustic radiation from an underwater structure[J].Journal of Sound and Vibration,2005(283): 853-873.

[2] 紀(jì)剛,張緯康,周其斗.靜水壓力作用的水下結(jié)構(gòu)振動(dòng)及聲輻射[J].中國(guó)造船,2006,47(3):37-44.

JI Gang, ZHANG Weikang, ZHOU Qidou.Vibration and radiation from underwater structure considering the effect of static water preload[J].Shipbuilding of China,2006,47(3):37-44.

[3] EVERSTIN G C, HENDERSON F M.Coupled finite element/boundary element approach for fluid-structure interaction[J].J.Acoust.Soc.Am., 1990,87:1936-1947.

[4] 宗福開.波傳播問題中有限元分析的頻散特性及離散化準(zhǔn)則[J].爆炸與沖擊,1984,3(4):16-23.

ZONG Fukai.Frequency dispersion characteristic and discretization of the finite element analysis in wave propagation problems[J].Explosion and Shock Waves,1984,3(4):16-23.

[5] 譚路.基于波數(shù)譜的細(xì)長(zhǎng)殼體圓柱殼振動(dòng)與聲輻射噪聲特性研究[D].武漢：海軍工程大學(xué)，2013.

[6] 田麗思,李相輝,馬越峰,等.MSC Nastran 動(dòng)力分析指南[M].北京：中國(guó)水利水電出版社,2012.

[7] KIM J K, IH J G.Prediction of sound level at high-frequency bands by means of a simplified boundary element method[J].Journal of the Acoustical Society of America, 2002,112(6): 2645-2655.

Principle to select element size of FE-BE algorithm for a ring stiffened cylindrical shell

HUANG Zhenwei, ZHOU Qidou, JI Gang, MENG Qingchang, DUAN Jiaxi

(Department of Ship and Engineering, Naval Univ.of Engineering, Wuhan 430033, China)

Investigating the principle to select element size of structural wet surface in the FE-BE algorithm has a great significance.The structural bending wavelength of the main vibration mode of a finite length ring stiffened cylinder was computed with the wave number spectrum method to provide some principle for selection of element size of a structural wet surface.Numerical computation showed the effectiveness of the principle.

ring stiffened cylindrical shell; FE-BE algorithm; vibration and sound radiation; element size

國(guó)家自然科學(xué)基金(51309230)

2015-10-27 修改稿收到日期：2016-02-19

黃振衛(wèi) 男，博士生，1986年生

周其斗 男，教授，博士生導(dǎo)師，1962年生

TH1

A

10.13465/j.cnki.jvs.2017.05.023