，

(大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024)

對于LNG運(yùn)輸船上的低溫液艙而言，外部熱量會通過絕熱維護(hù)系統(tǒng)進(jìn)入液艙內(nèi)部，會導(dǎo)致低溫液體蒸發(fā)，造成艙內(nèi)壓力升高及液貨損失。目前，國內(nèi)大多數(shù)的研究是從熱力學(xué)角度通過解析公式對低溫液艙蒸發(fā)率進(jìn)行了計算[1-3]，但沒有直觀揭示低溫液艙溫度、壓力及流場的分布規(guī)律；或者是利用公式計算了壁面熱流密度，并利用二維模型研究了儲罐內(nèi)部流場、溫度場以及壓力場的分布規(guī)律[4]。國外有學(xué)者通過4種方法從二維模型到三維模型數(shù)值計算了液艙絕熱維護(hù)系統(tǒng)的溫度場[5]，但是沒有研究液艙內(nèi)部流體相變及其流場的分布規(guī)律。為此，考慮采用三維模型，將絕熱維護(hù)系統(tǒng)及液艙內(nèi)部低溫流體蒸發(fā)耦合進(jìn)行計算，得到整個系統(tǒng)的流場、溫度場、壓力場等分布規(guī)律。

1 數(shù)值模型和數(shù)值方法

1.1 VOF方法

采用VOF方法來追蹤相間界面，求解一套守恒方程并通過求解控制單元內(nèi)單相或多相體積分?jǐn)?shù)的連續(xù)性方程來確定相間界面，對于第q相來說，方程如下。

▽·(αqρqVq)=

(1)

方程(1)中主項(xiàng)的體積分?jǐn)?shù)求解式為

(2)

1.2 蒸發(fā)模型

在液艙內(nèi)部飽和狀態(tài)下，克拉佩龍方程描述了壓強(qiáng)和溫度的相互關(guān)系。

(3)

式中：ρg、ρl分別表示氣、液相的密度。一般情況下，飽和狀態(tài)在交界面附近，而氣相部分的溫度高于飽和溫度，此時壓強(qiáng)的變化率與氣相溫度之間的關(guān)系如下[6]。

(4)

式中：R為氣體常數(shù);V為氣相體積;Pi、Tgi分別為初始壓強(qiáng)和初始溫度。

壁面漏熱量會引起液艙內(nèi)部流體能量的變化，同時引起傳熱傳質(zhì)過程，在氣液交界面處的傳熱傳質(zhì)過程可以描述為[7]

(5)

選擇Lee模型來確定蒸發(fā)過程，模型中氣相和液相的質(zhì)量轉(zhuǎn)移由下面的輸運(yùn)方程控制。

(6)

在引入飽和溫度時，質(zhì)量轉(zhuǎn)移用下式描述。

Tl>Tsat時，

(7)

Tg

(8)

1.3 控制方程

1)質(zhì)量守恒方程。

▽·(ρV)=Sm

(9)

式中：V是速度矢量；Sm為源相。

2)動量守恒方程。

▽·(ρVV)=

-▽p+▽·(τ)+ρg+F

(10)

式中：g為用在微元體上的重力；F為其他外部體積力，同時還包括其他模型的相關(guān)源相;τ為因分子粘性作用而產(chǎn)生的作用在微元體表面上的粘性應(yīng)力張量。

1.4 能量守恒方程

式中：keff為有效導(dǎo)熱系數(shù)；hj為j組分的焓；Jj為j組分的擴(kuò)散通量；τeff為粘性應(yīng)力張量；Sh為能量源相。

2 低溫物理模型實(shí)驗(yàn)

1.5 m3臥式獨(dú)立C型低溫液艙見圖1，液艙尺寸見表1，絕熱維護(hù)系統(tǒng)所用材料及其屬性見表2。

圖1中①、②、③、④處分別安裝溫度傳感器，在管路處安裝壓力傳感器。

