趙小西，衛(wèi)志軍，阮詩倫，岳前進，陳曉東（大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116024）

液體晃蕩荷載與液艙運動相位關(guān)系的數(shù)值研究

趙小西，衛(wèi)志軍，阮詩倫，岳前進，陳曉東
（大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116024）

該文從相位角度出發(fā)，采用CFD的數(shù)值方法研究簡諧激勵下液體晃蕩整體荷載與液艙運動的相位關(guān)系。首先，通過模型實驗驗證數(shù)值方法的可靠性，然后，以橫蕩激勵下的二維矩形艙晃蕩為例，研究不同載液率工況下，相位差和晃蕩整體荷載幅值隨激勵頻率的變化關(guān)系。結(jié)果表明：基于VOF法的數(shù)值模擬能較好地模擬液艙內(nèi)液體晃蕩；在自由面共振頻率附近（0.8 f0～1.2 f0），相位差隨激勵頻率增加而增大；頻率比相同時，相位差隨載液率的升高而增大；相位差90°時，相位差對應(yīng)頻率低于共振頻率，且與晃蕩荷載最大值對應(yīng)頻率基本相等。

液體晃蕩整體荷載；相位關(guān)系；數(shù)值模擬；局部沖擊荷載；模型實驗

0 引 言

帶液艙運載裝備在環(huán)境荷載下的運動激勵艙內(nèi)液體產(chǎn)生晃蕩，晃蕩引起的沖擊荷載作用于艙壁，進而影響工業(yè)結(jié)構(gòu)的整體運動姿態(tài)。特別在艙內(nèi)液體部分裝載時，劇烈晃蕩引起的抨擊壓力可能造成工業(yè)結(jié)構(gòu)的嚴重失穩(wěn)，但同時這種劇烈的晃蕩也可以加以合理利用，如船舶減搖水艙、調(diào)頻液體阻尼器（TLD）等的使用。液艙晃蕩對結(jié)構(gòu)整體運動產(chǎn)生不同的影響，與艙內(nèi)液體晃蕩整體荷載與結(jié)構(gòu)運動的相互作用關(guān)系密切相關(guān)。張力等[1]將調(diào)頻液體阻尼器（TLD）產(chǎn)生的晃蕩整體荷載等效成結(jié)構(gòu)附加阻尼，認為當(dāng)液體振蕩頻率與結(jié)構(gòu)運動的頻率接近或一致時，TLD產(chǎn)生與結(jié)構(gòu)運動相反的液體動壓力可以減小主體結(jié)構(gòu)的振動。馬維良、焦儂等[2]針對可控式被動減搖水艙，提出最佳相位PD控制策略，使艙內(nèi)液體的振蕩運動和船舶橫搖運動大致相差90°相位，滿足最佳減搖相位的要求，其減搖效果要優(yōu)于常規(guī)控制。Faltinsen和Timokha[3]在分析自由液面減搖艙時，提出晃蕩整體荷載與船體橫搖相位差接近90°時，具有較好的減搖效果，并據(jù)此通過已知的船體橫搖周期推出被動減搖艙的載液深。以上研究出于減振(搖)為目的，主要關(guān)注液體晃蕩整體荷載與結(jié)構(gòu)運動相位差90°情況，而對于相位差隨激勵頻率變化及載液率變化的規(guī)律卻缺乏細致深入的研究。針對這一問題，本文研究了簡諧運動激勵下，艙內(nèi)液體晃蕩整體荷載與液艙運動的相位差變化規(guī)律。

1 相位分析原理

考慮艙內(nèi)液體晃蕩影響的主體結(jié)構(gòu)運動，可簡化為在受到環(huán)境荷載作用的結(jié)構(gòu)上附加流體晃蕩整體荷載作用所產(chǎn)生的整體運動，結(jié)構(gòu)運動方程可以表達為：

式中：M、C和K分別為單自由度系統(tǒng)的質(zhì)量、阻尼和回復(fù)力（剛度）系數(shù)；Fw（t）為結(jié)構(gòu)外激激勵，F(xiàn)s（t）是運動過程中艙內(nèi)液體作用于結(jié)構(gòu)的晃蕩整體荷載。

結(jié)構(gòu)受到外界簡諧激勵x（t）=Acos（ω t）時，F(xiàn)altinsen等[3]將晃蕩整體荷載處理為：

