基于解析解和邊界元解的圓柱殼聲輻射對比研究

張俊杰，原春暉，劉 彥，吳武輝

(中國艦船研究設計中心 船舶振動噪聲重點實驗室，武漢 430064)

圓柱殼結構在工程中廣泛應用，其產生的振動和聲輻射也越來越受到關注。圓柱殼結構的分析方法分為解析法和數值法。由于解析法僅能研究結構形式簡單的圓柱殼，對于工程中復雜的圓柱殼，一般采用數值解法。在中低頻頻段，數值解法一般是采用有限元和邊界元聯合的方式。其計算流程一般是:① 建立復雜殼體結構有限元模型，施加邊界條件并計算;② 提取殼體表面網格作為邊界元網格;③ 提取殼體外表面的法向位移;④ 將邊界元網格和法向位移導入到邊界元軟件中計算，得到相關的聲學參量。詳細計算流程見文獻［1］。上述求解過程清晰，簡單易行，國內外學者運用有限元和邊界元聯合分析了圓柱殼結構的振動和聲輻射。

商德江等［2］運用有限元和邊界元分析了加肋雙層圓柱殼在受點力激勵作用的聲輻射特性，并運用有限元和解析解比較了彈性球殼在軸對稱點力激勵下表面法向位移和表面聲壓分布。徐張明等［3］運用有限元和邊界元分析了帶浮筏裝置的水中雙層加了圓柱殼振動和噪聲，文中對長方體聲腔和橢圓球殼簡單結構運用編程和邊界元軟件計算進行了聲振對比，而對于帶浮筏裝置復雜圓柱殼采用有限元軟件和邊界元軟件聯合處理方式。賀晨等［4］運用有限元、邊界元和統計能量分析方法對圓柱殼在流場中振動和聲輻射效率進行了計算，文中運用邊界元和解析解比較了空氣中圓柱殼的輻射效率，兩者得到的結果在中低頻率相差較大。其他作者如石煥文［5］、曹貽鵬［6］、陳海龍［7］和姚熊亮［8］也用有限元和邊界元對單層圓柱殼，雙層圓柱殼等殼體結構的振動和聲輻射進行了研究。

在上述文獻中，邊界元計算一般是假設每個波長六個單元線性單元或三個二次單元［9］以確定最大可以計算頻率，并沒有考慮不同加載方式—基于單元(Element)加載和基于節點(Node)加載對計算結果的影響，以及這兩種方式網格加密后收斂的問題。事實上，這兩種加載方式對計算結果以及網格加密的收斂性問題有較大影響。基于此，本文針對工程廣泛應用的圓柱殼結構，通過與解析法對比，探討了邊界元兩種加載方式對聲輻射影響，以及網格加密的收斂性問題，得到一些對工程有益的結論。

1 模型簡介和理論計算

1.1 模型簡介

流場中簡支在無限障板上有限長圓柱殼的幾何、材料參數和坐標系如圖1所示。圓柱殼半徑為R，長度為L;流場中聲速和密度分別為c0和ρ0;殼體楊氏模量、泊松比和密度分別為E1、v1和ρ1。分別用r，θ和 x表示柱坐標系的徑向、周向和軸向，w，v和u表示這三個方向的位移分量。

圖1 有限長圓柱殼幾何、材料參數和坐標系Fig.1 Geometry material parameters and coordinate system of finite cylindrical shell

1.2 圓柱殼運動方程

對于理想流體中受法向激勵力作用的單層圓柱殼，為了便于求出殼體的法向位移，給出用阻抗表示的法向位移的關系式為［10－11］:

1.3 流場聲壓推導

理想流體中，流場中的聲壓滿足Helmholtz方程，殼體和流場邊界處滿足位移的連續方程。運用傅立葉變換，最后得到時域中聲壓表達式:

若不考慮軸向模態耦合，把式(5)代入式(3)，運用三角函數的正交性得到(r=R)

為方便計算，無量綱化。設λ*=kR，k*=k0R，則式(6)可以變為:

為了加快計算速度，可以對式(7)進行變量替換，具體可以參考文獻［11］。

1.4 作用力

考慮點激勵作用，且激勵力作用在殼體內表面，沿著徑向。點激勵運用式(3)模態展開后的系數

式中:n=0時εn=1;n≥1時εn=2，F0外力幅值，x0和θ0分別是外力在軸向和周向的作用位置。

1.5 耦合方程求解

把式(6)和式(8)代入式(1)，便可以求出對應于模態(m，n)的對稱或者反對稱法向位移，接著根據式(2)的算出殼體的位移w，以及根據式(6)和式(3)算出殼體和聲場耦合面的聲壓p。