實(shí)驗(yàn)裝置包括實(shí)驗(yàn)液艙及支撐裝置、管路系統(tǒng)及控制閥件、溫度傳感器、壓力傳感器、壓力表，數(shù)據(jù)采集系統(tǒng)、地磅、空調(diào)、溫度濕度計、電腦、梯子、千斤頂?shù)容o助工具。實(shí)驗(yàn)介質(zhì)為液氮(標(biāo)準(zhǔn)大氣壓下，沸點(diǎn)：-195.8 ℃)。實(shí)驗(yàn)工況為充裝率49.13%。

表1 液艙尺寸 mm

表2 材料屬性

實(shí)驗(yàn)流程。在實(shí)驗(yàn)裝置準(zhǔn)備完畢后，進(jìn)行儀表標(biāo)定調(diào)試，并設(shè)置空調(diào)參數(shù)保證室內(nèi)溫度，然后開始對液艙進(jìn)行預(yù)冷，同時檢查閥門、儀表及采集系統(tǒng)是否正常工作。預(yù)冷工作就緒后，進(jìn)行加注，然后關(guān)閉所有閥門憋壓靜置。當(dāng)整個系統(tǒng)溫度場達(dá)到相對穩(wěn)定，打開氣相管子出口閥門泄壓，保持閥門一直開放，進(jìn)行0表壓(即基本與大氣壓強(qiáng)相等)蒸發(fā)過程的監(jiān)測實(shí)驗(yàn)。實(shí)驗(yàn)過程中用溫度、壓力傳感器進(jìn)行溫度、壓力監(jiān)測，并實(shí)時采集重量數(shù)據(jù)。

處理實(shí)驗(yàn)數(shù)據(jù)，提取充裝率49.13%工況下，①、②、③、④處溫度傳感器，以及管路處壓力傳感器對應(yīng)的測量值見表3。

表3 充裝率49.13%實(shí)驗(yàn)結(jié)果

3 數(shù)值計算

選用Fluent軟件進(jìn)行數(shù)值模擬，計算分為兩個階段。

1)固體域傳熱計算。只考慮固體熱傳導(dǎo)，對絕熱層、鞍座、內(nèi)殼及管路系統(tǒng)進(jìn)行計算，得到相對穩(wěn)定的溫度場，然后輸出各處所對應(yīng)的溫度值。

2)整個系統(tǒng)耦合計算。整理第1階段得到的數(shù)據(jù)，并編寫初始化UDF，以初始條件的方式將第1階段的計算結(jié)果加載到計算模型固體域上；液艙內(nèi)部流體，分別對氣相、氣液界面附近、液相3個區(qū)域給定不同的初始溫度，將絕熱維護(hù)系統(tǒng)以及低溫液艙內(nèi)的流體耦合起來進(jìn)行計算。

3.1 建模及網(wǎng)格劃分

計算模型具有對稱性，因此采用四分之一模型進(jìn)行數(shù)值模擬，見圖2，計算區(qū)域包括內(nèi)部流體域、內(nèi)殼及管路系統(tǒng)、絕熱層、鞍座。計算模型坐標(biāo)原點(diǎn)位于橫向半圓切面的圓心處，x軸沿模型縱向方向，y軸沿模型豎直方向，z軸沿模型橫向方向。采用ICEM軟件劃分網(wǎng)格，所有區(qū)域均采用結(jié)構(gòu)網(wǎng)格，整體網(wǎng)格見圖3；對內(nèi)殼及管路、鞍座進(jìn)行網(wǎng)格的細(xì)化，以滿足計算要求，見圖4。

3.2 數(shù)值模擬設(shè)置

第1階段。由文獻(xiàn)可知，獨(dú)立C型低溫液艙內(nèi)部氣相的溫度高于氣液界面處的溫度，氣液界面處的溫度高于液相的溫度。因此，在給定初始條件時，分別在內(nèi)殼及管路系統(tǒng)與氣相、氣液界面處、液相接觸處設(shè)置不同的溫度，具體見表4。絕熱層與空氣接觸表面設(shè)置為293.15 K的固定溫度壁面條件。

表4 初始溫度設(shè)置(充裝率49.13%) K

