大型補給艦船液艙晃蕩載荷計算方法

金 晶 薛鴻祥 唐文勇 張圣坤

上海交通大學海洋工程國家重點實驗室，上海 200240

0 引 言

補給艦船用于在海上航行或停泊中為已方艦艇提供燃料、淡水、食品等消耗品，以及魚雷、水雷、炮彈、導彈等武器，作戰艦艇通過航行補給不僅能延長活動半徑，提高在航率，增加海上活動時間，而且還避免了對固定基地的過分依賴。其中，大型綜合補給艦船還可為遠洋艦隊以及航母編隊提供強大的支援，受到各海軍強國的重視。隨著我國海軍戰略發展的需要，大型綜合補給艦船的研究與建造也勢在必行。

隨著補給艦船噸位的增大，裝載燃油等補給的液艙尺度也隨之增大，在某些裝載情況下，有可能會發生劇烈的液體晃蕩諧振運動，引起沖擊載荷，從而危及艙壁及艙室內部構件結構的安全。因此，在大型補給艦船的液艙結構設計中，應考慮可能產生的液體晃蕩沖擊載荷。液艙晃蕩沖擊載荷預報是當前的研究熱點，并已取得一定的成果，采用Fluent等商用軟件，可以計算液艙內流體自由液面的變化、流體質點的速度矢量和流場內的壓強分布等。例如，丁仕風等［1］基于 Fluent，采用VOF法對速度突變情況下液化天然氣船液艙內的晃蕩進行了模擬和分析；張書誼等［2］采用Fluent對矩形液艙的橫蕩流體載荷進行了數值模擬并和實驗結果進行了對比。但是，在初始設計階段或審圖階段，熟練掌握商用CFD軟件對于結構設計人員來說并不方便，因此，發展可靠且方便于工程應用的晃蕩載荷計算方法，并結合試驗數據進行修正仍是當前需重點研究的方向。

大型補給艦船的液艙截面可能為非矩形截面，邊界復雜，也可能含有內部構件，因此，其晃蕩載荷的計算方法應能處理非矩形邊界，并能準確、穩定地預報液艙內部的速度和加速度時間歷程，從而準確獲得艙壁以及構件上的晃蕩載荷。對于非矩形液艙，如棱形液艙的晃蕩模擬，傳統的VOF方法采用的是將矩形網格的折線近似處理為斜線，從而導致模擬不夠準確，引起一定的誤差［3］。沈猛等［4］通過引入部分單元參數（Partial Cell Parameter，PCP）的概念對VOF法進行了改進，用于處理棱形液艙，但自由表面重構采用的還是傳統的Hirt-Nichols算法。由于Hirt-Nichols算法為零階精度［5］，對晃蕩自由表面的模擬精度較低，從而導致速度、加速度預報結果離散度較大，給液艙內部構件的晃蕩載荷分析帶來了困難。端木玉等［6］采用Youngs方法對矩形液艙大幅晃蕩進行模擬，得出用Youngs法比用Hirt-Nichols法能更準確地處理液艙大幅晃蕩的結論，但該文的算法適于處理矩形截面的液艙晃蕩計算，對于棱形液艙等復雜邊界的情況還有待進一步的完善。本文將以棱形液艙為對象，采用改進的VOF法，基于精度較高的Youngs法重構表面，并結合部分單元參數概念［7］（以下簡稱“PCP-Youngs法”）模擬晃蕩運動，以有效改善自由表面的模擬效果，特別是提高艙內速度和加速度的預報精度以及數值穩定性，為在大型補給艦船液艙結構設計中加入對于晃蕩載荷沖擊影響的考慮奠定基礎。

1 液體晃蕩的計算模型

1.1 控制方程

對于如圖1所示的棱形液艙邊界處的部分不規則網格（網格內有部分障礙物），可采用網格的體積通度來描述網格中可以被流體通過的體積與整個網格的體積之比，用網格的面通度λ來描述網格面上能夠被流體通過的面積與該面積之比。基于PCP概念，考慮體積通度λ以后，流體控制方程中的連續性方程、動量方程和體積分數傳輸方程可分別表示如下。