其中：φ為晃蕩整體荷載滯后于結(jié)構(gòu)運動的相位差，ω為簡諧激勵的圓頻率。

（2）式可以變換為：

將（3）式代入結(jié)構(gòu)運動方程（1）可得：

從變形后的運動方程可以看出，可以將晃蕩整體荷載的影響，附加到結(jié)構(gòu)運動原來的阻尼和恢復(fù)力（剛度）系數(shù)上，令：

從（5）式可以看出，液體晃蕩附加阻尼Cs和附加恢復(fù)力（剛度）系數(shù)Ks的大小由液體晃蕩整體荷載幅值Fa和相位差φ決定，因此，確定Fa和φ的大小及變化規(guī)律是十分必要的。

2 液艙晃蕩數(shù)值模擬

2.1 VOF方法

與理論解析法相比，采用基于計算流體力學(xué)的數(shù)值模擬進行液艙晃蕩問題研究可以突破理論方法無法模擬波面翻轉(zhuǎn)等物理強非線性現(xiàn)象的限制，能夠模擬劇烈的流體晃蕩運動；與模型實驗相比，采用數(shù)值模擬方法更加方便，而且能夠?qū)问幀F(xiàn)象進行綜合分析。為了解決液體晃蕩中自由液面的變化描述，許多學(xué)者采用了不同的方法來捕捉和追蹤自由液面，通常采用的有VOF法[4-5]、MAC法[6]和LEVEL-SET法等[7]。其中VOF法的基本原理是通過研究網(wǎng)格單元中流體和網(wǎng)格體積比來確定自由面和追蹤流體的變化。VOF法使得流體的質(zhì)量得到守恒，并可跟蹤自由液面上發(fā)生翻轉(zhuǎn)、吞并、飛濺等復(fù)雜的自由面現(xiàn)象，并可以與各種離散方法結(jié)合，因而成為解決這一問題的重要工具[5，8]。

本文基于VOF法，采用流體計算商用軟件ANSYS FLUENT進行計算，數(shù)值計算設(shè)置如表1所示。

表1 數(shù)值計算設(shè)置Tab.1 Numerical calculation settings

2.2 液艙計算模型

本文采用的二維矩形液艙模型尺寸為L×H×B=970 mm×927 mm×158 mm，網(wǎng)格劃分間距為10 mm，共有144 336個網(wǎng)格，如圖1所示。

液艙做強迫橫蕩運動，橫向激勵為x（t）=Acos（ωt），其中：A=50 mm，ω為激勵圓頻率。ω表達式[3]為：

圖1 二維矩形液艙模型的網(wǎng)格剖分圖（單位：mm）Fig.1 Schematic diagram of 2D rectangular tank model(unit:mm)

3 液體晃蕩整體荷載與液艙運動

3.1 數(shù)值模擬方法驗證

液艙晃蕩的準確模擬是研究液體晃蕩荷載的前提。據(jù)此，通過開展液艙晃蕩模型實驗，驗證數(shù)值模擬結(jié)果的準確性。基于ANSYS FLUENT平臺的液體晃蕩整體荷載計算是局部沖擊壓強的積分，因此，通過驗證局部沖擊壓強同樣能反映液艙晃蕩整體荷載數(shù)值模擬的準確性。

驗證選取了20%H和50%H兩種載況(即液深分別為185 mm和464 mm)下局部沖擊壓強的時程進行對比。沖擊壓強監(jiān)測點及實驗液艙模型如圖2、圖3所示，監(jiān)測點的選取參考了衛(wèi)志軍等[10]面向晃蕩荷載研究的模型實驗，能夠反映液體劇烈晃蕩時的荷載響應(yīng)。

圖2 各監(jiān)測點位置（單位：mm）Fig.2 Locations of monitoring(unit:mm)

圖3 實驗液艙模型Fig.3 Tank model of experiment

圖4為載液50%H，激勵頻率f/f0為1.0時，p1（自由液面）和p2處沖擊壓強數(shù)值模擬與實驗值的對比。從圖中可以看出，不同監(jiān)測點處數(shù)值模擬結(jié)果與實驗值均吻合良好：盡管由于共振導(dǎo)致液體晃蕩的非線性使得晃蕩局部沖擊壓強在幅值上與實驗結(jié)果存在一定誤差，但壓強時程上的周期性變化趨勢卻基本同步。