第2階段。將第1階段所得到的數(shù)據(jù)以初始化UDF的形式加載到整體模型算例中，液艙內(nèi)氣相、氣液界面處、液相初始溫度設(shè)置與表4一致。初始壓強(qiáng)給定為0.2 MPa，艙內(nèi)部流體采用液氮及氮?dú)庾鳛槟M工質(zhì)，液氮在標(biāo)準(zhǔn)大氣壓下的物理屬性見表5。

表5 標(biāo)準(zhǔn)大氣壓下液氮物理屬性

采用壓力基瞬態(tài)計算，采用VOF多相流模型，并包含Implicit Body Force，采用標(biāo)準(zhǔn)k-ε雙方程湍流模型，以及標(biāo)準(zhǔn)壁面函數(shù)，同時計算Viscous Heating。選擇Lee模型作為蒸發(fā)冷凝模型，并設(shè)置飽和溫度為77.35 K，表面張力選擇Continuum Surface Stress模型，并設(shè)置表面張力為0.008 59 N/m。絕熱層與空氣接觸壁面設(shè)置為293.15 K的固定溫度壁面條件，其他區(qū)域接觸面設(shè)置為coupled熱邊界條件，內(nèi)殼、管路系統(tǒng)與流體接觸面設(shè)置為無滑移速度條件。管子出口設(shè)置為pressure outlet。壓力速度耦合采用coupled算法。

4 結(jié)果分析

圖1中①～④處傳感器的位置(相對于數(shù)值模型坐標(biāo)系)見表6。

表6 傳感器位置坐標(biāo) m

提取位置①～④處數(shù)值計算結(jié)果以及管路內(nèi)的壓強(qiáng)與實(shí)驗(yàn)值比較，見表7。

表7 實(shí)驗(yàn)與數(shù)值計算結(jié)果對比

由表7可以看出，數(shù)值計算結(jié)果與實(shí)驗(yàn)結(jié)果相對誤差均小于3%，進(jìn)而也說明了數(shù)值計算方法的可靠性。

圖5a)為設(shè)置完初始溫度并加載完UDF后的溫度場，可以看出初始化UDF已將溫度參數(shù)加載到整體模型上。圖5b)可以看出儲罐內(nèi)部流體存在著溫度梯度，整體表現(xiàn)是自下而上，溫度逐漸升高，對比圖5a)可以說明，氣相和液相通過氣液界面存在傳熱現(xiàn)象，且氣相內(nèi)部也存在著熱量交換，這也造成了初始 “過熱”氣體溫度逐漸降低，且越靠近氣液界面處氣相溫度降低越明顯。

圖6顯示了管路處溫度分布，由圖6a)可以看出管路的存在使得周圍的溫度明顯低于同一高度其他各處的溫度。這是由于蒸發(fā)氣體不斷地通過管路排出，而氣體相對于絕熱層及管路來說溫度更低，且管路的導(dǎo)熱系數(shù)相對較大，管路和氣相接觸，導(dǎo)致管路的溫度迅速降低，進(jìn)而影響了管路附近區(qū)域溫度分布。由圖6b)可以看出，在管路左側(cè)0.02 m處(圖6a)中黑線)，由于低溫儲罐及管路的影響，沿豎直方向在0.51～0.83 m范圍內(nèi)，溫度由104 K變化到111 K,溫度變化較小，而從0.83～0.91 m范圍內(nèi)，溫度由111 K變化到225 K,溫度變化明顯。

另外由于鞍座材料的導(dǎo)熱系數(shù)相對較大，所以鞍座溫度也會很快降低，進(jìn)而影響其周圍溫度分布，見圖7。

由圖7可見，鞍座的存在使得其附近區(qū)域溫度低于較遠(yuǎn)區(qū)域的溫度，且鞍座處溫度最低，大約為210 K，在距鞍座0.25 m范圍內(nèi)最大溫差為30 K。

圖8顯示了不同切面處溫度分布，從圖8a)可以直觀地看出液艙內(nèi)部溫度整體趨勢為上高下低，但從圖8c)、圖8d)可以看出在底部近壁面處，溫度稍高于其上部溫度。從圖8b)、圖8c)、圖8d)可以看出液體在壁面處的溫度都高于同一水平內(nèi)部流體溫度，這是因?yàn)橥獠總鳠崾贡诿嫣幜黧w溫度升高。