式（1）～式（3）中，u為流體速度矢量；p為流體壓力；F為流體的體積分數；ρ為流體密度；μ為運動粘性系數，f為流體微元體所受到的體積力，該力是引起流體晃蕩的主要激勵，采用相對坐標系，表達式為：

式中，ω為相應時刻橫搖角速度；rP為相對坐標中質點的位矢。

圖1 部分單元參數的變量位置與定義Fig.1 The variable distribution and definition in PCP method

控制方程的離散采用有限差分法，差分網格采用交錯形式。其中，網格的體積通度λ、壓力P和流體體積分數F定義于網格中心，而網格的速度u，v和網格的上、右、下、左4個面的通度 λt，λr，λb，λl定義于網格的4條邊上（圖1）。

考慮體積通度和面通度后的連續方程具體可參見文獻［8］。對于體積傳輸方程，本文采用了精度較高的Youngs法［9］。

Youngs法基于幾何學原理，在單個單元網格內是用斜線段近似模擬界面，其可以分兩種情況進行計算：一是計算單元為滿網格，F=1；二是計算單元為自由表面網格，0＜F＜1。其中，對于第2種情況涉及界面重構的問題，需要先計算網格內界面的法向n=，其中

根據法向和網格內F函數的大小，可將界面分為20種情況，除=0 和0 所對應的 4種情況外，其余16種可通過對稱和翻轉運算后簡化歸并為4大類。確定了界面類型后，便可計算出單元各條邊上的流體分數，進而根據n時刻的速度計算出通過上邊、右邊、底邊和左邊的流體體積分數的運輸量（以輸出為準），分別用 Ft，Fr，Fb和Fl表示。由4種情況的流體體積分數的運輸量，可以確定n+1時刻的流體體積分數，從而將傳統的Youngs的流體體積傳輸方程修改為：

1.2 自由表面條件

在自由表面上，動力學條件須滿足表面應力條件，忽略表面張力，則表面應力條件為：

式中，un為垂直于自由表面的法向速度；ut為切向速度；p0為艙室內氣體的壓力。

對于自由表面網格速度的求解，采用常數外插與連續性條件和自由表面條件相結合的方法。對于自由表面邊界速度，采用常數外插，但當判斷出網格轉換為不含自由表面邊界時，則要求網格滿足連續性條件，盡量避免壓力和速度的突變。

2 算例及對比分析

為了驗證基于PCP的Youngs法（PCP-Youngs法）對棱形液艙晃蕩數值模擬的可行性和計算精度，與文獻［10］中的棱形液艙模型晃蕩試驗數據進行了對比驗證，同時，還將本文的結果與基于PCP的Hirt-Nichols法的計算結果進行了對比。

選取文獻［10］中的30%和46%兩種裝載水平進行晃蕩計算，兩種裝載工況的橫搖激勵周期分別為1.517 s和1.207 s，艙室模型的橫搖幅值均為5.73°。

棱形液艙模型尺寸及各測試點的位置如圖2所示。圖中R1～R4分別為艙室邊界的壓力測試點，P1～P3，Q1～Q3分別為30%和 46%裝載率下液面高度、2/3液面高度和1/3液面高度處的監測點，這些監測點位于液艙的中間縱向剖面處，用于監測速度和加速度。

圖2 棱形液艙模型Fig.2 Prismatic liquid tank model

2.1 液面形態

裝載率為30%和46%，網格為114×78時的液體波面形態分別如圖3～圖6所示（其中圖3角隅處出現的毛刺是由于軟件畫圖設置方面的因素，因軟件在畫圖過程中采用的是簡化的矩形的網格劃分，并在整個網格內部按照網格的壓力填充了顏色）。由圖中可看出，采用Hirt-Nichols法會出現非物理性的飛濺，且自由表面也不連續，而采用Youngs法則未出現這些異常，自由液面相對光順、連續。

