三維菱形液艙劇烈晃蕩和共振頻率數值研究

辛建建, 方 田, 石伏龍

(1. 寧波大學 海運學院， 浙江 寧波 315211； 2. 中國船舶科學研究中心， 江蘇 無錫 214082；3. 武漢理工大學 交通學院， 武漢 430063)

風浪中航行的儲液船舶，如液化石油氣(LPG)、液化天然氣(LNG)船，在受到外部激勵時，液體會對容器壁產生強烈的砰擊作用，造成結構的變形甚至破壞，對船舶穩性和結構強度造成嚴重威脅.液體發生晃蕩時會伴隨產生波面破碎、波對艙壁的高速砰擊、液滴形成及氣泡卷入等強非線性現象.而且，晃蕩引起的砰擊載荷會對艙壁產生強烈的局部沖擊，影響船舶的穩定性，甚至破壞液艙的圍欄系統，導致沉船事故.因此，研究晃蕩現象和預報晃蕩砰擊載荷對液艙結構設計和船舶的安全運營至關重要.

在過去，模型實驗[1]和解析半解析方法[2]是船舶液艙晃蕩研究的常用手段.然而，模型實驗費用昂貴、條件苛刻(如低溫)，解析方法局限于簡單幾何形狀液艙中的晃蕩問題. 隨著計算機技術的飛速發展，數值方法處理晃蕩問題顯示出較大的優越性，可以克服解析方法無法解決的困難，諸如液面大變形、搖蕩運動影響等.與費用昂貴、時間較長的模型實驗相比，數值模擬方法的成本低廉、靈活方便.從自由表面跟蹤技術分，有流體體積分數(VOF)法[3]、水平集(Level Set)法[4]和光滑粒子流體動力學(SPH)法[5]等.根據計算網格可分為貼體網格法[6]、直角網格法[7]和無網格法[5].

在絕大多數的晃蕩流模擬中，對象為矩形液艙，邊界條件能直接施加，因為計算網格與液艙邊界重合.然而在實際工程應用中，液艙可能是任意形狀如菱形和橢圓形，另外液艙中也有可能存在內部結構如橫隔板.為了處理非矩形液艙或具有內部結構的情況，一個常用方法是貼體網格法[6].例如，Jiang等[8]基于 OpenFOAM模擬了菱形液艙的晃蕩問題，Zhao等[9]基于重疊網格法模擬了部分充液 LNG 液艙的三維晃蕩流.然而，對于多自由度激勵、三維復雜形狀液艙和多個內部結構存在的情況，貼體網格方法面臨挑戰.為了處理這些問題，直角網格方法是一個具有吸引力的替代方法.在直角網格方法中，液艙邊界并不與背景網格線重合，通過引入力源項或重構浸入邊界的插值模板以施加邊界條件.為了表示不規則邊界或內部結構物對晃蕩流的影響，Kim等[10]引入一個緩沖層以施加液艙斜邊界附近的自由表面邊界條件.Lee等[11]通過引入一個二次多項式內插格式以施加液艙不規則邊界的無滑移邊界條件，根據不規則邊界與網格線之間的相對位置關系，內插模板分為 13 種情況.

目前，三維復雜液艙劇烈晃蕩的數值預報仍面臨巨大的挑戰，文獻中也鮮見對菱形液艙共振晃蕩頻率的研究.為此，本文采用一個基于直角網格的三維多相流模型模擬了菱形液艙劇烈晃蕩問題，并預報了菱形液艙的共振頻率，重點專注于艙壁沖擊壓力和全局晃蕩流運動，忽略局部物理現象湍流和可壓縮的影響.在規則網格上求解不可壓縮黏性N-S方程，在長方體四角布置4個楔形體以表示菱形斜邊，采用徑向基函數虛擬網格法[12]，通過在固體區域內部布置合適的虛擬網格，以考慮不規則邊界對流場的影響.

另外，采用一個三維梯度增量Level Set(GALS)兩相流模型[13]捕捉強非線性自由表面.基于所開發的模型，模擬了菱形液艙的大幅晃蕩現象，并進行了數值收斂性驗證.進一步研究了不同充液水深下壁面沖擊壓力隨激勵頻率的變化規律，確定了菱形液艙的共振頻率，以期為LNG船舶的液艙結構設計提供理論和技術指導.

