有限長圓柱殼-流場部分耦合系統(tǒng)聲振建模及機理研究

郭文杰, 楊 舟

(華東交通大學(xué) 交通運輸工程學(xué)院， 南昌 330013)

圓柱殼-流場耦合系統(tǒng)在船舶與海洋工程、油氣儲運工程、航空航天工程等諸多領(lǐng)域有廣泛的應(yīng)用，如潛艇、輸油管路、含液體燃料的火箭貯箱等。早期的研究主要針對圓柱殼內(nèi)部充滿液體[1-3]或者浸沒于無限域[4-6]這類理想工況，由于結(jié)構(gòu)表面與流場完全耦合且不考慮自由液面聲邊界的影響，數(shù)學(xué)處理難度較低，可直接將聲壓在結(jié)構(gòu)位移場坐標(biāo)系展開，再根據(jù)流固耦合面速度連續(xù)條件得到聲場與位移場的數(shù)學(xué)聯(lián)系，已形成了一套非常清晰的研究思路。

實際的工程問題中，圓柱殼與流場兩種介質(zhì)也常處于部分耦合的狀態(tài)(圓柱殼軸線與自由液面平行)，比如海面系泊狀態(tài)(浮態(tài))下的潛艇、油氣混輸管道、未盛滿液體的圓柱型儲液容器等，其動力學(xué)特性關(guān)系到結(jié)構(gòu)的安全性、功能性，因此開展相關(guān)理論研究工作尤為必要。

(a)

(b)

(a)

(b)

為了克服上述不足，首先就是需要突破慣性思維的桎梏。針對部分耦合問題的建模難點，郭文杰等[19-20]在實踐中探索出利用聲場、結(jié)構(gòu)兩套不同坐標(biāo)系建模的方法，并成功對無限長圓柱殼-外流場部分耦合聲振問題進行建模求解。具體來講，結(jié)構(gòu)位移場建立在以殼體圓心為坐標(biāo)原點的柱坐標(biāo)系下，便于描述殼體運動，聲場坐標(biāo)系的原點建立在自由液面上，便于利用正弦函數(shù)得到聲壓函數(shù)表達式，接著利用Galerkin法處理聲固耦合面上振速連續(xù)條件，得到聲場和位移場系數(shù)向量的矩陣關(guān)系，進而可對耦合系統(tǒng)振動特性進行理論求解?；趦商鬃鴺?biāo)系建模方法，Zhao等[21]研究了簡諧點荷載作用下無限長部分浸沒圓柱殼的聲輻射特性，并開展試驗驗證。

在前期基礎(chǔ)上，本文進一步將圓柱殼模型由無限長拓展到有限長，即考慮軸向波數(shù)對聲固耦合系統(tǒng)的影響。此外，為深入探究內(nèi)、外聲場與圓柱殼結(jié)構(gòu)耦合振動問題的異同，本文對相同液面高度下部分充液與部分浸沒工況開展辨析研究，并從數(shù)學(xué)機理上揭示其異同的成因。

1 理論推導(dǎo)

1.1 模型描述

有限長圓柱殼參數(shù)定義如下：殼長為L，殼厚為h，橫截面中面半徑為Rs，結(jié)構(gòu)密度為ρ，楊氏模量為E，泊松比為μ，流體密度為ρf。殼體坐標(biāo)系原點O位于左端面圓心，相應(yīng)的直角坐標(biāo)(x,y,z)如圖3所示。實際理論推導(dǎo)中殼體坐標(biāo)系采用柱坐標(biāo)系(x,r,φ)，其中：x，r，φ分別表征軸向、徑向、周向角(起始點為y軸正方向，逆時針為正)。聲場坐標(biāo)系原點定義為y軸與自由液面交點Q，同樣采用柱坐標(biāo)系(x,R,θ)，R表征徑向，θ表征周向角，逆時針為正。定義Q點于y軸坐標(biāo)值的負數(shù)為液面特征高度H(位于y軸負半軸取正值)。α為液面對應(yīng)的夾角，滿足H/Rs=sinα。結(jié)構(gòu)坐標(biāo)系與聲學(xué)坐標(biāo)系夾角定義為β，當(dāng)自由液面在殼體圓心上方時β取值為正值。簡諧激勵力位于結(jié)構(gòu)坐標(biāo)系(Rs,x0,φ0)處，幅值為F0。