圖3 Hirt-Nichols法的液面形態（30%裝載率）Fig.3 Liquid surface form of Hirt-Nichols method（30%filling ratio）

圖4 Youngs法的液面形態（30%裝載率）Fig.4 Liquid surface form of Youngs method（30%filling ratio）

圖5 Hirt-Nichols法的液面形態（46%裝載率）Fig.5 Liquid surface form of Hirt-Nichols method（46%filling ratio）

圖6 Youngs法的液面形態（46%裝載率）Fig.6 Liquid surface form of Youngs method（46%filling ratio）

2.2 艙室邊界壓力

圖7和圖8所示分別為30%裝載率時測試點R1和R2的壓力歷程曲線對比，圖9和圖10所示分別為46%裝載率時測試點R1和R2的壓力歷程曲線對比。30%和46%裝載率下PCP-Youngs法和Hirt-Nichols法預報的R1，R2點沖擊壓力的最大值以及其與實驗數據相比的相對誤差如表1和表2所示，其中，30%裝載率的最大值為圖7和圖8所示的3個周期內最大幅值的平均值，46%裝載率與之類似。

圖7 R1處壓力（30%裝載率）Fig.7 Pressure of R1（30%filling ratio）

圖8 R2處壓力（30%裝載率）Fig.8 Pressure of R2（30%filling ratio）

圖9 R1處壓力（46%裝載率）Fig.9 Pressure of R1（46%filling ratio）

圖10 R2處壓力（46%裝載率）Fig.10 Pressure of R2（46%filling ratio）

表1 30%裝載率下PCP-Youngs法和Hirt-Nichols法預報的R1，R2點處的沖擊壓力和相對誤差Tab.1 The impact pressure and proportional error on R1and R2with PCP-Youngs method and H-N method（30%filling ratio）

表2 46%率下PCP-Youngs法和Hirt-Nichols法預報的R1，R2點處的沖擊壓力和相對誤差Tab.2 The impact pressure and proportional error on R1and R2with PCP-Youngs method and H-N method（46%filling ratio）

由表中可見，PCP-Youngs法的計算結果要明顯優于Hirt-Nichols法。從圖中可看出，30%裝載率下用兩種方法計算得到的壓力歷程曲線都呈現出了雙峰特征［11］，其中用PCP-Youngs法計算得到的壓力曲線的沖擊壓力時刻略微提前，其主要原因是本文為直接加載的周期性的運動激勵，而模型實驗的激勵則為由初始靜止狀態逐漸加至穩定周期的過程，因此，該現象符合實際情況。通過比較PCP-Youngs法和Hirt-Nichols法的結果發現，在30%裝載率時，PCP-Youngs法在R2處的沖擊壓力峰值（圖8）更接近于實驗數據，相對誤差也小一些（表1），而Hirt-Nichols法則僅在R1處所示的第2和第3個周期中更加接近于實驗數據；在46%裝載率時，PCP-Youngs法在R1和R2處的壓力峰值均較接近于實驗數據（圖9、圖10），預報的壓力曲線也相對較光順、連續。可見對于艙室邊界壓力的預報，從總體上而言，PCP-Youngs法要優于Hirt-Nichols法。

2.3 速度和加速度

圖11和圖12所示為30%裝載率情況下，中間縱向剖面上Q2點處的速度和加速度的時間歷程曲線對比。46%裝載率下以及其他各測試點處的速度和加速度歷程對比曲線均與之類似。由圖中可看出，采用PCP-Youngs法得到的速度和加速度曲線更加光順，擾動較小，而用Hirt-Nichols法得到的加速度曲線的穩定性則較差，存在較多震蕩。這是由于Hirt-Nichols算法的Donor-Acceptor型零階方法對晃蕩自由表面模擬的精度太低，從而直接導致自由表面網格流體的傳輸計算存在一定的誤差，且隨著計算中誤差的累積，內部流場速度場，特別是加速度場的數值模擬會存在較多的誤差與擾動。