1 數學模型

不可壓縮N-S方程由描述任意流場的質量和動量守恒方程組成.在三維直角坐標系下，非定常黏性不可壓縮N-S方程表達如下：

(1)

(2)

式中：速度向量u=[uvw]，u、v、w分別是直角網格x、y和z方向上的速度分量；t為時間；ρ為密度；p為壓力；υ為流體的運動黏性系數.對于外部激勵力fa，其在多數自由表面流動中為重力加速度g.然而，對于晃蕩問題，外部力包括重力和晃蕩運動引起的慣性力以考慮液艙壁面運動對艙內流場的影響，其表達式為[7]

(3)

式中：V和Ω分別為非慣性坐標系下的平移和旋轉速度；r和R分別為網格節點坐標和旋轉中心.式(3)右邊5項分別是非慣性坐標系的重力加速度、平移、旋轉、離心及科氏加速度.

對于自由表面的捕捉，通過求解GALS方程實時更新距離函數φ和梯度ψ以隱式捕捉兩相界面.GALS方程[14]表示為

(4)

(5)

在兩相流模式計算中，根據自由表面的位置變化，定義在計算網格單元中心密度ρ及黏性系數υ的計算公式為

ρ=Hε(φ)ρW+(1-Hε(φ))ρA

(6)

υ=Hε(φ)υW+(1-Hε(φ))υA

(7)

式中：Hε(φ)為Heaviside函數以區分流體和空氣界面的流體屬性；ρ及υ的下標W和A分別表示水和空氣.為了減緩因密度(或黏性)介質的跳躍對流動的影響，使得過渡層內物性參數光滑變化，抑制數值振蕩.

2 數值求解方法

2.1 N-S方程求解器

在本文研究中，基于自主開發的時間半隱式有限差分法，在交錯直角網格上求解非定常不可壓縮N-S方程.采用分步法解耦速度和壓力，TVD-RK2(Total Variation Diminishing Second-order Runge-Kutta)格式進行時間推進，對流項采用時間顯式處理，黏性項采用Crank-Nicolson格式半隱式處理.對于空間離散，對流項采用高階TVD-MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)格式離散，黏性項以中心差分格式離散.另外，對于壓力泊松方程離散形成的大型稀疏線性方程組，采用預處理的BICGSTAB(Bi-conjugate Gradients Stabilize)方法結合一種對角存儲格式求解.更多的細節請參考文獻[12].

2.2 GALS方法數值求解

GALS方程式(4)、(5)以特征形式表示為

(8)

求解過程分為兩步：① 采用Shu-Osher RK3方法從t+Δt到t進行向后時間積分求解特征曲線方程?x(τ)/?τ=u(x(τ),τ)；② 然后沿著特征曲線采用顯式梯度格式進行向前積分求解式(8)以更新φ和ψ.值得注意的是在任意位置的變量φ和ψ通過Hermite立方多項式插值得到以獲得三階精度，另外，通過對Hermite立方多項式格式進行微分以解析方式計算φ和ψ的一階和二階導數.詳細求解過程參考文獻[13].

2.3 不規則邊界處理

浸入邊界法的網格生成過程簡單，但由于浸入邊界本身并未包含在計算網格中，精確表示背景網格中浸入邊界的存在是個難點.本文采用基于徑向基函數(RBF)的隱式等值面表示方法[12]，光滑基函數擬合離散的物面控制點以構造任意復雜物體的等值曲面，并且根據等值面距離函數的正負號識別場點屬性.下一步就可重構虛擬網格的流體變量值以施加邊界條件.通過在浸入物體內部布置有限數量的虛擬網格以表示靜止或運動邊界的存在，以適當的插值模塊及技術(如雙線性插值方法)，施加浸入動邊界條件.

