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

劉東喜， 莊宿國(guó)， 王 晉， 尤云祥

(1. 上海海事大學(xué) 海洋科學(xué)與工程學(xué)院， 上海 201306； 2. 西安航天動(dòng)力研究所， 西安 710100； 3. 上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200240)

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

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

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

1 晃蕩問(wèn)題描述

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

根據(jù)線(xiàn)性勢(shì)流理論計(jì)算可得圖1中自由表面、C-W界面和W-D界面的共振晃蕩頻率.對(duì)于某二維矩形艙內(nèi)多層液體系統(tǒng)，矩形艙長(zhǎng)度為L(zhǎng)，艙內(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.令上述線(xiàn)性系統(tǒng)的無(wú)阻尼矩陣行列式等于0，則行列式可表示為固有頻率ω的6階多項(xiàng)式，即[3]

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

(1)

(2)

(3)

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

2 數(shù)值模型

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

(4)

(5)

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

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

(6)

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

H(φ)=

(7)

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

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

(8)

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

(9)

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

ni=?φ/?xi

(10)

(11)

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

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

(12)

函數(shù)δ(φ)可通過(guò)對(duì)H(φ)求導(dǎo)得到，即

(13)

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

3 自由衰減測(cè)試

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

圖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)和線(xiàn)性勢(shì)流理論計(jì)算值(ωcal).可以看出，其模擬值與理論值吻合良好，初步驗(yàn)證了本文所采用的數(shù)值模擬方法的準(zhǔn)確性.

4 結(jié)果與分析

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

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

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

圖4 左側(cè)艙壁處界面高度時(shí)程的計(jì)算值與試驗(yàn)結(jié)果對(duì)比Fig.4 Comparison of elevation time history at the left wall by simulation and test results

圖5 自由表面和液-液界面形狀的計(jì)算與試驗(yàn)結(jié)果對(duì)比Fig.5 Comparison of two interfaces obtained by simulation and test

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

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

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

5 結(jié)語(yǔ)

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