彈性液艙液體晃蕩數值模擬

周上然，朱仁慶

(江蘇科技大學船舶與海洋工程學院，江蘇 鎮江 212003)

彈性液艙液體晃蕩數值模擬

周上然，朱仁慶

(江蘇科技大學船舶與海洋工程學院，江蘇 鎮江 212003)

以流體的連續方程、N-S方程與結構動力學方程為控制方程，利用動力學軟件ANSYS14.0，模擬二維矩形液艙與三維棱形液艙液體晃蕩現象，用VC++語言編寫矩形艙壁變形運動的響應程序，分析二維矩形液艙艙壁變形對晃蕩壓強的影響；利用System Coupling 模塊實現三維棱形彈性液艙晃蕩的雙向耦合；通過提取監測點處壓強，對比壓強的時間歷程，分析了不同板厚、不同材質對晃蕩載荷的影響及不同監測點處壓強變化特點。

液體晃蕩；雙向耦合；晃蕩載荷

大型載液船舶在航行時容易引起液艙晃蕩。液艙晃蕩引起的晃蕩載荷、結構響應以及兩者之間的耦合效應，一直是大家研究和關注的焦點。研究液艙晃蕩耦合效應有3種途徑：理論分析、數值和實驗模擬。目前，數值模擬方法運用最為廣泛。李青[1](2004年)基于ALE方法，對粘性流體與彈塑性薄壁結構流固耦合問題進行求解。方智勇[2]等(2006年)將Level-set方法與通度系數方法結合，模擬了二維矩形液箱受迫晃蕩運動。祁江濤[3](2007年)利用VOF法描述晃蕩流場的自由液面運動，設定動網格參數，對不同幾何形狀的液艙晃蕩進行仿真計算。端木玉[4](2007年)采用改進的Youngs 法，計算了不同液艙結構形式液艙晃蕩效應。沈猛[5](2008年)采用部分單元參數法對傳統VOF數值方法進行改進，對棱形液艙模型在不同裝載水平于不同橫搖激勵頻率條件下進行仿真計算。侯玲[6](2009年)基于VOF模型模擬了三維液艙晃蕩，用Fluent軟件的自定義函數編寫程序，模擬了流體與結構之間的相互作用。本文以流體的連續方程、N-S方程與結構動力學方程為控制方程，利用動力學軟件ANSYS14.0，模擬二維矩形液艙與三維棱形液艙晃蕩現象，用VC++語言編寫艙壁變形運動的響應程序，分析二維矩形液艙艙壁變形對晃蕩壓強的影響；通過對比監測點壓強時間歷程，分析參數變化對晃蕩載荷的影響及不同監測點處壓強變化特點。

1 數值方程

1.1流體運動方程

若假設流體為不可壓縮粘性流，則在直角坐標系下流體的連續方程和動量守恒方程分別為：

連續方程：div(v)=0

(1)

(2)

式中：v為流場中任意一點的速度矢量；ρ為流體密度；μ為流體的粘性系數；Fb為作用在單位體積流體上質量力；P為流場壓強。

1.2結構運動的微分方程

(3)

1.3流體與結構交界面條件

流固耦合交界面應滿足運動學連續條件、動力學連續條件、能量守恒原理。運動學條件指在流固界面上，位移與速度一致，即：

ds=dfus=uf

(4)

動力學連續條件指在流固界面上，任意一點力的平衡，即：

σs·ns=σf·nf

(5)

能量守恒指耦合作用過程中，耦合界面上流體載荷與固體力在界面位移上所做虛功相等，即：

δus·fs=δuf·ff

(6)

式(1)～式(6)中：ds、df分別為結構、流體在耦合界面上的位移；us、uf分別為結構、流體在耦合界面上的速度；σs、σf分別為結構、流體在耦合界面上的應力；ns、nf分別為結構、流體在耦合界面上外法線方向矢量；fs、ff分別為耦合界面上固體、流體表面力；δus、δuf分別為耦合界面上固體、流體的虛位移。

2 二維液艙晃蕩耦合數值模擬

考慮艙壁變形對晃蕩壓力的影響，以二維艙室作為研究對象，為了便于對艙壁變形運動進行編程，將艙室的側壁簡化為梁結構，結構運動的微分方程(3)可以寫成以下形式：

(7)

利用ANSYS14.0建立二維矩形液艙模型，模型尺寸：長度0.7 m，高度0.5 m，裝載深度0.35 m，如圖1。流體介質選用水。模擬高裝載狀態下的橫蕩運動，設定液艙右艙壁為柔性艙壁，可以變形，左艙壁為剛性壁；將右艙壁變形運動簡化為單質點運動，運用Newmark-β法，忽略金屬的阻尼，采用C語言編寫質點運動的響應程序，選取適當的形函數，得出艙壁各點的運動響應，修正右艙壁網格位置，以UDF(自定義函數)為平臺，結合動網格技術，采用擴散光順方法，數值模擬液艙橫蕩。運動工況見表1；提取左右艙壁中間點處的壓強進行對比，如圖2所示。

圖1 二維矩形液艙模型

表1 模型的運動工況

