矩形艙內(nèi)三層液體晃蕩特性的數(shù)值分析

劉東喜， 莊宿國， 王 晉， 尤云祥

(1. 上海海事大學 海洋科學與工程學院， 上海 201306； 2. 西安航天動力研究所， 西安 710100； 3. 上海交通大學 海洋工程國家重點實驗室， 上海 200240)

隨著能源需求的不斷增長和近海油氣資源的日益減少，海洋油氣開采逐漸向深海區(qū)域發(fā)展，越來越多的浮式平臺投入使用.為進一步提高浮式生產(chǎn)儲卸油(Floating Production Storage and Offloading，F(xiàn)PSO)裝置的生產(chǎn)效率，國外石油公司提出在船體內(nèi)安裝大型清洗艙以作為傳統(tǒng)油-氣-水分離器的預處理設備.與分離器一樣，清洗艙內(nèi)油層與水層之間也存在一定厚度的乳化層.在風、浪、流等海洋環(huán)境載荷的共同作用下，F(xiàn)PSO裝置將會產(chǎn)生6個自由度的運動，從而導致平臺清洗艙內(nèi)多層液體的晃蕩.

目前，國內(nèi)外對多層液體晃蕩現(xiàn)象的研究主要采用數(shù)值模擬和模型試驗的方法.王日新等[1]采用標記網(wǎng)格(Marker and Cell，MAC)方法對縱搖激勵下矩形艙內(nèi)兩層液體的晃蕩情況進行了數(shù)值模擬；Xue等[2]采用流體體積(Volume of Fraction，VOF)方法和試驗方法研究了矩形艙內(nèi)油層與水層液體的晃蕩問題；Molin等[3]采用線性勢流理論方法和試驗方法研究了FPSO裝置矩形清洗艙內(nèi)三層液體的晃蕩情況；Kim等[4]等采用一種改進的移動粒子半隱式(Moving Particle Semi-implicit，MPS)方法模擬了矩形艙內(nèi)三層液體的晃蕩問題.但是，與傳統(tǒng)的單層液體晃蕩問題相比，多層液體晃蕩問題的相關研究成果還較少，有待進一步深入研究[5].

本文采用Sussman等[6]提出的耦合水平集和流體體積 (Coupled Level Set and Volume of Fluid, CLSVOF)界面捕獲方法對FPSO裝置矩形清洗艙內(nèi)三層液體的晃蕩特性進行數(shù)值模擬.開展了自由衰減測試試驗，得到了自由表面和兩個液-液界面的固有頻率，研究了液-液界面產(chǎn)生的Kelvin-Helmholtz不穩(wěn)定性，同時，通過對比模擬與試驗所得界面形狀和界面高度的變化情況，從而驗證了所用數(shù)值模擬方法的準確性.

1 晃蕩問題描述

本文采用的矩形艙內(nèi)三層液體晃蕩算例由Molin等[3]提供.矩形艙系統(tǒng)模型如圖1所示，其中，坐標x、y和z軸方向分別表示矩形艙長度方向、高度方向和寬度方向.矩形艙的長度、高度和垂直紙面寬度分別為 1.08、0.90 和 0.10 m.艙內(nèi)充裝三層互不相溶的液體，其物理性能參數(shù)見表1.表中：ρ為密度；ν為運動黏度；σ為表面張力系數(shù).

根據(jù)線性勢流理論計算可得圖1中自由表面、C-W界面和W-D界面的共振晃蕩頻率.對于某二維矩形艙內(nèi)多層液體系統(tǒng)，矩形艙長度為L，艙內(nèi)充裝三層液體，其厚度和密度分別為he和ρe(e=1，2，

圖1 矩形艙內(nèi)三層液體的晃蕩模型(m)Fig.1 Schematic model for three-layer liquid sloshing in a rectangular tank (m)

液體ρ/(g·cm-3)ν/(mm2·s-1)σ/(mN·m-1)二氯甲烷1.300.327.8水1.001.072.7環(huán)己烷0.781.324.7

3)，且ρ1>ρ2>ρ3.令上述線性系統(tǒng)的無阻尼矩陣行列式等于0，則行列式可表示為固有頻率ω的6階多項式，即[3]

Ψ6ω6+Ψ4ω4+Ψ2ω2+Ψ0=0

(1)

(2)

(3)

式中：g為重力加速度；n為模數(shù).

2 數(shù)值模型

矩形艙內(nèi)流體運動的控制方程為連續(xù)方程和動量方程，即

(4)

(5)

式中：下標i、j為x軸的方向；ui、uj分別為xi、xj方向的速度；p為壓力；fi為外部加速度；Fi為表面張力.