虛擬網格的插值重構過程可細分為3步：標記鏡像點、插值重構和滿足鏡像條件.首先，對于每個虛擬網格(GC)，沿著物體表面的外法向從虛擬網格單元中心延伸到流體域以確定唯一的一個鏡像點(MP)，如圖1所示.然后采用三線性插值格式計算每個鏡像點的流體變量.然而，存在一個或更多的模板點是虛擬網格本身的情況，為了避免繁瑣的內插模板分類過程，首先對封閉鏡像點(MP)的8個鄰近網格節點進行判斷.如果該網格節點為流體網格點(FC)，則將它確定為模板點，1、2、3、4即為二維插值格式中識別出的模板點，如圖1(a)所示.如果為虛擬網格點(GC)，則相應的邊界點(BP)作為模板點，而不是虛擬網格點本身，如圖1(b)所示.如果更多的鄰近網格點識別為虛擬網格點，同樣，相應的邊界點BP作為模板點(見圖1(c)).在求得鏡像點的流體變量后，虛擬網格點的流體變量可通過滿足物面的鏡像條件直接得到.當通過標準的七點中心差分格式在整個計算域求解N-S方程后，物面邊界條件就會得到隱式施加.更多的求解細節見文獻[12].

圖1 虛擬網格法的二維插值格式Fig.1 Two-dimensional interpolation format of virtual grid method

3 三維菱形液艙晃蕩分析

為研究復雜液艙系統下的晃蕩特性，數值模擬了不同水深和激勵頻率下橫搖激勵的三維菱形液艙晃蕩問題.根據Mikelis等[15]對不同充液水平下棱形液艙的晃蕩實驗研究，液艙模型的長、寬、高分別為L=0.741 m、W=0.399 m、H=0.405 m，其幾何模型如圖2所示.液艙遭受強迫橫搖激勵，橫搖中心Y1、Y2在x軸中心位置距離液艙底部0.194 m.液艙橫搖激勵的運動描述為θy=θ0sinωt.θy為繞橫搖中心Y1、Y2的橫搖角；橫搖幅值θ0=0.1 rad；ω為激勵頻率，ω=2π/T，T為橫搖周期.在右壁面寬度方向中心位置設置7個壓力傳感器R1～R7，距離液艙底部的距離分別為0.075，0.113 5，0.205，0.296 5，0.341，0.390，0.405 m，R7距離右壁面0.02 m.在寬度中心距離左邊界 0.135 6 m位置設置浪高儀WR.在液艙所有邊界設置速度無滑移邊界條件，在上壁面設置壓力為0，其余邊界設置Neumann邊界條件.由于采用規則的矩形計算區域，為了考慮菱形液艙的斜邊界對流場的影響，在液艙橫剖面的四角布置4個楔形體，采用虛擬網格法施加楔形體表面的物面邊界條件.

圖2 菱形液艙幾何模型和傳感器布置(m)Fig.2 Geometrical model of prismatic tank and placement of sensors (m)

3.1 數值驗證

首先模擬了充液比h/H=0.9 (h為水深)和周期T=0.97 s的工況以驗證本文多相流方法的精度和可靠性.分別采用80×50×60、100×60×80和120×80×100的3套網格離散計算區域以研究網格收斂性，根據CFL(Courant-Friedreichs-Lewy)條件決定變時間步長，其中時間松弛系數設置為ωCFL=0.3.本文方法3套網格下計算的在R3點的壓力和WR位置的波浪爬高與實驗結果進行了比較，如圖3所示.圖中η為波浪爬高.在3套網格下的壓力和波浪幅值與實驗數據吻合較好，盡管觀察到峰值壓力有略微差別.然而，隨著網格細化，峰值壓力逐漸趨近實驗數據，如圖3(a)所示，驗證了本文方法良好的網格收斂性.另外，波峰在η=0.4 m的位置發生截斷，因為波浪爬升到了液艙頂部.為了研究時間步長的收斂性，選擇3個不同的時間松弛系數ωCFL=0.1、ωCFL=0.3和ωCFL=0.6，分別在x、y和z方向選用100×60×80的計算網格.圖4展示了3個時間松弛系數下計算的壓力和波浪爬升與實驗數據的比較結果.總體而言，在3個時間松弛系數下本文計算結果與實驗數據吻合.隨著時間步長的縮小，壓力峰值會更加尖銳，更趨向于實驗結果，ωCFL=0.1和ωCFL=0.6時的壓力峰值相比實驗數據分別略小4.8%和7.1%，偏差在合理范圍內，再一次驗證了本文方法的精度和不同時間步長下的可靠性.為了在求解質量和計算成本間取得平衡， 在其他工況計算中時間松弛系數選為ωCFL=0.3，網格設置為100×60×80.