圖4 50%H載液，f/f0=1.0時，沖擊壓強數(shù)值與實驗結(jié)果時程對比：（a）p1；（b）p2Fig.4 50%H filling level,f/f0=1.0,comparison between numerical simulation and experiment on pressure:(a)p1;(b)p2

圖5為載液20%H時，激勵頻率f/f0分別為0.8和1.0時，同一監(jiān)測點p3（自由液面）處沖擊壓強數(shù)值模擬與實驗值的對比。從圖中可以看出，非共振與共振情況下數(shù)值模擬與實驗值均吻合良好，尤其在非共振時，在進入穩(wěn)態(tài)后，無論沖擊壓強幅值大小還是壓強周期性變化趨勢均基本一致。

圖5 20%H載液，p3處沖擊壓力數(shù)值與實驗結(jié)果時程對比：（a）f/f0=0.8；（b）f/f0=1.0Fig.5 20%H filling level,comparison between numerical simulation and experiment on pressure with p3:(a)f/f0=0.8;(b)f/f0=1.0

上述對比結(jié)果說明，采用VOF法模擬液體晃蕩荷載是比較可靠的：盡管在接近共振時，液體晃蕩局部沖擊壓強幅值存在一定誤差，但沖擊壓強周期性變化趨勢基本同步，并且在非共振時，數(shù)值與實驗結(jié)果吻合良好。對于本文而言，數(shù)值模擬晃蕩局部沖擊壓強帶來的幅值誤差，是可以忽略的，因為劇烈晃蕩時的局部沖擊壓強峰值是高頻荷載，對于晃蕩整體荷載而言，其所包含的高頻荷載對結(jié)構(gòu)整體運動的影響而言是較小的[11]。因此，在處理基于ANSYS FLUENT監(jiān)測的晃蕩整體荷載時，可以采用恰當(dāng)?shù)臄M合方法獲得有效的晃蕩整體荷載。

3.2 計算工況選取

基于ANSYS FLUENT模擬的液體晃蕩整體荷載，監(jiān)測的是橫向左右艙壁上的液體晃蕩整體荷載。根據(jù)液體晃蕩特點，在不同載液率時，分別呈現(xiàn)不同特性，因此選取了三種具有代表性的載液率工況20%H、30%H和50%H。

當(dāng)外激激勵頻率接近第一階固有頻率時，艙內(nèi)液體發(fā)生共振，此時晃蕩非常激烈，這種極端工況是工程界研究液體晃蕩對整體結(jié)構(gòu)運動影響的重點，因此本文選取的激勵頻率在0.6 f0～1.3 f0，接近共振頻率會進行適當(dāng)加密。

3.3 晃蕩整體荷載計算結(jié)果處理

圖6為載液20%H，激勵頻率f/f0為1.0時，即共振時，液艙晃蕩整體荷載的變化曲線，其中虛線表示通過數(shù)值模擬的結(jié)果，可以看出，晃蕩整體荷載雖然呈周期性變化趨勢，但并不完全規(guī)則，為了便于數(shù)值表達，利用MATLAB對原始數(shù)據(jù)進行擬合處理。基于MATLAB的Fourier擬合表達式如下：

其中：a0，an，bn為常數(shù)，k表示擬合階數(shù)，k=1，2，3，···，8。分別對液體晃蕩整體荷載進行各階擬合，對比可知：

其中：a（0k）和ω(k)分別代表進行k階擬合時的a0和ω值，f為液艙運動的激勵頻率。

由此可見，液體晃蕩一階整體荷載的響應(yīng)頻率與液艙運動頻率一致，超過二階（包括二階）晃蕩整體荷載的響應(yīng)頻率大于液艙運動的頻率?？紤]高頻（相對于結(jié)構(gòu)運動頻率）荷載對結(jié)構(gòu)整體運動的影響較小[11]，可以忽略不計，因此，對晃蕩整體荷載采用一階擬合，結(jié)合（7）式和（8）式，擬合后的晃蕩整體荷載表達式為：

其中：ω為液體晃蕩的圓頻率。

圖6中實線即為采用一階擬合后液體晃蕩整體荷載的變化曲線，擬合后的晃蕩整體荷載變化曲線能夠很好地反映晃蕩整體荷載的變化，且其變化頻率與液艙運動的頻率相等，為晃蕩整體荷載與液艙運動的相位差分析提供了可行性。

圖6 20%H載液，共振時，流體晃蕩整體荷載原始和擬合后變化曲線Fig.6 20%H filling level,raw and fitting global force with an resonance frequency