1.6 輻射聲功率

一旦求出殼體的位移w和殼體與聲場耦合面的聲壓p，輻射聲功率便可以求出:

輻射聲功率級定義為:

式中:W0=10－12W。

2 邊界元計算

2.1 邊界元基本方程

根據邊界積分方程，把輻射面離散為小單元，并設離散后節點數目為N個，則節點處的聲壓ps與法向振速向量ws的關系式可以表示為:

式中，矩陣［A］和［B］為影響矩陣，與振動面形狀、流場屬性和頻率相關，與法向速度分布無關。

通過矩陣變換，式(11)可得到:

根據式(12)，一旦給定節點速度便可計算得到節點處聲壓，進一步可以得到輻射聲功率。

2.2 輻射面的離散形式

對于工程復雜圓柱殼來說，聲振解析解無法求出，需要借助數值解法。數值解法的一般求解過程是先在有限元軟件中計算輻射面的振動速度，再將表面速度導入邊界元中插值進行聲場計算。有限元軟件可選擇性較多，如Ansys，MSC.Patran等，而聲場計算一般是運用軟件Sysnoise進行計算。Sysnoise軟件提供了基于單元(Element)和基于節點(Node)的兩種加載模式。

對于圖1的殼體結構，圓柱面展開示意圖如圖2(a)。根據基于單元(Element)和基于節點(Node)的兩種加載模式，以及考慮離散后節點數目大致相等的情況，圖2(b－d)給出了典型三種分區形式。圖中粗線表示單元的邊界，小黑點表示邊界元施加速度的位置。圖2(b)是基于單元(Element)的加載形式，速度邊界條件加載在單元中心位置上，單元速度分布均勻，等效為活塞運動。圖2(c)－圖2(d)是基于節點(Node)的加載形式，其中圖2(c)的單元與圖2(b)的單元位置和大小相同，只是速度邊界條件加載在單元的節點上;圖2(d)中單元節點的位置與圖2(b)中施加速度位置(除去左右兩側)相同。

圖2 殼體展開圖及網格分區圖Fig.2 Unfolded drawing and different meshes of cylindrical shell

3 數值計算及討論

3.1 殼體幾何、材料參數以及數值計算流程

3.1.1 殼體幾何和材料參數

計算中殼體幾何參數R=1 m，長度為L=3 m;流場中聲速和密度分別為c0=1 500 m/s和ρ0=1 000 kg/m3;殼體楊氏模量E1和泊松比v1以縱波cd1、橫波cs1和密度ρ1表示，其中縱波 cd1=cd(1+iζd)1/2，橫波 cs1=cs(1+iζs)1/2，cd=5 800 m/s，ξd=0.001，cs=3 100 m/s，ξs=0.01 以及ρ1=7 850 kg/m3，激勵力 F0=1 N。

3.1.2 解析解和邊界元計算流程

邊界元計算時殼體的表面位移從解析解中提取出來，解析解和邊界元對比計算流程為:

(1)運用式(1)求解出有限長圓柱殼在激勵力作用下的法向位移 Wαmn，運用式(9)求出殼體輻射聲功率;

(2)從解析解的位移中提取出對應于圖2(b－d)小黑點處的位移值;

(3)對應于圖2(b－d)，生成邊界元輸入文件并計算得到輻射聲功率;

(4)對比第(1)步和第(3)步的計算結果。

3.2 基于單元和基于節點的兩種加載方式結果與解析解對比分析

按照工程經驗，邊界元單元大小要滿足每波長六單元的基本要求。對于圖2(b－d)的分區形式，除非特殊說明，其邊界元網格的分區大小分別為20×40分區、20×40分區和21×40分區。按照這個要求計算出來的頻率是1 593 Hz。在下面的數值計算中，討論的頻率范圍為［0，1 000 Hz］，因此，按照每波長六單元的的要求結果是收斂的。下面分別就單個周向模態激勵和所有周向模態激勵(點激勵)時解析解和邊界元計算結果進行對比討論。

3.2.1 單個周向模態激勵情形

首先比較了點激勵的分量—某個周向模態激勵的情形，此時只需令式(8)中某個單個周向模態激勵有值，而其它周向模態激勵為零，且在計算時式(1)的軸向模態m截斷上限取為30。計算中點力作用位置為x0=3L/5和θ0=π/6，取周向模態n=1和n=2兩種情況。其中邊界元計算考慮了圖2(b－d)分區和加載的情況。解析解和邊界元的結果分別如圖3－4。由于三種分網方式得到的結果畫在一幅圖上不清晰，這里在保持坐標軸不變的情況下分別與解析解對比給出。需要說明的是，圖中曲線指向橫軸的陰影區域表示這個區間段算出的聲功率為負值，這與實際情況不符，稱該陰影區域為該方法“不可計算頻率范圍”，該“不可計算頻率范圍”是由于邊界元計算中奇異性引起的。