圖3 h/H=0.9時三套網格下本文方法和實驗得到的壓力和波浪爬高時間歷程Fig.3 Time history of pressure and wave elevation by present method on three grids and experiment at h/H=0.9

圖4 h/H=0.9時三個時間松弛系數下本文方法和實驗得到的壓力和波浪爬高時間歷程Fig.4 Time history of pressure and wave elevation by present method on three time coefficients and experimental data at h/H=0.9

分別模擬了h/H=0.3、T=1.517 s和h/H=0.61、T=1.112 s的兩個工況以進一步驗證本文方法的可靠性.圖5、6分別為兩個工況下波高和3個測量點壓力與Mikelis等實驗數據的比較結果.當h/H=0.3時，在垂向壁面R2點的壓力脈沖值略高于實驗結果，自由表面曲線在波谷與實驗值吻合，但波峰略高于實驗值(見圖5).當水深繼續增加到h/H=0.61時，壓力變化趨勢和形態與實驗結果吻合，自由表面曲線的波高略大于實驗值(見圖6).本文結果與實驗值的偏差可能是由于湍流或空氣可壓縮的影響所致，考慮到三維晃蕩的復雜性、非線性和瞬態性，偏差在合理范圍內.從圖5可看出，晃蕩波周期性地砰擊壁面產生脈沖形壓力曲線.同時，在每個周期內出現連續的雙波峰現象，其由強烈的非線性波浪砰擊產生.當液艙運動到一側時，流體被高速運動液艙驅動到同一側，快速運動的水體砰擊垂向壁面，導致第1個壓力峰值產生，盡管液艙達到最大橫搖角，但水體在慣性作用下繼續向前運動，導致第2個壓力峰值產生.當水深增加到h/H=0.61時，由于此時晃蕩流的非線性程度降低，壓力曲線呈現周期性的單壓力峰值現象，自由表面曲線呈現較為規則的正弦曲線變化趨勢，R2點的壓力呈現周期性的簡諧振蕩特性，均衡壓力在 1 250 Pa左右，該點在靜水面以下，砰擊壓力由靜水壓力和動壓共同主導.

圖6 h/H=0.61時的的壓力和自由表面時間歷程Fig.6 Time history of pressure and free surface at h/H=0.61

圖7、8分別為兩個工況下在兩個瞬時時刻的三維波面圖，其中4個藍色楔形體表示菱形液艙的幾何外形，用來考慮不規則邊界與流場的相互作用.h/H=0.3時，整個水體隨著液艙搖蕩向一側產生劇烈運動，水體上層流體速度大于底部流體，形成大的涌浪，如圖7(a)所示.當液艙達到最大橫搖角時，涌浪砰擊壁面，隨后，液艙向反方向運動，驅動水體向反方向運動，并形成相同的涌浪，如圖7(b)所示.h/H=0.61時，水體跟隨著液艙運動并一直爬升到液艙頂部斜邊，生成陡峭的、與上斜邊貼合的非線性波浪.當波浪爬升到斜邊最大位置時，其在重力作用下下落，形成波浪翻卷，同時，后方的波浪繼續向前運動，形成皺褶形的波浪紋，如圖8所示.

圖7 h/H=0.3時兩時刻的瞬時自由表面Fig.7 Instantaneous free surface at two moments at h/H=0.3

圖8 h/H=0.61時兩時刻的瞬時自由表面Fig.8 Instantaneous free surface at two moments at h/H=0.61

3.2 共振頻率預報

在載液船舶液艙設計中，液艙共振頻率(共振頻率≈自然頻率)預報對于避免危險工況的發生、保障航行安全具有重要的意義.對于矩形液艙，液艙自然頻率可解析計算為

(9)

m,n=0, 1, 2, …

圖9 4個充液水深下艙壁壓力幅值與激勵頻率之間的關系Fig.9 Pressure amplitude at tank walls versus excitation frequency at four filling water levels