3.4 相位差計算

圖7為載液30%H，激勵頻率f/f0分別為0.9、0.95、1.0和1.2時，液體晃蕩整體荷載與液艙運動位移的時程對比曲線。從圖中可以看出，激勵頻率不同，晃蕩整體荷載與液艙運動在時域內(nèi)的相對變化也不相同，這種不同可以通過液艙晃蕩整體荷載滯后于液艙運動的相位差直觀反映。相位差的計算如下：

其中：ΔT為晃蕩整體荷載（擬合）滯后于液艙運動的時間之差，f為激勵頻率。

3.5 計算結(jié)果分析

圖8為三種不同載液率下，晃蕩整體荷載與液艙運動位移相位差隨頻率變化的曲線。從圖中可以看出，相位差均隨頻率呈遞增變化，變化范圍為0°～180°，且變化主要發(fā)生在共振區(qū)間段（0.8 f0～1.2 f0）；在共振區(qū)間，同一激勵頻率比對應(yīng)的相位差隨載液率的升高而增大；相位差90°對應(yīng)的頻率低于共振頻率：載液20%H時，發(fā)生在頻率0.95f0附近；載液30%H，發(fā)生在頻率0.9f0附近；載液 50%H時，發(fā)生在0.85 f0附近。

圖9為三種載液率下，晃蕩整體荷載大小隨頻率變化的曲線。從圖中可以看出，整體荷載呈單峰變化，峰值對應(yīng)的頻率隨載液的不同而不同：20%H載液時，發(fā)生在共振頻率附近，載液30%H和50%H時，則分別發(fā)生在頻率0.85 f0和0.9 f0附近。

圖9 載液不同，流體晃蕩整體荷載隨頻率變化曲線Fig.9 The relationship between sloshing-induced global force and oscillation frequency

綜合圖8、圖9可知，液體晃蕩整體荷載與液艙運動相位及整體荷載幅值變化主要發(fā)生在共振區(qū)間；相位差90°對應(yīng)頻率低于共振頻率；且相位差90°對應(yīng)頻率與晃蕩整體荷載峰值對應(yīng)頻率：20%H載液率時，前者對應(yīng)頻率低于后者頻率約0.05 f0，30%H和50%H載液率時，兩者頻率相等。

4 結(jié) 論

本文采用基于CFD的數(shù)值模擬方法，對二維矩形艙橫蕩激勵下的液體晃蕩整體荷載與液艙運動的相位差及荷載幅值隨外激激勵頻率的變化關(guān)系進行研究。通過對比數(shù)值模擬與模型實驗所監(jiān)測的局部沖擊荷載值時程，驗證了基于VOF法的數(shù)值模擬方法的可靠性。然后利用數(shù)值模擬方法，對比不同載液率下相位差及晃蕩整體荷載的數(shù)值模擬結(jié)果，獲得以下結(jié)論：

（1）各載液率下，相位差均隨激勵頻率增加而增大，在共振區(qū)間（0.8 f0～1.2 f0）對頻率變化敏感，變化范圍為0°～180°；

（2）載液不同，同一頻率比對應(yīng)的相位差也不相同，隨著載液的升高呈增大趨勢；

（3）相位差為90°時，相位差對應(yīng)頻率低于共振頻率，并隨著載液率的升高，偏離程度呈增大趨勢：20%H，30%H和50%H載液時，分別低于共振頻率約0.05 f0，0.1 f0和0.15 f0。

（4）相位差90°時，相位差對應(yīng)頻率與晃蕩整體荷載最大值對應(yīng)頻率基本相等（20%H載液時，前者低于后者約0.05 f0）。

根據(jù)本文研究表明，液艙晃蕩整體荷載與結(jié)構(gòu)運動的相位關(guān)系及整體荷載大小與結(jié)構(gòu)運動頻率及載液率存在對應(yīng)關(guān)系。通過分析相位差和晃蕩整體荷載大小隨結(jié)構(gòu)運動頻率及載液率的變化規(guī)律，然后結(jié)合結(jié)構(gòu)運動方程，能夠用于進一步分析液艙晃蕩對結(jié)構(gòu)整體運動姿態(tài)影響。而且，通過對液艙晃蕩產(chǎn)生的附加阻尼Cs和附加剛度Ks等參數(shù)進一步研究，還可以直接為減搖艙及TLD設(shè)計提供參考。

