金 晶 薛鴻祥 唐文勇 張圣坤

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

0 引 言

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

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

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

1 液體晃蕩的計算模型

1.1 控制方程

對于如圖1所示的棱形液艙邊界處的部分不規(guī)則網(wǎng)格（網(wǎng)格內(nèi)有部分障礙物），可采用網(wǎng)格的體積通度來描述網(wǎng)格中可以被流體通過的體積與整個網(wǎng)格的體積之比，用網(wǎng)格的面通度λ來描述網(wǎng)格面上能夠被流體通過的面積與該面積之比?；赑CP概念，考慮體積通度λ以后，流體控制方程中的連續(xù)性方程、動量方程和體積分?jǐn)?shù)傳輸方程可分別表示如下。

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

式中，ω為相應(yīng)時刻橫搖角速度；rP為相對坐標(biāo)中質(zhì)點的位矢。

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

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

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

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

根據(jù)法向和網(wǎng)格內(nèi)F函數(shù)的大小，可將界面分為20種情況，除=0 和0 所對應(yīng)的 4種情況外，其余16種可通過對稱和翻轉(zhuǎn)運算后簡化歸并為4大類。確定了界面類型后，便可計算出單元各條邊上的流體分?jǐn)?shù)，進(jìn)而根據(jù)n時刻的速度計算出通過上邊、右邊、底邊和左邊的流體體積分?jǐn)?shù)的運輸量（以輸出為準(zhǔn)），分別用 Ft，F(xiàn)r，F(xiàn)b和Fl表示。由4種情況的流體體積分?jǐn)?shù)的運輸量，可以確定n+1時刻的流體體積分?jǐn)?shù)，從而將傳統(tǒng)的Youngs的流體體積傳輸方程修改為：

1.2 自由表面條件

在自由表面上，動力學(xué)條件須滿足表面應(yīng)力條件，忽略表面張力，則表面應(yīng)力條件為：

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

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

2 算例及對比分析

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

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

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

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

2.1 液面形態(tài)

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

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

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

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

圖6 Youngs法的液面形態(tài)（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法預(yù)報的R1，R2點沖擊壓力的最大值以及其與實驗數(shù)據(jù)相比的相對誤差如表1和表2所示，其中，30%裝載率的最大值為圖7和圖8所示的3個周期內(nèi)最大幅值的平均值，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法預(yù)報的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法預(yù)報的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法的計算結(jié)果要明顯優(yōu)于Hirt-Nichols法。從圖中可看出，30%裝載率下用兩種方法計算得到的壓力歷程曲線都呈現(xiàn)出了雙峰特征［11］，其中用PCP-Youngs法計算得到的壓力曲線的沖擊壓力時刻略微提前，其主要原因是本文為直接加載的周期性的運動激勵，而模型實驗的激勵則為由初始靜止?fàn)顟B(tài)逐漸加至穩(wěn)定周期的過程，因此，該現(xiàn)象符合實際情況。通過比較PCP-Youngs法和Hirt-Nichols法的結(jié)果發(fā)現(xiàn)，在30%裝載率時，PCP-Youngs法在R2處的沖擊壓力峰值（圖8）更接近于實驗數(shù)據(jù)，相對誤差也小一些（表1），而Hirt-Nichols法則僅在R1處所示的第2和第3個周期中更加接近于實驗數(shù)據(jù)；在46%裝載率時，PCP-Youngs法在R1和R2處的壓力峰值均較接近于實驗數(shù)據(jù)（圖9、圖10），預(yù)報的壓力曲線也相對較光順、連續(xù)?？梢妼τ谂撌疫吔鐗毫Φ念A(yù)報，從總體上而言，PCP-Youngs法要優(yōu)于Hirt-Nichols法。

2.3 速度和加速度

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

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

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

對速度和加速度的時間歷程曲線進(jìn)行兩次中值、兩次均值濾波，并將速度和加速度曲線與濾波后的曲線相減，得到速度和加速度曲線的數(shù)值波動曲線如圖13和圖14所示。表3所示為速度和加速度波動曲線的標(biāo)準(zhǔn)差。由表3可看出，PCP-Youngs法的數(shù)值波動約減小到了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 速度和加速度曲線波動的標(biāo)準(zhǔn)差Tab.3 The standard deviation of velocity wave and acceleration wave

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

3 網(wǎng)格敏感性分析

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

本文的方法可用于大型補給艦船液艙結(jié)構(gòu)晃蕩設(shè)計載荷的快速、準(zhǔn)確預(yù)報。

［1］丁仕風(fēng)，唐文勇，張圣坤.速度突變情況下液化天然氣船液艙內(nèi)晃蕩問題的仿真［J］.上海交通大學(xué)學(xué)報，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數(shù)值模擬［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，F(xiàn)EKEN 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］沈猛，王剛，唐文勇.基于改進(jìn)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］劉儒勛，王志峰.數(shù)值模擬方法和運動界面追蹤［M］.合肥：中國科學(xué)技術(shù)大學(xué)出版社，2001.

［6］端木玉，朱仁慶.基于Youngs法的液艙晃蕩大幅晃蕩數(shù)值模擬［J］.船舶力學(xué)，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］朱仁慶.液體晃蕩及其與結(jié)構(gòu)的相互作用［D］.無錫：中國船舶科學(xué)研究中心，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］.江蘇科技大學(xué)學(xué)報：自然科學(xué)版，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.