(a)

(b)

1.2 結(jié)構(gòu)運動方程

圓柱殼運動方程建立在結(jié)構(gòu)坐標(biāo)系(x,r,φ)下，采用Flügge殼體理論描述

(1)

式中：u，v和w分別為軸向、周向和徑向的位移；fr為外力載荷；fp為圓柱殼內(nèi)表面的流體負載，[L]為經(jīng)典的Flügge殼體理論微分算子，詳見文獻[22]。

為了研究的方便，本文假設(shè)結(jié)構(gòu)兩端的邊界條件均為簡支，由此位移場及載荷可展開為三角級數(shù)形式(時間項略去)

(2)

式中：Umn，Vmn和Wmn分別為3個方向位移的幅值；fmn為流體負載fp的幅值；Fmn為外力載荷fr的幅值；m為軸向半波數(shù)；n為周向函數(shù)展開序列數(shù)；ω=2πf為角頻率，f為頻率；軸向波數(shù)km=mπ/L。

基于函數(shù)的正交性，可將式(2)代入式(1)，并進行正交化處理，即可得到如下的控制方程

(3)

式中，矩陣[T]中的具體算子詳見Guo等的研究。

根據(jù)式(3)可進一步化簡，消去Umn，Vmn后可得到形式更為簡潔的振動控制方程

(4)

式中，Imn=(T11T22-T12T21)/det(T)。

點激勵力可以表示為如下的函數(shù)形式，國際單位為N

(5)

式中，δ(·)為狄利克雷函數(shù)。

同樣對式(5)進行正交化處理，可以得到激勵力幅值Fmn的表達式

(6)

1.3 聲場建模

聲場建模是理論研究的一大重點，而聲邊界的處理是聲場建模的主要難點。要對部分耦合系統(tǒng)進行理論研究，就必須建立滿足所有聲邊界條件的數(shù)學(xué)模型。首先，聲壓p必須滿足聲學(xué)Helmholtz方程

(7)

式中：kf=ω/cf為聲波波數(shù)，ω=2πf為角頻率，f為頻率；?2為拉普拉斯算子。

然后，聲壓p還要進一步滿足無窮遠處Sommerfeld輻射條件

(8)

最后，聲壓p還需于自由液面處滿足聲壓釋放條件(空氣密度相比水的密度極小，且聲速也相對較小，故而空氣的特性阻抗遠小于水，因此水中聲波入射到液面時可作為絕對軟邊界處理，即認為聲壓為零

p=0, 自由液面處

(9)

本文將聲壓函數(shù)建立在獨立的聲場坐標(biāo)系(x,R,θ)下，然后通過采用正弦三角級數(shù)來自動滿足液面處聲壓釋放條件，具體形式如下

(10)

由于聲壓負載在結(jié)構(gòu)表面部分存在，可以寫成分段函數(shù)的形式

(11)

根據(jù)式(2)中載荷的展開形式，基于正交化處理可以得到聲壓負載fp幅值fmn的表達式

(12)

1.4 控制方程求解

根據(jù)控制方程式(4)可知，在確定聲壓負載幅值fmn的表達式后，控制方程求解的關(guān)鍵在于構(gòu)建聲壓幅值和位移幅值的數(shù)學(xué)聯(lián)系。

根據(jù)聲固耦合面處位移連續(xù)條件，聲壓與徑向位移有如下關(guān)系

(13)

由于式(13)難以直接解析求解，可選擇合適的數(shù)值方法進行求解。由于Galerkin法具有計算精度較高、適用性較廣的優(yōu)點，本文選擇該方法對式(13)進行數(shù)值求解。構(gòu)成位移場和聲場的基函數(shù)主要有兩類，對應(yīng)的權(quán)函數(shù)也有兩類，一類基于位移周向展開函數(shù)；另一類基于聲壓周向展開函數(shù)

(14)

因此可將方程式(13)轉(zhuǎn)化為二重積分的弱形式，且積分域為聲固耦合面

(15)

式中：n取值為-N～N；j取值為1～2N+1，N為截斷項數(shù)。

ρfω2[Vs]{Wmn}=[Vp]{Amj}

(16)

式中：[Vs]和[Vp]均為2N+1階方陣，{Wmn}=[Wm,-N,Wm,-N+1, …,Wm,N-1,Wm,N]T為位移幅值向量；{Amj}=[Am,1,Am,2, …,Am,2N,Am,2N+1]T為聲壓幅值向量；上標(biāo)T為對矩陣進行轉(zhuǎn)置。

由于聲固耦合面上(r,φ)坐標(biāo)系和(R,θ)坐標(biāo)系之間存在幾何對應(yīng)關(guān)系，根據(jù)式(15)，[Vs]，[Vp]中的元素可表示為如下積分形式(假設(shè)權(quán)函數(shù)基于徑向位移的周向展開函數(shù))

(17)

exp[i(a-1-N)φ]dφ

(18)

式中，β=3π/2-θ-φ,R和θ的取值滿足如下幾何關(guān)系

(19)

同時，式(12)中{fmn}也可表示為矩陣的形式

{fmn}=[Tp]{Amj}

(20)

式中，矩陣[Tp] 中每一個元素的具體表達式如下

N)φ]Kb(krR)sin(bθ)dφ