[1]張 力.導(dǎo)管架海洋平臺冰激振動控制的實驗研究[D].大連:大連理工大學(xué),2008. Zhang Li.Experimental study on ice-induced vibration mitigation of jacket offshore platforms[D].Dalin:Dalian University of Technology,2008.

[2]馬維良,焦 儂.可控式被動減搖水艙最佳相位PD控制策略研究[J].船舶工程,2012,1(34):36-39. Ma Weiliang,Jiao Nong.Research on optimal phase PD control strategy of passively controlled anti-rolling tank[J].Ship Engineering,2012,1(34):36-39.

[3]Faltinsen O M,Timokha A N.Sloshing[M].Cambridge:Cambridge University Press,2009.

[4]Hirt C W,Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries[J].Journal of Computational Physics,1981,39:201-225.

[5]劉楨兵.基于VOF法的液艙晃蕩數(shù)值模擬及載荷計算[J].重慶理工大學(xué)學(xué)報（自然科學(xué)）,2011,25(3):24-29,61. Liu Zhenbing.Numerical simulation of liquid tank sloshing and pressure calculation based on Volume of Fluid method[J]. Journal of Chong Qing University of Technology(Natural Science),2011,25(3):24-29,61.

[6]Mikells N E,Journee J M J.Experimental and numerical simulations of sloshing behavior in liquid cargo tanks and its effect on ship motions[C]//National Conference on Numerical Methods for Transient and Coupled Problems,July 1984. Venice,Italy,1984:274-278.

[7]Zhu Renqing,Wu Yousheng,Atilla N.Numerical simulation of liquid sloshing-a review[J].Shipbuilding of China,2004,45 (2):14-27.

[8]沈 猛,王 剛,唐文勇.基于改進VOF法的棱形液艙液體晃蕩分析[J].中國造船,2009,50(1):1-8.Shen Meng,Wang Gang,Tang Wenyong.Numerical analysis of liquid tank sloshing based on Volume of Fluid method[J]. Shipbuilding of China,2009,50(1):1-8.

[9]Thiagarajan K P,Rakshit D,Repalle N.The air-water sloshing problem:Fundamental analysis and parametric studies on excitation and fill levels[J].Ocean Engineering,2011,38:498-508.

[10]衛(wèi)志軍,岳前進,張文首,等.大尺度儲艙內(nèi)流體晃蕩砰擊壓力的測量方法研究[J].中國科學(xué),2014,44(07):746-758. Wei Zhijun,Yue Qianjin,Zhang Wenshou,et al.Experimental investigation of violent slamming pressure in large-scaled tank[J].SCIENTIA SINICA Phys,Mech&Astron,2014,44(07):746-758.

[11]Huang Z J,Danaczko M A,Esenkov O E,Martin C B,O’Donnell B J,Yung T W.Coupled tank sloshing and LNG carrier motions[C]//The International Society of Offshore and Polar Engineers Conference(ISOPE),June 21-26,2009.Osaka, Japan,2009:92-104.

Numerical study on phase relation between sloshing-induced global force and tank motion

ZHAO Xiao-xi,WEI Zhi-jun,RUAN Shi-lun,YUE Qian-jin,CHEN Xiao-dong
(State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology,Dalian 116024,China)

From phase analysis point of view,this paper presents the phase relation between the sloshing global load and tank motion under harmonic forced excitation based on the CFD method.And the sloshing model tests were carried out to verify the reliability of the numerical method.Then,the phase differences between the sloshing global load and the motion of a 2D-rectangular tank and the amplitudes of sloshing global load with forced-frequencies are obtained for different filling levels.The results show that the tank sloshing can be simulated well by VOF method;the phase differences raise when increasing forced-frequencies around resonant frequencies of free surface(0.8 f0～1.2 f0);the phase differences also increase for higher filling levels at same frequency ratio;the frequency for 90-degree phase difference is less than the resonant frequency and it is roughly equal to that for the maximum sloshing global load.

sloshing global load;phase difference;numerical simulation;pressure;model test

U661.3

A

10.3969/j.issn.1007-7294.2015.08.005

1007-7294（2015）08-0919-08

2015-03-14

國家重大科技子課題專項（2011ZX05026-006-06）；創(chuàng)新研究群體基金資助（50921001）

趙小西（1987-），男，碩士研究生；

阮詩倫（1976-），男，副教授，碩士生導(dǎo)師，E-mail:ruansl@dlut.edu.cn。