由圖9可以看出，在近壁面處，液相溫度由77.35 K逐步增加到近壁面處的81.5 K。由圖8d)可以看出，在出口處氣體溫度有所升高，這是由于管路在此處直接與外界接觸，且氣體也直接與外界接觸，傳熱過程加劇所致。

從圖10可以具體地看出，在液相整體溫度相對均勻，溫度梯度非常小，而對于氣液界面附近及氣相區(qū)域溫度梯度相對較大。

圖11顯示出了內(nèi)殼及管路系統(tǒng)溫度場。

由圖11a)可以看出內(nèi)殼與管路接觸處溫度明顯低于其周邊溫度，這是由于低溫氣體在此處排出，會帶走部分熱量，致使其溫度降低；圖11b)顯示了內(nèi)殼與氣液界面接觸處，由液相溫度過渡到氣相溫度的溫度梯度；由圖11c)可以看出，內(nèi)殼與鞍座接觸處溫度高于周邊溫度，但溫度差異不明顯。

圖12顯示了氣體的體積分?jǐn)?shù)。

可以看出自氣液界面向下，氣體的體積分?jǐn)?shù)在逐漸減小，且在與壁面接觸處氣體的體積分?jǐn)?shù)并無顯著變化，這說明蒸發(fā)現(xiàn)象主要發(fā)生在氣液界面附近。這是由于壁面附近的邊界層流對流循環(huán)會將熱量攜帶到氣液界面處，再加上氣液界面處的壁面漏熱，以及氣相與液相之間的傳熱，使得氣液界面處液體溫度升高相對較大。因此，蒸發(fā)現(xiàn)象明顯，而在壁面附近由于液體及固體壁面溫差相對來說并不大，漏熱量并不顯著，同時壁面邊界層的對流循環(huán)也會將壁面漏入的熱量帶走，因而蒸發(fā)并不顯著。此外，分析可以看出，在正常的低溫液體存儲下，流體與固體壁面之間的溫差相對較小，此狀態(tài)時的漏熱量并不會引起低溫液體的沸騰，而傳質(zhì)過程以表面蒸發(fā)的方式發(fā)生。

壓強(qiáng)分布見圖13。

由圖13可見，當(dāng)處于自由蒸發(fā)狀態(tài)時，液艙內(nèi)壓強(qiáng)上低下高，且管路內(nèi)部壓強(qiáng)基本處于大氣壓強(qiáng)，而氣相部分壓強(qiáng)稍高于大氣壓強(qiáng)。由于液體壓強(qiáng)的存在導(dǎo)致液艙底部壓強(qiáng)有明顯的增大，相對氣相來說壓強(qiáng)增大4013.3 Pa，這在結(jié)構(gòu)設(shè)計時應(yīng)予以注意。

整體速度矢量見圖14(管路內(nèi)氣體速度明顯高于其他部分，為顯示出速度分布差異，對顯示范圍作了相應(yīng)的調(diào)整得到圖14)，可以看出氣相流體的速度高于液相流體速度，且壁面附近流體速度高于內(nèi)部流體速度，氣體的比熱容比液體小，所以傳入相同的熱量，氣相溫度升高的更快，溫差增大，因而氣相流體流動更劇烈。

由圖15和圖16速度矢量圖可以看出壁面處的流體沿壁面向上運(yùn)動，這是由于壁面處的流體隨著熱量的不斷流入，溫度升高，形成熱邊界層，且在浮升力作用下沿著壁面向上流動。

流線見圖17，局部放大見圖18。

由圖17、18可見，隨著外界熱量的進(jìn)入，在壁面附近形成的邊界層流沿壁面向上流動，在氣液界面處向中心流動，在此過程中一部分由于蒸發(fā)進(jìn)入到氣相，一部分滯留在氣液界面附近區(qū)域，在中心區(qū)域形成中心射流向下運(yùn)動，繼續(xù)參與液體內(nèi)部的對流循環(huán)運(yùn)動。液體的這一對流運(yùn)動過程，將熱量不斷地帶入自由液面區(qū)域，使得自由液面處的液體流速增大，對流強(qiáng)度增大。同理，在氣相也會發(fā)生類似的對流循環(huán)過程，沿壁面運(yùn)動的邊界層流一部分從管路排出，一部分向下運(yùn)動，形成回流區(qū)域。圖19所示的整體流線圖，也表明了在氣液界面處及氣相靠近壁面處，對流強(qiáng)度較大，換熱明顯。