(21)

最終可以得到僅與徑向位移幅值相關(guān)的耦合系統(tǒng)控制方程

(22)

式中：矩陣[J]為2N+1階單位矩陣；[G]矩陣為對角矩陣；[G]j, j=Im, j-1-N；{Fmn}=[Fm,-N,Fm,-N+1, …,Fm,N-1,Fm,N]T為激勵力向量，詳見式(6)。直接求解控制方程式(22)可得到位移幅值向量{Wmn}，進而可得到結(jié)構(gòu)任意點振動位移響應(yīng)。

當(dāng)求解自由振動時，激勵力為零，式(22)轉(zhuǎn)化為特征值求解問題

{0}

(23)

式中，{0}為零向量，具體求解時可通過定義式(23)中系數(shù)矩陣行列式值為零來求解角頻率ω，再根據(jù)f=ω/(2π)，可以得到各階固有頻率。

2 數(shù)值分析

本文參數(shù)選擇如下：結(jié)構(gòu)密度ρ=7 850 kg/m3，彈性模量E=206 GPa，泊松比μ=0.3，殼體長度L=1.284 m，中面半徑R=0.18 m，殼厚h=0.003 m，流體密度ρf=1 025 kg/m3，流體聲速cf=1 500 m/s。當(dāng)計算受迫振動響應(yīng)時結(jié)構(gòu)阻尼因子η取為0.01，復(fù)楊氏模量為E′=E(1+iη)。

2.1 方法收斂性分析

由圖4可以看出，隨著截斷項數(shù)M，N的增加，徑向均方振速級VML很快趨于穩(wěn)定；當(dāng)激勵頻率在1 200 Hz以內(nèi)時，M和N均取12時結(jié)果已足夠收斂。后文計算中M，N取值均為12。

(a)

2.2 方法適用范圍

接著進一步分析本文方法特征高度的取值范圍，定義無量綱特征高度為H/Rs，取值為-0.8～0.8，間隔0.1。分別計算兩種不同權(quán)函數(shù)下首階固有頻率值，并與有限元軟件Nastran中仿真(虛擬質(zhì)量法)計算結(jié)果進行對比，如表1所示。

由表1可知，權(quán)函數(shù)無論是位移基函數(shù)還是聲壓基函數(shù)，在保證計算精度較高的情況下，無量綱特征高度H/Rs取值范圍為-0.8～0.8。由Amabili的研究可知，Amabili提出的斜邊近似法適用范圍約為-0.37～0.37，換句話說，本文方法適用范圍相比斜邊近似法要大得多。

2.3 方法準(zhǔn)確性驗證

接下來，本文開展了自由振動和受迫振動的方法準(zhǔn)確性驗證，且權(quán)函數(shù)選擇位移基函數(shù)，有限元仿真對比繼續(xù)采用虛擬質(zhì)量法。

2.3.1 自由振動準(zhǔn)確性分析

