周上然，朱仁慶

(江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003)

彈性液艙液體晃蕩數(shù)值模擬

周上然，朱仁慶

(江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003)

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

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

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

1 數(shù)值方程

1.1流體運(yùn)動方程

若假設(shè)流體為不可壓縮粘性流，則在直角坐標(biāo)系下流體的連續(xù)方程和動量守恒方程分別為：

連續(xù)方程：div(v)=0

(1)

(2)

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

1.2結(jié)構(gòu)運(yùn)動的微分方程

(3)

1.3流體與結(jié)構(gòu)交界面條件

流固耦合交界面應(yīng)滿足運(yùn)動學(xué)連續(xù)條件、動力學(xué)連續(xù)條件、能量守恒原理。運(yùn)動學(xué)條件指在流固界面上，位移與速度一致，即：

ds=dfus=uf

(4)

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

σs·ns=σf·nf

(5)

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

δus·fs=δuf·ff

(6)

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

2 二維液艙晃蕩耦合數(shù)值模擬

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

(7)

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

圖1 二維矩形液艙模型

表1 模型的運(yùn)動工況

圖2 左艙壁與右艙壁監(jiān)測點處壓強(qiáng)時間歷程

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

(8)

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

ω越接近ωn，動態(tài)效應(yīng)越顯著。計算ωn時，要考慮附加質(zhì)量的影響，通常附加質(zhì)量會降低結(jié)構(gòu)的固有頻率。

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

3.1棱形液艙模型的建立

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

圖3 液艙簡化模型

3.2不同模型材質(zhì)對晃蕩壓強(qiáng)的影響

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

表2 模型的運(yùn)動工況

表3 監(jiān)測點在模型中位置

圖4 模型自由液面波形圖

表4監(jiān)測點P1至P4壓強(qiáng)結(jié)果比較

材質(zhì)P1點的壓強(qiáng)幅值P2點的壓強(qiáng)幅值P3點的壓強(qiáng)幅值P4點的壓強(qiáng)幅值鋁5600Pa4100Pa3950Pa3150Pa船用鋼6250Pa4750Pa4800Pa3900Pa

圖5 鋼質(zhì)模型監(jiān)測點(P1至P4)的壓強(qiáng)時間歷程

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

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

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

3.3不同板厚對艙室結(jié)構(gòu)應(yīng)力的影響

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

表5 監(jiān)測點P1、P3、P5的結(jié)果比較

圖6 鋁質(zhì)模型監(jiān)測點(P1至P4)的壓強(qiáng)時間歷程

圖7 模型自由液面波形圖

圖8 1.2 s時工況3、工況4、工況5模型應(yīng)力圖

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

(9)

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

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

(10)

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

圖9 工況3、工況4、工況5模型P5點應(yīng)力時間歷程

4 結(jié)語

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

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

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

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

[2]方智勇.基于Level-set法的液體晃蕩與彈性結(jié)構(gòu)耦合作用研究[D].鎮(zhèn)江：江蘇科技大學(xué)，2006.

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

[4]端木玉.液艙內(nèi)三維液體非線性晃蕩的數(shù)值模擬[D]. 鎮(zhèn)江: 江蘇科技大學(xué), 2007.

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

[6]侯玲.液艙晃蕩與彈性防晃結(jié)構(gòu)的相互耦合作用的研究[D]. 鎮(zhèn)江: 江蘇科技大學(xué), 2009.

[7]韓芳,鐘冬望,蔡路軍.考慮瞬態(tài)反應(yīng)影響的結(jié)構(gòu)動力放大系數(shù)研究[J].力學(xué)與實踐,2013，35，(4):64-65.

2014-05-26

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

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

U661.1

A