圖2 左艙壁與右艙壁監測點處壓強時間歷程

從圖2可以看出，右艙壁監測點處的壓強峰值低于左艙壁監測點處，艙壁的運動變形吸收了沖擊能量，降低了晃蕩壓強峰值；右艙壁監測點處的壓強峰值之間的過渡呈現一定的“振蕩”，因為作用在右艙壁上的液動壓力引起了艙壁的動態響應，沖擊壓力的峰值衰減之后，艙壁運動產生的慣性載荷為主要作用載荷。在相同外激頻率，艙壁的幾何尺寸相同的情況下，剛性結構不吸收沖擊能量，不發生運動變形；而彈性結構吸收了液體的部分沖擊能量，會發生運動和變形，轉化為結構的動能與彎曲應變能，體現了一定的動態效應。動態效應可以用動力放大系數Rd[7]來表示，計算方法如式(8)：

(8)

式中：ωn為結構的固有頻率；ω為液動壓力的變化頻率；ξ為結構的阻尼系數，一般情況下，金屬材料的取值較小，可以忽略。

ω越接近ωn，動態效應越顯著。計算ωn時，要考慮附加質量的影響，通常附加質量會降低結構的固有頻率。

3 三維棱形液艙晃蕩耦合模擬

3.1棱形液艙模型的建立

計算模型為按照實尺度1：55進行縮比的液艙模型,如圖3所示。圖中,l、b、h為液艙模型長度、寬度、高度，分別為0.8、0.7、0.5 m;ru、rl為斜板上傾角與下傾角，大小取135°；hu、hl為斜頂板與斜底板的高度，大小分別為0.13、0.09 m。模型的材質選用船用鋼與鋁2種材質。利用ANSYS14.0 System Coupling模塊實現液艙與晃蕩的雙向耦合。動網格光順方法采用擴散光順，液艙網格單元采用四邊形單元，流體介質選用水。流體網格采用六面體單元。運動工況見表2。液艙模型設立5個監測點，位置見表3。

圖3 液艙簡化模型

3.2不同模型材質對晃蕩壓強的影響

取表2中工況1、工況2進行數值模擬，計算時間取5個周期：6 s。自由液面波形圖見圖4。分別提取鋼質模型與鋁質模型監測點壓強時間歷程，如圖5與圖6所示；對比監測點壓強峰值，見表4。

表2 模型的運動工況

表3 監測點在模型中位置

圖4 模型自由液面波形圖

表4監測點P1至P4壓強結果比較

材質P1點的壓強幅值P2點的壓強幅值P3點的壓強幅值P4點的壓強幅值鋁5600Pa4100Pa3950Pa3150Pa船用鋼6250Pa4750Pa4800Pa3900Pa

圖5 鋼質模型監測點(P1至P4)的壓強時間歷程

從圖4可以看出，在高裝載的液艙晃蕩過程中，超過一定深度以后，在較深處的液體隨著艙室一起運動；隨著深度的增加，波能不斷衰減。

從圖5與圖6可以看出，底部區域監測點的壓強變化周期性特征明顯，基本上與外激頻率接近。自由液面附近的監測點壓強變化頻率較大，與艙壁的作用時間短，壓強雙峰特征明顯，第1個峰值為液體沖擊力，第2個峰值為沖擊之后。因為結構的運動慣性產生的作用力，與結構的剛度、液體的沖擊速度、結構的動態響應相關；壓強周期性變化以振蕩形式過渡；鋁材質艙壁的受力明顯低于鋼材質艙壁，底部區域監測點的壓強相差值比例低于自由液面附近的監測點的差值比例。因為隨著深度的增加，艙壁對液體運動的“限制”明顯增加，液體運動受限，運動沒有被完全激發，所以壓力以靜壓為主；在自由液面處，液體運動激烈，壓力主要成分為沖擊壓強。

從表4中可以看出，2種不同的材質，相同的板厚，鋁材質模型監測點的壓強幅值要低于鋼材質模型監測點的壓強幅值。相同尺寸的相況下，2種材質的屈曲強度相差在一定的范圍，彈性模量較小的材質能夠體現出一定彈性效應。剛性結構的固有頻率大于彈性結構的固有頻率；越接近結構的固有頻率，則結構的動態效應越顯著，所以在同等情況下，彈性結構動態效應較大，吸收能量較多。

3.3不同板厚對艙室結構應力的影響

取表2中工況3、工況4、工況5進行模擬。計算時間取3個周期：3.6 s。自由液面波形圖如圖7所示。提取1.2 s時刻3種工況下模型應力云圖與P5點一個周期內時間歷程，如圖8、圖9所示。對比監測點P1、P3壓強峰值與P5一個周期內的應力峰值，見表5。

表5 監測點P1、P3、P5的結果比較

圖6 鋁質模型監測點(P1至P4)的壓強時間歷程

圖7 模型自由液面波形圖

圖8 1.2 s時工況3、工況4、工況5模型應力圖