首先針對自由振動，本文分別取浸沒深度H=-0.09 m，H=0和H=0.09 m，計算前10階固有頻率，并與有限元軟件Nastran仿真計算結(jié)果進行對比，如表2所示。定義固有頻率的相對誤差Error=(f1-f2)f2×100%。

表1 各無量綱特征高度下首階固有頻率對比

表2 當(dāng)H=-0.09 m時部分浸沒圓柱殼前10階固有頻率對比

由表2～表4可知，對于部分浸沒圓柱殼，本文方法計算得到的前10階固有頻率與Nastran仿真計算結(jié)果十分吻合，最大相對誤差的絕對值不超過1%，這也說明本文方法計算自由振動是準(zhǔn)確可靠的。

2.3.2 受迫振動準(zhǔn)確性分析

緊接著，本文進一步分析受迫振動時方法準(zhǔn)確性。假設(shè)點激勵力位于坐標(biāo)x0=L/2，φ0=0處，激勵力幅值F0=1 N，激勵頻率取值范圍為2～500 Hz，間隔2 Hz。觀測點位于坐標(biāo)x1=L/4，φ1=0處。分別計算特征深度H=-0.09 m和H=0.09 m時徑向位移響應(yīng)級頻譜曲線，并與有限元軟件Nastran仿真(虛擬質(zhì)量法)結(jié)果進行對比，如圖5所示。本文定義徑向位移響應(yīng)級[23]RDL=20×lg(w/w0)，其中w為測點徑向位移絕對值，w0=1×10-12m。

表3 當(dāng)H=0時部分浸沒圓柱殼前10階固有頻率對比

表4 當(dāng)H=0.09 m時部分浸沒圓柱殼前10階固有頻率對比

從圖5可以看出，本文方法頻響曲線和Nastran仿真結(jié)果吻合良好，進一步驗證了本文方法的準(zhǔn)確性。此外，本文方法計算效率遠高于仿真方法(虛擬質(zhì)量法)計算效率。以圖5(b)中頻響曲線計算為例，當(dāng)計算機配置均為4 GHz CPU與24.0 GB RAM時，虛擬質(zhì)量法在2 400個單元的仿真計算規(guī)模下耗時約2.2 h。而基于本文方法采用MATLAB2014a編程計算約耗時4 min，計算效率提高數(shù)十倍。

2.4 液面特征高度H對固有頻率的影響

液面特征高度的變化會引起附連水質(zhì)量的變化，進而影響耦合系統(tǒng)自振特性。為研究自振特性隨特征高度的變化規(guī)律，取無量綱特征高度H/Rs為-0.8～0.8，計算前4階固有頻率隨特征高度變化規(guī)律，如圖6所示。

圖5 不同方法下部分浸沒圓柱殼徑向位移級對比

(a) 第一階

從圖6可以發(fā)現(xiàn)，同階次固有頻率隨著特征高度增大而減小，主要原因是由于特征高度增大導(dǎo)致濕表面積增加，因而使得附連水質(zhì)量增加，從而增加了耦合系統(tǒng)的等效質(zhì)量，導(dǎo)致同階次固有頻率降低。此外，同階次固有頻率隨特征高度增加而下降的速率無顯著規(guī)律，這是因為自由液面對聲波的反射也會對附連水質(zhì)量產(chǎn)生影響。因此在兩類因素的綜合影響下，雖然同階次固有頻率整體上隨液面增高而降低，但是下降速率無明顯規(guī)律，忽快忽慢。

3 部分充液工況

3.1 部分充液工況建模及求解特點

部分充液工況與部分浸沒工況在理論建模及數(shù)學(xué)推導(dǎo)方面是基本相似的，如圖7所示。但由于流場由外變內(nèi)，主要存在兩點差異：一是聲學(xué)Helmholtz方程的基本解有區(qū)別；二是聲固耦合面振速連續(xù)條件有方向差異。

3.2 部分充液工況下本文方法的準(zhǔn)確性驗證

為驗證部分充液工況下本文方法的準(zhǔn)確性，本節(jié)繼續(xù)采用有限元軟件Nastran的仿真計算結(jié)果作為對比。特征高度分別取H=-0.127 3 m，0，0.127 3 m，兩種方法的前10階固有頻率對比如表5～表7所示。

