韓越洋， 朱翔,3， 郭文杰， 李天勻,3

(1.華中科技大學(xué) 船舶與海洋工程學(xué)院，湖北 武漢 430074；2.船舶與海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430074；3.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心， 上海 200240；4.華東交通大學(xué) 鐵路環(huán)境振動(dòng)與噪聲教育部工程研究中心，江西 南昌 330013)

圓柱殼作為一種基本結(jié)構(gòu)單元，被廣泛應(yīng)用于液體儲(chǔ)存及運(yùn)輸方面。對于許多液體散裝運(yùn)輸設(shè)備，如遠(yuǎn)洋運(yùn)輸船、公路鐵路油罐車和太空運(yùn)輸飛船等，其方向穩(wěn)定性、安全性以及結(jié)構(gòu)的可靠性很容易受到內(nèi)部部分填充的液體的影響。在接近共振的情況下，由液體晃蕩產(chǎn)生的局部力和局部力矩可能會(huì)導(dǎo)致結(jié)構(gòu)的不穩(wěn)定，甚至圓柱殼結(jié)構(gòu)和內(nèi)部流體相互作用影響到結(jié)構(gòu)設(shè)備性能，可能成為一些事故的原因[1-4]。學(xué)者通過數(shù)值方法等對殼體中的液體運(yùn)動(dòng)進(jìn)行了研究[5-10]。其中Wang[5]利用幾何比例邊界有限元法，研究了水平圓柱體內(nèi)部的液體晃動(dòng)。Liu等[6]采用運(yùn)動(dòng)網(wǎng)格結(jié)合流體體積方法來捕捉晃動(dòng)過程中液態(tài)氫罐自由表面的波動(dòng)，研究了在正弦激勵(lì)下晃蕩的液體的水動(dòng)力性能。文獻(xiàn)[11-23]利用解析法、半解析法研究殼體中液體的晃蕩。Evens等[11]利用伽遼金法構(gòu)造矩陣求解特征值，對剛性邊界下水平半充液圓柱殼的內(nèi)部液體的晃蕩頻率進(jìn)行了計(jì)算。McIver[12]采用二維雙極坐標(biāo)系并且將特征值問題轉(zhuǎn)化積分方程求解問題，計(jì)算了水平二維充液圓柱殼的液體晃蕩頻率。Hasheminejad等[13]利用連續(xù)的共形坐標(biāo)變換提出了精確的二維流體動(dòng)力學(xué)分析，分析了內(nèi)部有擋板的水平部分充液圓柱殼的液體橫向晃蕩特性。Papaspyrou[14-15]根據(jù)Evans的半解析法，先后分析了橫向、縱向外部激勵(lì)下，水平半充液圓柱殼的晃蕩響應(yīng)。文獻(xiàn)[24-26]采用實(shí)驗(yàn)手段對水平充液殼體中液體的晃蕩運(yùn)動(dòng)進(jìn)行了研究。

針對二維充液殼體中液體晃蕩的解析解或者半解析解大多只針對半充液圓柱殼或者矩形截面。由于建模和求解上的難度，針對充液深度變化的圓柱殼截面自由液面晃蕩的解析方法，涉及文獻(xiàn)較少。本文針對水平任意充液深度的二維圓柱殼中液體晃蕩運(yùn)動(dòng)求解，提出一種半解析法。首先在圓柱殼和自由液面上分別建立2個(gè)極坐標(biāo)系，根據(jù)2個(gè)坐標(biāo)系的變換關(guān)系將圓柱殼和自由液面的2個(gè)坐標(biāo)系聯(lián)系起來。根據(jù)拉普拉斯方程得出自由液面的速度勢函數(shù)，進(jìn)而根據(jù)邊界條件及坐標(biāo)關(guān)系，得出自由液面的運(yùn)動(dòng)控制方程，并通過伽遼金法解出自由液面晃蕩頻率。通過改變充液深度，對二維圓柱殼自由液面晃蕩頻率的特性進(jìn)行了分析。

1 二維任意深度充液剛性圓柱殼晃蕩理論模型

如圖1所示，二維圓柱殼上建立以殼體圓心O為坐標(biāo)原點(diǎn)的極坐標(biāo)系(r,φ)，且φ取值范圍為[-α,α]。自由液面上則建立以液面中心Q為坐標(biāo)原點(diǎn)的極坐標(biāo)系(R,θ)，θ取值范圍為[-π/2,π/2]。H為自由液面坐標(biāo)系到圓柱殼結(jié)構(gòu)坐標(biāo)系的距離，H值為正則液面高度小于半徑；H值為負(fù)則液面高度大于半徑。2個(gè)坐標(biāo)系在空間上任意一點(diǎn)的夾角定義為β=φ-θ。當(dāng)H值為正時(shí)，∣β∣=θ-φ；當(dāng)H值為負(fù)時(shí)，∣β∣=φ-θ。rs為圓柱殼半徑。H值與充液深度h的關(guān)系為h=rs-H。