圖11 Q2處水平速度（30%裝載率）Fig.11 Velocity of Q2（filling ratio 30%）

圖12 Q2處水平加速度（30%裝載率）Fig.12 Acceleration of Q2（filling ratio 30%）

對速度和加速度的時間歷程曲線進行兩次中值、兩次均值濾波，并將速度和加速度曲線與濾波后的曲線相減，得到速度和加速度曲線的數值波動曲線如圖13和圖14所示。表3所示為速度和加速度波動曲線的標準差。由表3可看出，PCP-Youngs法的數值波動約減小到了Hirt-Nichols法的1/3～1/4。

圖13 Q2處水平速度波動（30%裝載率）Fig.13 Velocity wave of Q2（filling ratio 30%）

圖14 Q2處水平加速度波動（30%裝載率）Fig.14 Acceleration wave of Q2（filling ratio 30%）

表3 速度和加速度曲線波動的標準差Tab.3 The standard deviation of velocity wave and acceleration wave

通過上述對比可以看出，與Hirt-Nichols法相比，PCP-Youngs法在自由液面、艙壁壓力、艙室內部的速度和加速度預報等方面均有不同程度的改善，特別是有效改善了加速度預報的數值穩定性，為更加準確地預報液艙內部構件的晃蕩載荷提供了有力保障。

3 網格敏感性分析

為了驗證本文方法的可靠性和穩定性，選擇30%裝載率工況進行了網格敏感性分析，設定了3種網格大小情況：57×39，114×78和171×117。

選取底邊艙斜板上R1處的壓力進行對比分析，計算結果如圖15和圖16所示。由圖可見，隨著網格尺度的加密，沖擊峰值時刻也略微提前，但相差不大。選取液艙中間縱向剖面上Q2處進行速度和加速度對比，其結果如圖17和圖18所示。由圖可見，隨著網格的細化，速度和加速度的歷程曲線基本一致，只是在峰值處有細微變化。表4所示為不同網格密度下加速度和速度的峰值。由表可見，對于用不同網格密度計算出的速度和加速度，僅相差約1%～3%。由此可見，網格的粗細對于PCP-Youngs法的計算結果影響較小，特別是對速度和加速度的影響更小，其在實際工程計算中不需要采用很精細的網格便可得到較準確的數值模擬結果。

圖15 R1處壓力網格敏感性分析Fig.15 The pressure cell sensibility analysis of R1

圖16 R1處壓力網格敏感性分析局部放大圖Fig.16 Zoomed picture in pressure peak of R1

圖17 Q2處速度網格敏感性分析Fig.17 Cell sensibility analysis for velocity of Q2

圖18 Q2處加速度網格敏感性分析Fig.18 Cell sensibility analysis for acceleration of Q2

表4 不同網格密度下速度和加速度的峰值Tab.4 The maximum and minimum value of velocity and acceleration in different cell numbers

本文采用 57×39，114×78和 171×117這 3種網格，用PCP-Youngs法對30%的裝載工況進行了15個周期的計算，所耗費的時間分別約為8，40和150 min，可節省的計算資源非常可觀。因此，在實際工程應用中，只需采用網格尺度為液艙邊界長度2%左右的較粗網格，約10 min即可快速、準確地計算出液艙內部構件位置處的速度和加速度，從而快速確定液艙內部構件的晃蕩載荷。

4 結 論

本文針對大型補給艦船液艙結構的晃蕩沖擊問題，考慮液艙為非矩形截面的復雜邊界，以及可能含有內部構件的情況，基于改進的VOF法，采用Youngs法重構自由液面，并結合PCP概念，提出了一種精度更高且數值穩定的處理非矩形復雜艙室邊界晃蕩載荷的計算方法。分別對不同裝載情況下的艙室邊界壓力、液艙中間縱向剖面處的速度、加速度以及網格敏感性等重要參數進行了對比分析，得到以下結論：