單一周向模態激勵時，殼體表面的速度變化平緩，如圖5是周向模態n=1頻率為300 Hz時殼體表面速度實部沿軸向和周向分布圖。此時在周向模態n=1和n=2，運用解析解和三種分區方式計算的輻射聲功率如圖3和圖4所示。由圖可知，三種分區計算結果除了“不可計算頻率范圍”，其它頻率段與解析解的結果趨勢保持一致。通過仔細比較曲線間的相對誤差可知，圖2(b)分區基于單元(Element)加載的結果要優于基于節點(Node)加載的結果。事實上，這不是個例，接下來討論點激勵時也是這個規律。

圖3 周向模態n=1時解析和不同分區結果比較圖Fig.3 The comparisons of analytical solution and results of different meshes with n=1

圖4 周向模態n=2時解析和不同分區結果比較圖Fig.4 The comparisons of analytical solution and results of different meshes with n=2

圖5 周向模態n=1時速度實部沿軸向和周向分布Fig.5 The distribution of real part of velocity in axial and circumferential with n=1

3.2.2 所有周向模態激勵(點激勵)情形

點激勵展開式式(3)包含所有周向和周向模態疊加，殼體表面的速度變化較大，如圖6是單點激勵情況下300 Hz時殼體表面速度實部沿軸向和周向分布圖。

點激勵計算時式(1)的軸向和周向模態m和n截斷上限均取為30。不失一般性，給出了單點激勵和兩點同時激勵時殼體輻射的聲功率如圖7和圖8。單點激勵位置為x0=3L/5和θ0=π/6;兩個點力同時作用時激勵位置分別為x0=3L/7和θ0=π/4，x0=3L/5和θ0=π/6。由圖7－8可知，無論是從三種分區結果與解析解的相對誤差，還是從“不可計算頻率范圍”的個數，圖2(b)基于單元(Element)加載的結果明顯要好于圖2(c－d)基于節點(Node)加載的結果。對比圖3－8可以推斷，殼體表面速度分布越不均勻，圖2(c－d)基于節點(Node)加載的結果可能與解析值差別越大，而圖2(b)基于單元(Element)加載的結果與解析值接近程度與殼體表面速度分布基本無關。

圖6 點力作用下速度實部沿軸向和周向分布Fig.6 The distribution of imaginary part of velocity in axial and circumferential with a point force excitation

圖7 一個點力作用下解析和不同分區結果比較圖Fig.7 The comparisons of analytical solution and results of different meshes with a point force excitation

圖8 兩個點力作用下解析和不同分區結果比較圖Fig.8 The comparisons of analytical solution and results of different meshes with two point force excitation

3.3 點激勵下基于單元和基于節點的兩種加載方式網格加密收斂性討論

加密收斂性是討論在已經滿足了每波長六單元所要求的頻率以內，分析圖2(b－d)分區網格繼續加密結果收斂性問題，也就是與解析值逼近的問題。選取的分區數目分別為16×32、20×40和30×60，三種分區方式每波長六單元對應的頻率分別為1 275 Hz、1 593 Hz和2 388 Hz。計算中考慮的最高頻率為1 000 Hz，是滿足每波長六單元要求的。根據圖7－8結論可知，基于節點加載的分區方式圖2(c)和圖2(d)的結果類似，因此這里僅選取圖2(b－d)具有代表性的基于單元(Element)加載的圖2(b)和基于節點(Node)加載的圖2(c)。事實上，這兩幅圖的網格是一致的，只是加載方式不同。

圖9 圖2(b)分區網格加密收斂性Fig.9 The convergence of Fig.2(b)when refining the mesh

圖10 圖2(c)分區網格加密收斂性Fig.10 The convergence of Fig.2(c)when refining the mesh

圖9和圖10分別給出了圖2(b)和圖2(c)分區網格加密收斂性。由圖9可知，隨著網格細分，圖2(b)基于單元(Element)加載的聲功率與解析解越來越逼近，在網格30×60時幾乎重合，這說明基于單元(Element)加載的方法在網格加密后具有很好的收斂性。而對于圖10，在基于節點(Node)加載在網格加密后，在頻段［0－150 Hz］是收斂的。在大于150 Hz的頻段，不同分區數目出現了不同的“不可計算頻率范圍”。根據這個性質，可以試著劃分不同的網格數量以避開“不可計算頻率范圍”，但是總的來說網格加密后不一定得到收斂結果。