采用Sussman等[6]提出的CLSVOF方法分析三層液體晃蕩過程中的自由表面和液-液界面的運動特性.其中：VOF方法是通過在每個計算單元上定義一個體積分數(shù)α來實現(xiàn)界面追蹤的，且α=0表示該單元被上層流體充滿，α=1 表示該單元被下層流體充滿，0<α<1 表示界面跨過該單元；水平集方法是通過定義等值函數(shù)φ來實現(xiàn)對液-液界面追蹤的，且φ=0處為界面，某一點的φ值為該點到界面的距離，φ<0表示界面上層流體，φ>0表示界面下層流體.函數(shù)φ和α的輸運方程為

(6)

為避免流體物理性能在界面出現(xiàn)間斷，引入一個光滑的Heaviside函數(shù)H(φ)以對界面很小范圍內(nèi)的流體物理性能進行連續(xù)化處理，H(φ)的形式為

H(φ)=

(7)

式中：h為網(wǎng)格尺寸.密度和黏度的控制方程為

ρ(φ)=ρaH(φ)+ρb[1-H(φ)]

(8)

μ(φ)=μaH(φ)+μb[1-H(φ)]

(9)

式中：下標b和a分別為上層和下層流體. 界面法向向量n和曲率κ(φ)的計算公式如下：

ni=?φ/?xi

(10)

(11)

式中：xk表示k方向分量.

Funada等[7]發(fā)現(xiàn)，表面張力對界面不穩(wěn)定性以及界面短波的影響顯著.本文采用CSF(Continuum Surface Force)模型[8]計算表面張力，即

(12)

函數(shù)δ(φ)可通過對H(φ)求導得到，即

(13)

為了精確求解，必須始終保持φ為距離函數(shù)，即|?φ/?xk|=1，但是經(jīng)過幾步的計算后，φ將不再滿足距離函數(shù)的條件.在保持界面位置不變的情況下，需要重新構造等值函數(shù)φ.

3 自由衰減測試

為得到自由表面和兩個液-液界面固有頻率的數(shù)值模擬值，本文進行了自由衰減測試，并通過模擬值與理論值的對比來初步驗證數(shù)值模擬方法的準確性.根據(jù)圖1中矩形艙模型的尺寸和三層液體的充裝高度，采用網(wǎng)格劃分軟件ICEM生成一組單元數(shù)為 77.76 萬的三維精細網(wǎng)格系統(tǒng).計算時間步長設定為 0.5 ms.計算開始之前，將矩形艙內(nèi)自由表面和兩個液-液界面初始化為傾斜狀態(tài)，傾斜幅度為5°；計算開始之后，三層液體將在重力作用下做自由振蕩；計算過程中，實時記錄艙壁處3個界面高度的變化數(shù)據(jù).由于根據(jù)線性勢流理論方法(式(1)～(3))得到的固有頻率沒有考慮流體黏性，所以本文的自由衰減測試也假設流體為無黏性流體.圖2示出了距離矩形艙左側艙壁5 mm處3個界面的高度變化曲線對應的頻譜密度P.圖中，譜密度峰值對應的頻率即為固有頻率.

圖2 自由表面和液-液界面高度的頻譜密度Fig.2 Power spectrums of free decay of sloshing waves

Tab.2 Simulated and theoretical values of natural frequencies of the first four sloshing modes

模態(tài)階數(shù)ωsim/(rad·s-1)自由面C-W面W-D面ωcal/(rad·s-1)自由面C-W面W-D面15.1561.9090.9555.2041.8381.05027.6392.8651.9097.5502.9491.98639.3573.6282.8649.2523.6082.775410.6944.2013.62810.6804.0803.436

表2列出了前4階模態(tài)晃蕩固有頻率的自由衰減數(shù)值模擬值(ωsim)和線性勢流理論計算值(ωcal).可以看出，其模擬值與理論值吻合良好，初步驗證了本文所采用的數(shù)值模擬方法的準確性.

4 結果與分析

艙體繞z軸做受迫橫搖運動，其運動方程為θ(t)=θesin(ωet)，θ為矩形艙的橫搖角度，θe=1° 為橫搖角的幅值，ωe=1.83 rad/s為激勵頻率.由表2可知，該激勵頻率等于中部C-W界面的最低階固有頻率.