從圖7可以看出，在高裝載深度的情況下，液體晃蕩極易表現出大幅駐波現象，液體沖頂現象明顯。從圖8可以看出，對棱型液艙來說，處于同一水平位置上的內應力越是靠近橫縱艙壁交界處，其值越大。因為板交界處與角隅處受到的約束較多，所以相比模型其他部位應力較大；其次，液體晃蕩沖擊壓力作用于艙頂，艙頂受到的沖擊壓力作用要顯著，所以艙頂及各艙壁的交界處的應力明顯大于其他部位。隨著厚度的增加，模型的內應力隨之降低。因此，在設計制造棱型液艙時，要對艙室的舷側結構進行加強，重點關注角隅處及縱橫艙壁交界處。矩形板條梁的內應力近似計算：

(9)

式中：Mw為矩形板內彎矩；I為截面慣性矩；Z為截面上點距中性軸的距離。

整理式(9)，得到應力近似計算公式：

(10)

根據監測點壓強圖可以看出，隨著板厚的增加，監測點的壓強峰值略微增大。模型中舷側板的長寬之比大于2.5，可以按照板條梁模型計算板的柔度與應力。3塊板的模型的柔度均大于100，在柔性范圍之內。對于矩形板，相同的外部激勵，相同的長度與寬度材料屬性不變，板厚超過一定的范圍，板體現出剛性效應，板厚對沖擊壓強的影響不大，見表5。因為隨著板厚的增加，板的固有頻率不斷增大，固有周期不斷減小，動力放大系數減小；當液動壓力的增長時間超過結構的固有周期時，沖擊載荷對結構的作用表現為準靜態加載,此時結構響應與靜力分析結果相似。這種情況下，舷側板在失穩之前，內應力σ可以按照式(10)計算，且與厚度的平方t2成反比；如果板的厚度過小，部分區域承受的內應力達到失穩應力，變成塑性板，則喪失了彈性效應。

圖9 工況3、工況4、工況5模型P5點應力時間歷程

4 結語

(1)在相同外激頻率，艙壁的幾何尺寸相同的情況下，剛性結構的固有頻率大于彈性結構的固有頻率；越接近結構的固有頻率，則結構的動態效應越顯著，所以在同樣的液動壓力的變化頻率的情況下，彈性結構的動態效應較大，吸收了液體的部分沖擊能量，會發生運動和變形，轉化為結構的動能與彎曲應變能，降低了沖擊作用力；沖擊壓力的峰值衰減之后，艙壁運動產生的慣性載荷為主要作用載荷。

(2)對比各監測點壓強歷程，距離自由液面較深的點壓強變化相對平緩，靜壓起主要作用，壓強增長時間較長，變化周期與液艙的運動周期接近；艙室頂點主要承受沖擊壓強作用，具有脈沖性特點，壓強的增長時間短，能夠引起結構動態響應，沖擊壓強的幅值具有隨機性的特點；對于棱形艙室，處于同一水平位置上的晃蕩壓力越是靠近橫縱艙壁交界處其值越大，因此要對交界處加強結構強度。高載液深度下，液體晃蕩沖擊壓力作用于艙頂，艙頂受到的沖擊壓力作用要顯著。

(3)對于矩形板，一定柔度范圍內，相同的外部激勵，相同的長度與寬度材料屬性不變，板厚超過一定的范圍，則板厚對沖擊壓強的影響不大。因為隨著板厚的增加，板的固有頻率不斷增大，固有周期不斷減小，則結構的動態效應不顯著，不能吸收液體的沖擊能；當晃蕩載荷的作用時間超過結構固有周期，晃蕩載荷對結構的作用可以等效為靜載荷作用，這種情況下，舷側板在失穩之前，內應力σ與厚度的平方t2成反比。鋁材質模型監測點的壓強幅值要低于鋼材質模型監測點的壓強幅值。相同尺寸的情況下，2種材質的屈曲強度相差在一定的范圍，彈性模量較小的材質能夠體現出一定的彈性效應。

[1]李青.大晃動粘性流體與儲液容器的相互作用數值分析[D]. 北京：清華大學，2004.

[2]方智勇.基于Level-set法的液體晃蕩與彈性結構耦合作用研究[D].鎮江：江蘇科技大學，2006.

[3]祁江濤,顧民. LNG船液艙晃蕩的數值模擬[J].中國造船, 2007，48 (B11): 541-548.

[4]端木玉.液艙內三維液體非線性晃蕩的數值模擬[D]. 鎮江: 江蘇科技大學, 2007.

[5]沈猛,王剛,唐文勇. 基于改進VOF法的棱形液艙液體晃蕩分析[J]. 中國造船, 2009，50(1): 1-8.

[6]侯玲.液艙晃蕩與彈性防晃結構的相互耦合作用的研究[D]. 鎮江: 江蘇科技大學, 2009.

[7]韓芳,鐘冬望,蔡路軍.考慮瞬態反應影響的結構動力放大系數研究[J].力學與實踐,2013，35，(4):64-65.

2014-05-26

國家自然科學基金(10472032，50879030，51179077)江蘇高校優勢學科建設工程資助項目。

周上然(1987-)，男，碩士研究生，主要從事船舶力學研究；朱仁慶(1965-)，男，博士，教授，主要從事船舶與海洋工程載荷與響應研究。

U661.1

A