圖1 二維剛性充液圓柱殼坐標(biāo)系示意

假設(shè)液體為無粘性且不可壓縮的理想流體，則對于小幅晃蕩的理想流體，其速度勢函數(shù)應(yīng)滿足拉普拉斯方程：

(1)

式中φ(r,θ)為液體速度勢函數(shù)。

由于液面為小幅晃蕩,還需要滿足一維線性波的表面邊界條件:

(2)

式中：θ=±π/2 ；K=ω2/g；ω為液體晃蕩的圓頻率。

假設(shè)殼體為剛性殼體，則速度勢函數(shù)φ(r，θ)應(yīng)滿足流體與剛性壁面的邊界條件，即：

(3)

其中r=rs，rs為圓柱殼半徑。

通過上述約束條件，自由液面晃蕩的速度勢函數(shù)的解可利用分離變量法寫為對稱/反對稱2種形式的解。首先，速度勢函數(shù)滿足拉普拉斯方程(1)，則速度勢函數(shù)反對稱與對稱模態(tài)的解分別為：

(4)

(5)

對于自由液面反對稱模態(tài)形式，為后續(xù)方便計(jì)算可將式(4)改寫為：

(6)

將式(6)代入式(2)中，可得：

(7)

從式(7)中可得A2n-1與A2n之間的關(guān)系為A2n=-K/(2n)A2n-1，代入式(6)中，可得：

(8)

將式(3)從結(jié)構(gòu)坐標(biāo)系轉(zhuǎn)換至液面坐標(biāo)系上，即：

(9)

由圖1可知結(jié)構(gòu)坐標(biāo)系和液面坐標(biāo)系的關(guān)系：

R2=H2+r2-2Hrcosφ

(10)

Rsinθ=rsinφ

(11)

則可得液面坐標(biāo)系中坐標(biāo)R和θ對結(jié)構(gòu)坐標(biāo)系r求導(dǎo)關(guān)系式：

(12)

(13)

將式(8)、(12)以及(13)代入式(9)可將邊界條件為：

2)θ+φ]-KR2n-1sin[(2n-1)θ+φ]]}=0

(14)

(D-KT){A2n-1}=0

(15)

晃蕩頻率的求解可利用齊次線性方程組系數(shù)矩陣行列式為零進(jìn)行求解，即：

|D-KT|=0

(16)

利用搜根法即可解出晃蕩頻率。同理，對于自由液面對稱模態(tài)，計(jì)算方法與反對稱模態(tài)計(jì)算方法類似，首先將速度勢函數(shù)改寫：

(17)

速度勢函數(shù)需要滿足約束條件式(2)，可得；

(18)

則速度勢函數(shù)可寫為：

(19)

將式(12)、(13)以及(19)代入式(9)中可得：

KR2ncos[2nθ+φ]]}=0

(20)

(L-KY){B2n}=0

(21)

其中矩陣L=[lmn]和矩陣Y=[ymn]的元素分別為：

此時(shí)晃蕩頻率也可以通過令式(21)中系數(shù)矩陣行列式為零求解得到。

2 收斂性分析及正確性驗(yàn)證

通過本文的理論推導(dǎo)即可以得出自由液面晃蕩頻率。由于勢函數(shù)是通過無窮級數(shù)表達(dá)，因此有必要針對級數(shù)的截?cái)囗?xiàng)數(shù)進(jìn)行收斂性分析。給定圓柱殼半徑rs=1 m, 液面高度h為0.5 m，分別計(jì)算m，n的截?cái)鄶?shù)N=5,10,20,30,40時(shí)，殼體自由液面無量綱晃蕩頻率Ω=Krs，如表1所示。很明顯在N=20時(shí)晃蕩頻率已經(jīng)穩(wěn)定。從表中得出N=40時(shí)可以得到更精確的值，但計(jì)算時(shí)間也會(huì)大大延長。所以本文采用截?cái)鄶?shù)N=20，既能保證精度，又可以節(jié)省時(shí)間。

表1 二維剛性充液圓柱殼自由液面的無量綱晃蕩頻率Ω收斂性分析

為了驗(yàn)證本文理論模型和算法的正確性，計(jì)算二維剛性充液圓柱殼不同充液深度下自由液面的晃蕩頻率，并同文獻(xiàn)結(jié)果進(jìn)行對比分析。計(jì)算參數(shù)如下：圓柱殼半徑為rs=1 m，H分別為0.5,0,-0.5 m，對應(yīng)充液深度h分別為0.5,1,1.5 m。由方程(15)和(21)可以得出不同充液深度下二維剛性充液圓柱殼自由液面的晃蕩頻率，如表2所示。其中q為模態(tài)階數(shù)。