圖3示出了矩形艙左側艙壁處頂部自由表面、中部C-W界面和底部W-D界面高度隨時間的變化曲線.可以看出：當激勵頻率等于中部C-W界面的最低階固有頻率時，C-W界面和W-D界面在振蕩初期出現(xiàn)了共振，約20個振蕩周期后，其運動趨于穩(wěn)定；C-W界面的運動幅度明顯大于另外兩個界面，該結果符合預期；W-D界面與C-W界面的振蕩頻率相同，但其振蕩幅值小于C-W界面；自由表面的運動幅度較弱，振蕩幅值遠小于另外兩個內(nèi)部界面，故其振蕩相位也稍微偏離另外兩個內(nèi)部界面；另外，W-D界面高度的變化曲線呈現(xiàn)出峰寬谷窄的特征，且峰值明顯小于谷值，該現(xiàn)象與非線性自由表面波的特征相反.

圖3 左側艙壁處自由表面和液-液界面的高度變化曲線Fig.3 Elevation time history of three interfaces at left wall of the tank

圖4 左側艙壁處界面高度時程的計算值與試驗結果對比Fig.4 Comparison of elevation time history at the left wall by simulation and test results

圖5 自由表面和液-液界面形狀的計算與試驗結果對比Fig.5 Comparison of two interfaces obtained by simulation and test

圖4所示為左側艙壁處3個界面高度時程的計算結果與Molin等[3]的試驗結果對比.可以看出，計算值與試驗值吻合良好，從而進一步驗證了本文所用數(shù)值模擬方法的準確性.需要說明的是，試驗所得 W-D 界面高度的變化曲線中出現(xiàn)的上、下抖動可能是由于傳感器失效等測量誤差造成的.圖5示出了兩個不同時刻的自由表面和液-液界面形狀的計算結果和試驗結果.可見兩者吻合良好.內(nèi)部界面晃蕩運動的幅度遠大于頂部自由表面.另外，由于上述激勵頻率等于中間層液體的最低階固有頻率，且接近于底層液體的2階固有頻率，所以由圖5可以觀察到中部C-W界面的1階模態(tài)和底部W-D界面的2階模態(tài).

Funada等[7]指出，當兩種流體的交界面出現(xiàn)切向速度間斷時，流體界面是不穩(wěn)定的，將會產(chǎn)生Kelvin-Helmholtz不穩(wěn)定性.在這種不穩(wěn)定性的作用下，流體界面產(chǎn)生的小擾動經(jīng)過線性和非線性的增長后會在強的非線性作用下發(fā)展成湍流混合.當激勵頻率接近于矩形艙多層液體系統(tǒng)某一液-液界面的最低階固有頻率時，界面的大幅度共振運動將會導致界面出現(xiàn)Kelvin-Helmholtz不穩(wěn)定性.圖6示出了t=96.3 s時由Kelvin-Helmholtz不穩(wěn)定性所致C-W界面出現(xiàn)的鋸齒狀短波，該現(xiàn)象進一步表明多層液體晃蕩問題比單層液體晃蕩問題復雜得多.同時，圖6對比了C-W界面Kelvin-Helmholtz(K-H)不穩(wěn)定性的數(shù)值模擬結果與Molin等[3]的試驗結果.由圖6中界面鋸齒波的形狀以及相應的短波長度來看，本文的模擬結果與試驗結果非常接近，表明本文所用數(shù)值模擬方法不僅能夠準確預測液-液界面的整體晃蕩行為(界面晃蕩運動)，而且能夠準確預測液-液界面的小尺度物理現(xiàn)象(K-H不穩(wěn)定性).

圖6 C-W界面K-H不穩(wěn)定性的數(shù)值模擬與試驗結果對比Fig.6 Comparison between numerical solution and test results of K-H instability at the C-W interface

5 結語

本文采用CLSVOF界面捕獲方法對FPSO裝置矩形清洗艙內(nèi)三層液體的晃蕩特性進行了數(shù)值分析，并將模擬結果與試驗結果進行對比.結果表明，多層液體儲卸艙系統(tǒng)存在不止一種晃蕩固有頻率.當清洗艙外部激勵頻率接近于液-液界面的最低階固有頻率時，將產(chǎn)生共振效應并導致液-液界面的晃蕩幅度遠大于自由表面，表明即使清洗艙的受迫運動幅度不大，但只要發(fā)生內(nèi)部共振，大幅度的內(nèi)部晃蕩運動也會明顯降低其分離性能.因此，在設計多層液體儲卸艙系統(tǒng)時，不僅要考慮頂部自由表面的運動特性，而且要考慮內(nèi)部界面的運動特性.另外，共振效應也會導致液-液界面出現(xiàn)復雜的鋸齒狀短波.所用CLSVOF方法不僅能夠預測自由表面和液-液界面的大尺度整體晃蕩行為，而且能夠表征液-液界面出現(xiàn)的小尺度物理現(xiàn)象，適用于多層液體晃蕩特性的數(shù)值分析.