從表5～表7可以看出，部分充液工況下本文方法計算得到的前10階固有頻率與Nastran仿真計算結(jié)果十分吻合，最大相對誤差的絕對值不超過1%，這也說明本文方法計算部分充液圓柱殼振動特性是準(zhǔn)確可靠的。

(a)

(b)

表5 當(dāng)H=-0.127 3 m時部分充液圓柱殼前10階固有頻率對比

表6 當(dāng)H=0時部分充液圓柱殼前10階固有頻率對比

表7 當(dāng)H=0.127 3 m時部分充液圓柱殼前10階固有頻率對比

為進一步說明本文方法的準(zhǔn)確性，將本文方法計算的固有頻率值與Amabili等的研究試驗結(jié)果進行對比，如表8所示。(Amabili等研究模型參數(shù)如下：L=0.664 m，R=0.175 m，h=0.001 m，ρ=7 680 kg/m3，E=206 GPa，μ=0.3，ρf=1 000 kg/m3。)

表8 不同充液高度下本文方法與試驗結(jié)果頻率對比

從表8可以看出，不同充液高度下解析法和試驗方法得到的數(shù)據(jù)大體吻合良好，進一步說明本文方法的準(zhǔn)確性。

4 部分充液與部分浸沒工況對比分析

4.1 全充液與全浸沒工況的對比

對于全充液圓柱殼或者無限域圓柱殼，當(dāng)所研究的頻率較低時，其同階次的固有頻率往往是比較接近的。主要原因是由于頻率較低時對應(yīng)的無量綱徑向波數(shù)krRs比較小，根據(jù)貝塞爾函數(shù)小宗量展開的性質(zhì)，內(nèi)、外流場對應(yīng)的流體負載在數(shù)學(xué)表達式上可認為近似相等[24]。

(24)

式中：m，n分別為軸向半波數(shù)和周向波數(shù)；Cmn為與流體相對位置(內(nèi)流場或外流場)無關(guān)的量。此外式(24)中貝塞爾函數(shù)也可能取其他形式(比如第一類漢克爾函數(shù)和第一類貝塞爾函數(shù))，此時式(24)依然可簡化為約等號右側(cè)形式。

本處簡單說明一下式(24)中約等號右側(cè)表達式是如何近似得到的。以內(nèi)流場為例，當(dāng)無量綱波數(shù)較小時，可對貝塞爾函數(shù)進行小宗量近似展開

(25)

式中：Γ(· )為伽馬函數(shù)，且Γ(n+1)=n！。

根據(jù)式(25)可以得到如下表達式

(26)

據(jù)此可判斷全充液或全浸沒(無限域)時，二者對應(yīng)的前幾階固有頻率值應(yīng)該十分接近。比如以第一階固有頻率為例，全充液時為97.9 Hz，全浸沒時為98.9 Hz，差異很小。

4.2 部分充液與部分浸沒工況的對比

為了進一步了解部分耦合時內(nèi)、外流場是否也有類似的性質(zhì)，同時也為分析部分充液時固有頻率隨特征高度的變化規(guī)律，將無量綱特征高度H/Rs取值為-0.8～0.8，計算對應(yīng)的第一、第三階固有頻率值，并與圖6(部分浸沒工況)中部分數(shù)據(jù)進行對比，結(jié)果如圖8所示。

(a) 第一階

(b) 第三階

從圖8可以看出，部分充液工況下固有頻率隨特征高度的變化規(guī)律與部分浸沒時類似，即同階次固有頻率隨液面的增高而降低且下降的速率無明顯規(guī)律。此外，值得注意的是，半充液工況與半浸沒工況下(H/Rs=0)同階次固有頻率值幾乎相同。為了進一步證實這一規(guī)律，本文將半充液工況及半浸沒工況下前10階固有頻率繪圖比對，如圖9所示。

從圖9可以明顯看出，無論是本文方法還是數(shù)值仿真，半浸沒工況及半充液工況下前10階固有頻率值是幾乎一致的。

4.3 數(shù)學(xué)機理分析

為了從數(shù)學(xué)機理上揭示和解釋圖9中規(guī)律的成因，可依據(jù)本文理論推導(dǎo)及借鑒貝塞爾函數(shù)的小宗量展開性質(zhì)進行分析。