表2 不同充液深度下二維充液圓柱殼自由液面的無量綱Ω晃蕩頻率對比

表3 二維剛性充液圓柱殼自由液面的晃蕩頻率f有限元和本文方法結(jié)果對比

進(jìn)一步將本文計(jì)算的晃蕩振型和有限元結(jié)果進(jìn)行對比分析。圖2為H=0.5 m時(shí)流體有限元單元模型。圖3為本文方法和有限元法計(jì)算得到的自由液面前4階模態(tài)振型對比，其中有限元計(jì)算出的振型根據(jù)模型中自由液面節(jié)點(diǎn)數(shù)據(jù)繪制得出。從圖3可見本文計(jì)算得到的自由液面模態(tài)振型和有限元計(jì)算的模態(tài)振型也吻合很好。由此可見，和已有文獻(xiàn)晃蕩頻率以及和有限元法晃蕩頻率和振型對比，都表明本文理論方法的準(zhǔn)確性及可靠性。

圖2 二維剛性充液圓柱殼中流體有限元網(wǎng)格模型

圖3 二維剛性充液圓柱殼自由液面的晃蕩振型

3 不同充液深度下自由液面晃蕩頻率分析

對殼體半徑為rs=1 m，不同充液深度下(H=0.7 m，0.5 m，0 m，-0.5 m，-0.7 m)二維剛性充液圓柱殼自由液面晃蕩頻率進(jìn)行計(jì)算并對比分析。

圖4(a)給出了不同充液深度下晃蕩頻率隨模態(tài)階數(shù)q的變化。圖4(b)則給出了晃蕩頻率隨充液深度的變化曲線。由圖4(a)可知，充液深度一定時(shí)，充液圓柱殼自由液面晃蕩頻率隨著模態(tài)階數(shù)q的增大而增加。根據(jù)文獻(xiàn)[27]提出的晃動(dòng)流體等效力學(xué)模型分析，隨著模態(tài)階數(shù)增加，內(nèi)部液體晃蕩的固定質(zhì)量(即不參與晃蕩的流體質(zhì)量)占比逐漸增大，而晃蕩質(zhì)量占比則逐漸減小，晃蕩固有頻率逐漸增大。對于不同的充液深度，以上結(jié)論也都成立。

圖4 二維剛性充液圓柱殼自由液面晃蕩頻率

結(jié)合圖4可知，在第1階模態(tài)的時(shí)候，晃蕩頻率隨著深度的增大而增大。這是由于在1階模態(tài)時(shí)，隨著液深的增加，內(nèi)部流體的固定質(zhì)量占比逐漸增大，晃蕩程度減小，晃蕩質(zhì)量占比逐漸降低，因此晃蕩頻率呈增大的規(guī)律。

但是在較高階模態(tài)時(shí)，如第5階模態(tài)時(shí)，晃蕩頻率則是隨著充液深度的增大先減小后增大，且在半充液的時(shí)候出現(xiàn)轉(zhuǎn)折。這可以解釋為在高階晃蕩模態(tài)中，液面晃蕩的節(jié)點(diǎn)增多，而在圓柱殼中液深關(guān)于殼體圓心上下的變化會(huì)導(dǎo)致液體自由表面區(qū)域的變小，相對剛度增大，使得激發(fā)高階晃蕩難度變大，則晃蕩頻率增大。

并且，在H=0.7 m時(shí)的晃蕩頻率分別在第2階模態(tài)時(shí)超過H=0.5 m時(shí)的晃蕩頻率，第3階模態(tài)的時(shí)候超過H=0 m的晃蕩頻率，第5階模態(tài)的時(shí)候超過H=-0.5 m的晃蕩頻率。綜上得出充液深度越小，晃蕩頻率隨著模態(tài)階數(shù)遞增的速率也越快，圖4(a)中曲線的斜率中也可直觀反映出這一結(jié)論。

4 結(jié)論

1)特定深度下，自由液面晃蕩頻率隨著模態(tài)階數(shù)的增大而增大；在模態(tài)階數(shù)q=1時(shí)，自由液面晃蕩隨充液深度的增大而增大；在高階模態(tài)中，自由液面晃蕩頻率隨著充液深度的增大先減小后增大，且在半充液狀態(tài)時(shí)發(fā)生轉(zhuǎn)折；不同充液深度下，自由液面晃蕩頻率增大速率隨著充液深度的減小而增大。

2)本文的半解析法與有限元法相比建模簡便，參數(shù)分析時(shí)無需重復(fù)建模，計(jì)算結(jié)果準(zhǔn)確、可靠。與其他解析方法相比，具有良好的準(zhǔn)確性和適應(yīng)性，能適應(yīng)不同充液深度下圓柱殼的晃蕩頻率求解，并且能進(jìn)一步拓展到三維圓柱殼中的自由液面晃蕩問題以及彈性圓柱殼的流固耦合振動(dòng)分析。