4 結論

本文針對工程廣泛應用的圓柱殼結構，在單元數目大致相等的情況下，通過與解析法的結果進行對比，探討邊界元了基于單元(Element)加載方式和基于節點(Node)加載方式對聲輻射影響，以及這種加載方式網格加密的收斂性問題，得到如下結論:

(1)當殼體速度變化不大時，比如單個模態激勵情形，基于單元(Element)加載方式和基于節點(Node)加載方式的計算結果與解析解計算結果基本一致，基于單元(Element)加載方式的結果要優于基于節點(Node)加載結果。

(2)當殼體速度變化較大時，比如單點激勵或多點激勵情形，基于單元(Element)加載方式的計算結果與解析解計算結果基本一致，而基于節點(Node)加載方式的結果與解析結果在部分頻段相差較大，而且可能出現多個“不可計算頻率范圍”。

(3)殼體表面速度分布越不均勻，基于節點(Node)加載的結果可能與解析值差別越大，而基于單元(Element)加載的結果與解析值接近程度與殼體表面速度分布無關。

(4)當改變加密網格大小確定收斂性時，基于單元(Element)加載方式在網格加密后具有很好的收斂性;而基于節點(Node)加載方式并不一定得到收斂結果。

［1］王 晶，商德江.Ansys和Sysnoise之間的數據接口技術研究［J］.應用科技，2004，31(8):32－34.

WANG Jing，SHANG De-jiang.Study of the interface between ANYSYS and SYSNOIS［J］. Applied Science and Technology，2004，31(8):32 －34.

［2］商德江，何祚鏞.加肋雙層圓柱殼振動聲輻射數值計算分析［J］.聲學學報，2001，26(3):193－201.

SHANG De-jiang，HE Zuo-yong.The numerical analysis of sound and vibration from a ring－stiffened cylindrical double－ shell by FEM and BEM［J］.Acta Acustica，2001，26(3):193－201.

［3］徐張明.帶浮筏裝置的水中殼體模型振動、聲輻射及其結構參數靈敏度研究［D］.上海:上海交通大學，2003.

［4］賀 晨，盛美萍，石煥文，等.圓柱殼體振動聲輻射效率數值計算分析［J］.噪聲與振動控制，2006，26(4):51－54.

HE Chen，SHENG Mei-ping，Shi Huanwen，et al.The numerical analysis of vibration and sound radiation efficiency from a cylindrical shell［J］.Noise and Vibration Control，2006，26(4):51 －54.

［5］石煥文，盛美萍，孫進才，等.加縱肋平底圓柱殼振動和聲輻射的FEM/BEM研究［J］.振動與沖擊，2006，25(2):88－92.

SHI Huan-wen，SHENG Mei-ping，SUN Jin-cai，et al.On FEM/BEM forthe problemsofvibration and acoustic radiation from an axially stiffened cylindrical shell with two end plates［J］.Journal of Vibration and Shock，2006，25(2):88－92.

［6］曹貽鵬，張文平.軸系縱振對雙層圓柱殼體水下聲輻射的影響研究［J］.船舶力學，2007，11(2):293－299.

CAO Yi-peng;ZHANG Wen-ping.A study on the effects of the longitudinal vibration of shafting on acoustic radiation from underwater double cylindrical shell［J］.Journal of Ship Mechanics，2007，11(2):293 －299.

［7］陳海龍，錢德進，計 方，等.靜水壓與阻尼覆蓋層對圓柱殼聲輻射的影響［J］.傳感器與微系統，2009，28(3):42－44.

CHEN Hai-long，QIAN De-jin，JI Fang，et al..Effects of hydrostatic pressure and damping materialon acoustic radiation of cylindrical shell［J］.Transducer and Microsystem Technologies，2009，28(3):42 －44.

［8］姚熊亮，楊娜娜，陶景橋.雙層殼體水下振動和聲輻射的仿真分析［J］.哈爾濱工程大學學報，2004，25(2):136－140.

YAO Xiong-liang，YANG Nana，TAO Jing-qiao.Numerical research on vibration and sound radiation of underwater double cylindrical shell［J］.Journal of Harbin Engineering University，2004，25(2):136 －140.

［9］李增剛.Sysnoise Rev 5.6詳解［M］.北京:國防工業出版社，2005.

［10］何祚鏞.結構振動與聲輻射［M］.哈爾濱:哈爾濱工程大學出版社，2001.

［11］張俊杰.基于不同理論的流場中圓柱殼振動能量流和輻射聲功率研究［D］.武漢:華中科技大學，2010.