5 結(jié)論

1)對于固體域溫度分布而言，鞍座、管路與液艙內(nèi)殼接觸處，溫度分布會有變化，這在設(shè)計建造時應(yīng)予以重視；對于液艙內(nèi)部的流體而言，整體溫度呈下低上高分布，同時壁面附近處溫度高于內(nèi)部區(qū)域溫度，且相對液相來說，氣液界面處以及氣相區(qū)域溫度梯度較大。

2)相對壁面附近而言，蒸發(fā)現(xiàn)象在氣液界面處更加明顯。

3)自由蒸發(fā)時，液艙內(nèi)部壓強(qiáng)呈現(xiàn)上低下高分布，氣相區(qū)域壓強(qiáng)基本均勻，且稍高于大氣壓強(qiáng)，而液相由于液體壓強(qiáng)的存在，隨著深度的增加，壓強(qiáng)逐漸增大。

4)氣體的速度高于液體，且在氣相近壁面處及管路內(nèi)尤為明顯，液相沿壁面向上流動的邊界層流一部分蒸發(fā)進(jìn)入氣相，一部分沿界面向中心流動且在界面下方形成回流區(qū)域。

5)液艙內(nèi)部氣相近壁面處及氣液界面附近自然對流現(xiàn)象相對于其他區(qū)域更為明顯。

相對于目前的研究而言，本文所述將低溫液艙絕熱維護(hù)系統(tǒng)及艙內(nèi)低溫流體相變耦合進(jìn)行計算的方法更符合實(shí)際情況，更便于分析整個系統(tǒng)的相應(yīng)規(guī)律，但隨著低溫液艙的大型化，前期初始條件加載需要處理的溫度數(shù)據(jù)會顯著增多，在進(jìn)一步應(yīng)用時可以考慮將數(shù)據(jù)分組回歸出函數(shù)表達(dá)式的形式進(jìn)行加載。本研究所得到的蒸發(fā)特性規(guī)律，對低溫液艙的設(shè)計建造具有重要的參考價值，同時在得到整個系統(tǒng)的溫度場分布時，可以進(jìn)一步提取溫度數(shù)據(jù)，對結(jié)構(gòu)的熱應(yīng)力進(jìn)行深入分析，從而更好地優(yōu)化設(shè)計。

[1] 劉文華,陸晟.中小型LNG船C型獨(dú)立液貨艙蒸發(fā)率計算[J].船舶設(shè)計通訊,2012(1):25-28.

[2] 時光志,盛蘇建.獨(dú)立C型液貨艙的傳熱分析及蒸發(fā)率計算[J].船海工程,2013,42(1):65-69.

[3] 田竹剛,尹奇志,王冰.船用LNG燃料儲罐蒸發(fā)量與維持時間計算方法[J].交通信息與安全,2016,34(3):57-63.

[4] 姚壽廣,閆興武,葉勇,等.某型船用LNG儲罐熱響應(yīng)仿真研究[J].江蘇科技大學(xué)學(xué)報(自然科學(xué)版),2016,30(2):136-141.

[5] MIANA M, LEGORBURO R, DEZ D, et al. Calculation of boil-off rate of liquefied natural gas in mark III tanks of ship carriers by numerical analysis[J]. Applied Thermal Engineering,2016,93:279-296.

[6] HOPFINGER E J, DAS S P. Mass transfer enhancement by capillary waves at a liquid vapour interface[J]. Experiments in Fluids,2009,46(4):597-605.

[7] GROTLE E J, et al. Non-isothermal sloshing in marine liquefied natural gas fuel tanks[C]. The Twenty-sixth International Ocean and Polar Engineering Conference,2016:879-885.