然而，對于不規則形狀液艙如菱形液艙，其自然頻率無法解析計算得到.鑒于此，本文數值研究了不同水深下液艙激勵頻率與艙壁局部位置壓力幅值之間的關系，進一步確定了不同水深下菱形液艙的共振頻率.圖9展示了4個充液水深下艙壁壓力測量點上的壓力幅值隨液艙激勵頻率的變化趨勢，壓力幅值取穩定狀態時5個周期壓力峰值的平均值，以避免初始擾動和瞬時波動的影響.可以看出，隨著激勵頻率的增加，艙壁各點處的壓力幅值迅速增加到最大值，當激勵頻率繼續增加時，壓力幅值迅速減小，而且壓力幅值減小的速度大于增加的速度.壓力幅值最大時對應的激勵頻率即為共振頻率，也就是自然頻率.由此可得到不同水深下菱形液艙的自然頻率，如表1所示.為了提供參考，計算了相同主尺度下矩形液艙的最低階自然頻率，此時m=1,n=0.由表1可看出，相同水深下菱形液艙的最低階自然頻率小于相同主尺度下矩形液艙的最低階自然頻率，原因可能是液艙與晃蕩流體強烈的非線性相互作用.應該指出的是，本文預報激勵頻率的精度為小數點后一位，因為采用的最小激勵頻率間隔是0.1 rad/s，為了提高預報精度，需要取更小的激勵頻率步長.

表1 4個充液水深下矩形液艙和菱形液艙的自然頻率Tab.1 Natural frequency of rectangular and prismatic tanks at four filling water levels

隨著水深的增加，自然頻率增加，與式(9)一致.而且，隨著水深的增加，晃蕩波浪產生的沖擊壓力幅值逐漸減小，從h/H=0.46時P1點的 2 300 Pa到h/H=0.46時P1點的 1 200 Pa.從3.1節也可看出，隨著充液比的增加，壁面沖擊壓力的非線性變化特征逐漸減弱，從h/H=0.3時的雙壓力峰值現象過渡到h/H=0.61時的單壓力峰值現象，最后在h/H=0.9時呈現規則的正弦曲線變化特征.因為隨著充液比的增加，晃蕩液體的運動幅值會更多受到周圍液體黏性阻尼作用的約束，而且在高充液比如h/H=0.75或0.9時，從圖9(c)、9(d)中點P5的壓力幅值可看出，晃蕩波浪砰擊到液艙頂部，頂部砰擊的阻尼作用也會對波浪的非線性產生抑制作用.另外，在充液比h/H=0.46或0.61時，底部斜邊點P1和壁面下部P2點遭受巨大的沖擊壓力，超過 2 000 Pa，壁面中上部的沖擊壓力也較大，接近 1 500 Pa，如圖9(a)、9(b)，說明此水深下整個垂向壁面都遭受較大的砰擊載荷，液艙壁的整體強度需要加強.在充液比增加到h/H=0.75或0.90時，艙壁的整體沖擊壓力顯著降低，但艙壁上斜邊和頂部會遭受顯著的沖擊壓力.由此可知，載液船舶在海上航行時，理想的充液水深是接近滿載時.當充液水深近似為液艙高度一半時，液艙處于非常不利的工作環境.

4 結論

采用基于直角網格法的三維多相流模型模擬了橫搖激勵下菱形液艙大幅晃蕩問題，研究了不同水深下艙壁局部位置壓力幅值與激勵頻率之間的關系，得出菱形液艙的共振頻率，主要結論有：

(1) 高充液比下晃蕩研究表明本文計算模型具有良好的網格和時間收斂性，不同充液比下的沖擊壓力和波浪爬高也與實驗數據吻合，驗證了本文方法的精度.

(2) 隨著充液比增加，壁面沖擊壓力的非線性變化特征逐漸減弱，從充液比h/H=0.3時的雙壓力峰值現象過渡到h/H=0.61時的單壓力峰值現象，最后在h/H=0.9時呈現規則的正弦曲線變化特征.

(3) 同一充液比下菱形液艙自然頻率低于同主尺度矩形液艙的自然頻率，隨著充液比增加，最大的壓力幅值逐漸減小，因為有更多的黏性阻尼作用，而且高充液比時，也存在頂部砰擊的阻尼作用.

(4) 當充液水深近似為液艙高度一半時，整個液艙垂向艙壁都遭受較大的砰擊載荷，液艙處于非常不利的工作環境.