(a) 本文方法計算結(jié)果

(b) Nastran仿真計算結(jié)果

本處以半充液工況為例，滿足α=β=0，且θ=φ+π/2。由于內(nèi)、外流場的主要區(qū)別在于貝塞爾函數(shù)類型不同，且控制方程式(23)中僅有矩陣[Vp]和[Tp]與貝塞爾函數(shù)有關(guān)，故重點在于分析式(18)和式(21)。這里假設(shè)權(quán)函數(shù)為聲壓的周向展開函數(shù)，由此式(18)和式(21)可以改寫為如下形式

(27)

(28)

當(dāng)a≠b時，式(27)、式(28)中正弦三角函數(shù)積分為零，由此可判斷矩陣[Vp]和[Tp]是對角矩陣，表達式如下

(29)

(30)

根據(jù)振動控制方程式(21)，僅與貝塞爾函數(shù)相關(guān)的矩陣可表示為[Tp] [Vp]-1，由此[Tp][Vp]-1也是對角矩陣

(31)

同理，半浸沒工況時[Tp][Vp]-1也有類似性質(zhì)，可推導(dǎo)出相同的近似表達式，這就從數(shù)學(xué)機理上闡述了為何半充液工況及半浸沒工況下同階次固有頻率相近。換句話說，當(dāng)圓柱殼處于半充液或半浸沒狀態(tài)時，流體負載可以近似簡化為相同的數(shù)學(xué)表達式，從而導(dǎo)致同階次固有頻率十分接近。

此外，由于聲壓負載關(guān)于特征高度是連續(xù)函數(shù)，固有頻率關(guān)于特征高度的函數(shù)曲線理論上也是光順的，因此可以推斷，在半充液或半浸沒附近時，相同液面高度下部分充液與部分浸沒對應(yīng)的同階次固有頻率也是非常接近的。從圖8可以大概看出，當(dāng)H/Rs在-0.2～0.2內(nèi)時，上述規(guī)律依舊存在。

為了更直觀、充分地了解該規(guī)律，分別取H/Rs=-0.2和0.2，計算部分充液工況和部分浸沒工況下前10階固有頻率，并開展對比研究，如圖10所示。

(a) H/Rs=-0.2

(b) H/Rs=0.2

由圖10可知，當(dāng)H/Rs=-0.2和0.2時，部分充液和部分浸沒工況下固有頻率十分接近。說明在半浸沒或半充液附近，內(nèi)、外流場在等液面高度時同階次固有頻率具有一致性。這一規(guī)律對相關(guān)的工程問題有一定的指導(dǎo)價值，例如采用試驗手段開展部分浸沒(H/Rs在-0.2～0.2內(nèi))圓柱殼結(jié)構(gòu)的模態(tài)試驗時，可以通過內(nèi)部充入等液面高度的相同流體來取得相近的試驗效果。且由于外流場試驗受到試驗場地的限制，其難度會遠高于內(nèi)流場試驗，采用內(nèi)流場試驗替代外流場試驗即經(jīng)濟性又便利。

5 結(jié) 論

本文基于兩套坐標(biāo)系的建模思路以及Galerkin方法，提出了有限長部分浸沒圓柱殼耦合振動的求解方法，并將該方法拓展到研究部分充液圓柱殼的耦合振動。此外，本文還對比分析了相同液面高度下，圓柱殼與內(nèi)、外流場耦合振動頻率的異同。具體結(jié)論如下：

(1) 本文方法的計算精度高且適用范圍廣，無量綱特征高度H/Rs可從-0.8取到0.8；此外，本文方法還具有計算量小的優(yōu)勢。

(2) 圓柱殼與流場部分耦合時，液面高度增大會導(dǎo)致濕表面增大，從而導(dǎo)致增加附連水質(zhì)量，使得同階次固有頻率下降。

(3) 當(dāng)圓柱殼處于半充液或半浸沒附近時，相同液面高度下同階次固有頻率十分接近。這一規(guī)律對于實際的工程問題有指導(dǎo)意義，如可以通過部分充液圓柱殼模態(tài)試驗來替代外流場試驗。