1）通過與模型試驗對比，顯示本文建立的PCP-Youngs法可有效改進自由表面的模擬效果；通過與實驗數據進行對比，發現與Hirt-Nichols法相比，用PCP-Youngs法得到的壓力預報相對誤差可減小約2%～12%。

2）與Hirt-Nichols法相比，本文建立的PCPYoungs法不僅具有很好的數值穩定性，速度和加速度曲線有較強的周期性且精度更高，而且數值波動也減小到了Hirt-Nichols法的1/3～1/4左右。

3）網格敏感性分析顯示，當網格尺度不大于液艙邊界長度的2%左右時，數值計算的網格密度對計算結果影響較小，特別是對速度和加速度的時間歷程曲線的影響更小，曲線幾乎一致，不同網格之間的誤差約為1%～3%，在實際工程中的誤差許可范圍之內。因此，在實際工程應用中不必采用極細的網格便可得到較滿意的預報結果。

本文的方法可用于大型補給艦船液艙結構晃蕩設計載荷的快速、準確預報。

［1］丁仕風，唐文勇，張圣坤.速度突變情況下液化天然氣船液艙內晃蕩問題的仿真［J］.上海交通大學學報，2008，42（6）：919-923.DING S F，TANG W Y，ZHANG S K.Simulation of liquid sloshing in the cabin caused by liquefied natural gas ship’s variable speed［J］.Journal of Shanghai Jiao Tong University，2008，42（6）：919-923.

［2］張書誼，段文洋.矩形液艙橫蕩流體載荷的Fluent數值模擬［J］.中國艦船研究，2011，6（5）：73-77.ZHANG S Y，DUAN W Y.Numerical simulation of sloshing loads on rectangular tank based on fluent［J］.Chinese Journal of Ship Research，2011，6 （5）：73-77.

［3］KKEEFSMAN K M T，FEKEN G，VELDMAN A E P，et al.A volume of fluid based simulation method for wave impact problems［J］.Journal of computational Physics，2005，206（1）：363-393.

［4］沈猛，王剛，唐文勇.基于改進VOF法的棱形液艙液體晃蕩分析［J］.中國造船，2009，50（1）：1-9.SHEN M，WANG G，TANG W Y.Liguid sloshing analysis on prismatic tanks based on improved VOF method［J］.Shipbuilding of China，2009，50（1）：1-9.

［5］劉儒勛，王志峰.數值模擬方法和運動界面追蹤［M］.合肥：中國科學技術大學出版社，2001.

［6］端木玉，朱仁慶.基于Youngs法的液艙晃蕩大幅晃蕩數值模擬［J］.船舶力學，2009，13（1）：9-18.DUAN M Y，ZHU R Q.Numerical simulation of violent liquid sloshing in a tank based on Youngs method［J］.Journal of Ship Mechanics，2009，13（1）：9-18.

［7］LIN P Z.A fixed-grid model for simulation of a moving body in free surface flows［J］.Computers and Fluids，2007，36（3）：549-561.

［8］朱仁慶.液體晃蕩及其與結構的相互作用［D］.無錫：中國船舶科學研究中心，2002.ZHU R Q.Time domain simulation of liquid sloshing and its interaction with flexible structure［D］.Wuxi：China Ship Scientific Research Center，2002.

［9］端木玉，朱仁慶.流體體積方程的求解方法［J］.江蘇科技大學學報：自然科學版，2007，21（2）：10-15.DUAN M Y，ZHU R Q.Method for solving fluid volume equation［J］.Journal of Jiangsu University of Science and Technology：Natural Science Edition，2007，21（2）：10-15.

［10］MIKELIS N E，MILLETR J K，TAYLOR K V.Sloshing in partially filled liquid tank and its effect on ship motions：numerical simulations and experimental verification［J］.The Royal Institution of Naval Architects Transaction，1984，126：267-281.

［11］KIM Y，SHIN Y S，LEE K H.Numerical study on slosh-induced impact pressures on three-dimensional prismatic tanks［J］.Applied Ocean Research，2004，26（5）：213-226.