王敏，邵倩倩，楊曉帆，李敬法

（1 北京計算科學(xué)研究中心，北京100193； 2 廣東石油化工學(xué)院石油工程學(xué)院，廣東茂名525000；3 北京師范大學(xué)地理科學(xué)學(xué)部，北京100875； 4 北京石油化工學(xué)院機械工程學(xué)院，北京102617）

引 言

浮頂油罐是最重要的儲油設(shè)施。含蠟原油在儲存過程中會因為向外界傳熱而冷卻、析蠟和膠凝[1]。為保證罐內(nèi)含蠟原油能夠隨時外輸，通常需要定期、快速地對浮頂油罐進(jìn)行加熱。而盤管組是最常用的油罐加熱融化方式之一。盤管組的布置方式將對罐內(nèi)含蠟原油融化過程產(chǎn)生較大的影響[2-6]。因此，開展不同盤管組布置方式下，浮頂油罐內(nèi)含蠟原油融化規(guī)律研究，對于科學(xué)設(shè)計盤管組的布置方式、合理制定浮頂油罐加熱融化方案、削減成本等具有重要的意義。

盤管組加熱的浮頂油罐內(nèi)含蠟原油加熱融化問題是當(dāng)前油罐研究的熱點之一，而且這方面研究目前已經(jīng)取得了一定的成果[2-8]。針對油輪燃料艙重油加熱問題，Magazinovic[2]提出了一種緊湊型盤管組布置方式。為了考察新方案的加熱效果，建立了盤管組加熱的燃料艙傳熱物理數(shù)學(xué)模型，采用數(shù)值計算方法對艙內(nèi)重油加熱溫度場、流場以及盤管組對流傳熱系數(shù)進(jìn)行計算，評估新型盤管組布置方式的加熱效果。針對1.0×105m3浮動油罐，Sun 等[3]分別采用實驗和數(shù)值模擬方法研究了盤管加熱條件下，罐內(nèi)含蠟原油的溫升規(guī)律。此外，從盤管傳入油罐的有效能利用角度出發(fā)，對盤管的管徑和長度進(jìn)行優(yōu)選。在考慮外界氣溫、太陽輻射、盤管內(nèi)蒸氣壓力以及原油物性的基礎(chǔ)上，Sun 等[4]建立了盤管組加熱的原油儲罐傳熱物理數(shù)學(xué)模型，研究了盤管組結(jié)構(gòu)、液高以及盤管組熱流量對罐內(nèi)原油溫升過程的影響。Zhao 等[5]采用數(shù)值計算方法對盤管組加熱的1.0×104m3浮頂油罐內(nèi)原油區(qū)域流動與傳熱過程進(jìn)行計算，對罐內(nèi)油溫隨時間的變化進(jìn)行分析，并對罐內(nèi)縱向和橫向油溫隨時間的演化規(guī)律進(jìn)行研究。孫巍[6]建立了大型浮頂油罐內(nèi)盤管加熱的原油溫升物理數(shù)學(xué)模型，研究了不同盤管排布方式、油品液位、加熱蒸汽溫度、外界環(huán)境等條件下，罐內(nèi)原油溫度場與流場的變化規(guī)律。陸雅紅等[7]利用Fluent 軟件對原油儲罐蒸氣盤管加熱溫度場進(jìn)行模擬，以熱流量為主要指標(biāo)，對蒸氣盤管的結(jié)構(gòu)進(jìn)行優(yōu)化設(shè)計。劉鳳榮[8]以大慶地區(qū)某2 萬方浮頂油罐為研究對象，建立了原油儲罐管式加熱物理數(shù)學(xué)模型，研究了冬季浮頂油罐加熱作用下罐內(nèi)原油溫度場和流場分布規(guī)律。Saifed-Din 等[9]分別采用實驗和數(shù)值模擬方法研究了臥式儲罐內(nèi)水平圓管加熱的罐內(nèi)熱水流場和溫度場變化規(guī)律，分析了圓管數(shù)量對圓管傳熱系數(shù)、流場和溫度場以及加熱效率的影響?？紤]到盤管組布置條件下，浮頂油罐內(nèi)加熱融化區(qū)域的復(fù)雜性，一些研究人員將罐底盤管組簡化為一個恒熱流的加熱平面，以此來研究罐內(nèi)介質(zhì)的溫升規(guī)律[10-12]。針對盤管加熱的10×104m3浮頂油罐溫升過程，劉佳等[10]將盤管組簡化為水平加熱面，通過調(diào)整水平面的熱通量來模擬不同的盤管組熱流量，基于該模型研究水平面加熱下罐內(nèi)原油的加熱溫度場、流場的演化規(guī)律，在此基礎(chǔ)上，對不同液位、不同蒸汽流量對罐內(nèi)原油溫度場、流場的影響進(jìn)行分析。以上這些研究主要是針對盤管組加熱時罐內(nèi)介質(zhì)的溫升過程展開，并未考慮介質(zhì)可能發(fā)生的融化相變問題。

在盤管組加熱的儲罐介質(zhì)融化研究方面，Wang等[13]采用層流模型對小尺寸模擬罐內(nèi)單根水平圓管恒溫加熱的含蠟原油融化過程進(jìn)行計算，研究了罐內(nèi)含蠟原油的融化規(guī)律，分析了含蠟物性對水平圓管熱流量及罐內(nèi)含蠟原油融化液相分?jǐn)?shù)的影響。相變儲能也是儲罐內(nèi)材料融化研究的重要方面[14-17]。Ahmadi 等[14]建立了螺旋狀加熱盤管罐內(nèi)石蠟融化儲能物理數(shù)學(xué)模型，研究了盤管直徑、螺旋管直徑、盤管流量、進(jìn)口溫度等對罐內(nèi)石蠟融化過程的影響。然而，以上研究均針對小尺寸油罐的層流融化展開，研究結(jié)果與實際浮頂油罐內(nèi)含蠟原油湍流融化的結(jié)果存在較大差異。

綜上所述，目前在盤管組加熱的浮頂油罐內(nèi)含蠟原油融化研究方面存在以下幾方面的不足：（1）在實際儲罐加熱融化研究方面，現(xiàn)有的研究主要關(guān)注油罐溫升過程，較少涉及罐內(nèi)含蠟原油的融化相變過程；（2）僅有的一些油罐融化研究基本都是針對小尺寸油罐內(nèi)層流融化開展的，該過程與實際浮頂油罐內(nèi)含蠟原油湍流融化過程差異較大；（3）在盤管組結(jié)構(gòu)對油罐加熱效果的影響研究方面，現(xiàn)有的研究主要關(guān)注盤管組結(jié)構(gòu)對溫升速率、罐內(nèi)原油溫度場均勻程度的影響，并未從罐內(nèi)流場、溫度場以及盤管組熱流量演化角度來分析盤管組布置方式對油罐加熱和融化過程的影響。基于此，針對1000 m3雙盤式浮頂油罐，在考慮含蠟原油融化過程中形態(tài)和物性變化的基礎(chǔ)上，本文建立了盤管組加熱的浮頂油罐內(nèi)含蠟原油湍流融化物理數(shù)學(xué)模型。考慮到罐內(nèi)含蠟原油區(qū)域的不規(guī)則性，引入浸入邊界法，在結(jié)構(gòu)化網(wǎng)格上處理浮頂油罐內(nèi)各區(qū)域耦合傳熱問題，實現(xiàn)浮頂油罐的一體化耦合求解。利用該模型，對不同盤管組布置方式下罐內(nèi)含蠟原油的融化規(guī)律進(jìn)行研究，并對盤管組熱流量及其分配情況進(jìn)行分析。

1 物理數(shù)學(xué)模型

1.1 物理模型

浮頂油罐內(nèi)傳熱過程包括罐內(nèi)氣-液-固耦合傳熱（空氣層、含蠟原油與盤管組、鋼板層、保溫層之間的傳熱）、浮頂油罐與外界大氣的強制對流換熱以及浮頂油罐與罐底土壤層的導(dǎo)熱，傳熱過程復(fù)雜。融化過程中，含蠟原油的形態(tài)和物性還將發(fā)生顯著變化。具體而言，初始時刻罐內(nèi)油溫較低，蠟晶多孔介質(zhì)形成。此時，含蠟原油表現(xiàn)出非牛頓流體特性。加熱時，盤管附近蠟晶多孔介質(zhì)溫度快速升高，蠟晶多孔介質(zhì)逐漸破壞，固液分散體系形成。此時含蠟原油仍表現(xiàn)出非牛頓流體特性。繼續(xù)加熱時，蠟晶開始融化，純液態(tài)含蠟原油逐漸形成。具體過程見文獻(xiàn)[13]。

圖1 盤管組加熱的浮頂油罐內(nèi)含蠟原油融化物理模型Fig.1 Physical model of melting of waxy crude oil in a floating roof tank heated by coils

此外，浮頂油罐內(nèi)含蠟原油流動處于湍流狀態(tài)（對于1000 m3浮頂油罐，當(dāng)特征溫差為50℃時，Rayleigh 數(shù)高達(dá)3.0×1013），計算時間過長?？紤]到油罐的軸對稱特性，本文將三維油罐簡化為二維油罐，如圖1 所示。此外，罐底部布置有半徑為r0盤管，壁面溫度恒為T0，間距恒為δr，盤管組中心連線與r軸正方向夾角為θ，如圖1(a)所示。

1.2 數(shù)學(xué)模型

盤管組加熱的浮頂油罐內(nèi)含蠟原油湍流融化和溫升控制方程如式（1）～式（3）所示。其中，罐內(nèi)湍流由大渦模擬方法計算[18-22]。

源項Sui和ST計算式為

其中，Kd為蠟晶多孔介質(zhì)滲透率，由Kozeny-Carman 公式計算[13,22]。fi和q分別為由浸入邊界法引起的分布力密度和分布熱密度、ΔH為相變潛熱，本文中ΔH=(1-gs)L，其中，L為含蠟原油相變總潛熱。

1.2.1 融化相變過程 式（6）等號右邊第一項代表了含蠟原油融化相變過程中蠟晶與含蠟原油之間的相互作用。-ui和-ui,s分別為含蠟原油和蠟晶的速度向量。本文計算中假設(shè)，當(dāng)油溫高于或等于含蠟原油顯觸點Tt時，ui,s等于ui；而當(dāng)油溫低于含蠟原油顯觸點時，ui,s為0[13]。本文中采用達(dá)西定律描述該束縛作用。此外，含蠟原油融化過程中還會吸收融化潛熱，融化潛熱表現(xiàn)為能量方程的內(nèi)熱源項ST。

1.2.2 非牛頓流體 當(dāng)油溫高于含蠟原油反常點Ta時，含蠟原油表現(xiàn)出牛頓流體的特性，其黏溫曲線可以由實驗數(shù)據(jù)擬合得到。當(dāng)油溫降低到含蠟原油反常點Ta以下時，含蠟原油表現(xiàn)出非牛頓性，本文中采用冪律方程描述含蠟原油的非牛頓性，如式（8）所示。此外，當(dāng)油溫降低到含蠟原油顯觸點Tt以下一定范圍時，含蠟原油的表觀黏度仍可以由冪律方程很好地描述[13]。因此，本文采用冪律方程描述含蠟原油的非牛頓特性。

1.2.3 浸入邊界法 由圖1 可知，罐內(nèi)含蠟原油區(qū)域是不規(guī)則的，而空氣層、保溫層和鋼板層都是規(guī)則的。如果采用非結(jié)構(gòu)化網(wǎng)格離散含蠟原油區(qū)域，勢必破壞地上浮頂油罐一體化求解過程，從而降低計算效率。基于此，本文引入浸入邊界法，在結(jié)構(gòu)化網(wǎng)格上對不規(guī)則的含蠟原油區(qū)域進(jìn)行離散，在一套結(jié)構(gòu)化網(wǎng)格上統(tǒng)一處理浮頂油罐內(nèi)規(guī)則區(qū)域（空氣層、鋼板層、保溫層）和非規(guī)則區(qū)域（含蠟原油區(qū)），從而確保了油罐部分的一體化求解過程。浸入邊界法的核心思想是：將固體區(qū)域簡化為固體邊界，將其對流體的影響以分布力和分布熱的方式作用于流體上。實際處理中，采用拉格朗日節(jié)點描述固體邊界，采用歐拉節(jié)點描述流體區(qū)域，二者之間通過狄拉克函數(shù)傳遞數(shù)據(jù)[23-25]。浸入邊界法的實施過程包括兩步：預(yù)測步和校正步。其中，預(yù)測步為動量方程的求解[式（6）、式（7）中fi和q為0]。本節(jié)重點介紹校正步的實施過程。式（9）給出了校正步的控制方程

離散控制方程式（9），并將式（10）代入，整理得

其中，A、b和X的表達(dá)式分別為

式中，h為網(wǎng)格尺寸，UjB為拉格朗日節(jié)點j的速度。求解式（12），獲得拉格朗日節(jié)點分布力FkB，并修正歐拉節(jié)點速度uj，如式（13）所示

如此就完成了速度場的校正。溫度場校正與此類似，具體過程可見文獻(xiàn)[23-25]。

1.3 邊界條件和初始條件

1.3.1 邊界條件 如圖1 所示，浮頂油罐融化模型包括地上油罐部分和罐底土壤部分，其邊界條件包括油罐邊界條件和土壤邊界條件。對于油罐，罐頂和罐壁為第三類邊界條件。油罐中心軸為軸對稱邊界條件。罐底與土壤層左上部為耦合邊界條件。對于土壤，土壤層左右邊界分別為軸對稱和絕熱邊界條件。下邊界溫度恒定為10℃。上邊界包括左側(cè)耦合邊界和右側(cè)第三類邊界。盤管壁面為第一類邊界條件。除此之外，本文邊界條件中包含了外界氣溫。計算中取融化初始時刻氣溫為1℃，融化12 h 時氣溫為12℃，其余時刻氣溫由余弦函數(shù)插值計算。

1.3.2 初始條件 浮頂油罐內(nèi)含蠟原油的融化初場是罐內(nèi)含蠟原油儲存溫度場、速度場等?？紤]到實際生產(chǎn)中不允許罐內(nèi)含蠟原油完全膠凝，本文以浮頂油罐內(nèi)含蠟原油儲存一定時間后的罐內(nèi)未完全融化的含蠟原油溫度場、速度場等作為融化初場。

2 數(shù)值計算方法

基于有限容積法，本文采用Fortran 語言開發(fā)了盤管組加熱的浮頂油罐內(nèi)含蠟原油融化和溫升數(shù)值計算程序，流程如圖2 所示。程序中非穩(wěn)態(tài)項采用全隱格式離散，時間步長取0.01 s。擴(kuò)散項和對流項分別采用中心差分格式和GAMMA 格式離散。采用同位網(wǎng)格SIMPLE 算法耦合壓力和速度[26-28]。本文的收斂標(biāo)準(zhǔn)為連續(xù)性方程余量小于1.0×10-6。本文中以1000 m3浮頂油罐為例，對罐內(nèi)含蠟原油湍流融化規(guī)律進(jìn)行研究。其中，1000 m3浮頂油罐的尺寸x1、x2、x3、x4、x5、r1、r2和r3分別為-5.0、0.01、8.01、8.02、8.52、6.0、6.11和11.11 m[25]。

本文對罐內(nèi)三種盤管組布置方式下（表1），罐內(nèi)含蠟原油的湍流融化過程進(jìn)行計算，研究盤管組布置方式對罐內(nèi)融化過程的影響，其中算例2 為對照算例?？紤]到盤管組內(nèi)加熱循環(huán)水的順利排放，工程設(shè)計中通常采用傾斜方式布置盤管組。

圖2 本文計算流程Fig.2 Flow chart of this research

表1 三種盤管組傾角Table 1 Three inclinations of heating coils at tank bottom

罐內(nèi)含蠟原油、空氣、蠟晶、鋼板、土壤和保溫層的熱物性如表2所示。

考慮融化過程中含蠟原油形態(tài)和流變性變化，含蠟原油黏度由式（14）計算。

表觀黏度μa由式（8）計算得到。其中，稠度系數(shù)K和流動特性指數(shù)n通過實驗數(shù)據(jù)擬合得到。本文中K和n的計算式為

表2 含蠟原油、空氣、蠟晶、鋼板、土壤和保溫層熱物性參數(shù)Table 2 Physical properties of waxy crude oil，air layer，wax crystal，steel layer，soil and insulating layer

析蠟量gs與油溫-T的關(guān)系由DSC 實驗數(shù)據(jù)擬合得到，如式（17）所示。

本文中以中國西部地區(qū)生產(chǎn)的含蠟原油為基礎(chǔ)開展研究，該含蠟原油的析蠟點Tw、反常點Ta和顯觸點Tt分別為27、25和23℃。

3 模型驗證

在計算之前，首先對本文模型進(jìn)行驗證。本節(jié)模型驗證主要包括浮頂油罐傳熱模型驗證、浸入邊界法驗證以及網(wǎng)格無關(guān)解和時間步長無關(guān)解的驗證三部分內(nèi)容。

需要說明的是，在模型驗證方面，目前文獻(xiàn)中還沒有實際浮頂油罐內(nèi)含蠟原油加熱融化的實驗數(shù)據(jù)。而在浮頂油罐內(nèi)含蠟原油溫降膠凝規(guī)律實驗研究方面，目前已經(jīng)有一些研究成果，且文獻(xiàn)中也給出了相關(guān)的實驗數(shù)據(jù)[29-30]。考慮到浮頂油罐內(nèi)含蠟原油的加熱融化過程與罐內(nèi)含蠟原油的溫降膠凝過程屬于同一個過程的不同發(fā)展方向，因此，本節(jié)將采用浮頂油罐溫降實驗數(shù)據(jù)對本文傳熱模型進(jìn)行驗證。本研究中還采用了浸入邊界法描述罐底盤管組與罐內(nèi)含蠟原油的耦合傳熱過程，本節(jié)將利用文獻(xiàn)中的數(shù)據(jù)對本文浸入邊界法進(jìn)行驗證。

3.1 浮頂油罐傳熱模型驗證

本節(jié)將利用文獻(xiàn)中的實驗數(shù)據(jù)對本文數(shù)值程序的準(zhǔn)確性進(jìn)行驗證。驗證中，本文選取中國西部地區(qū)和中國東部地區(qū)某油罐內(nèi)量油孔（r=38.7 m）處實測縱向油溫的平均值作為比較對象，從而確定模型的準(zhǔn)確性。其中，油罐半徑為40 m，罐內(nèi)液高分別為17.7 m和19.44 m，其他數(shù)據(jù)見文獻(xiàn)[29-30]。計算中，罐頂空氣層、含蠟原油區(qū)、罐底土壤層和保溫層的網(wǎng)格數(shù)（x×r）分別取97×241、160×241、82×302、266×13，總網(wǎng)格數(shù)為95504 個。時間步長取5.0 s。圖3 比較了實測結(jié)果與本文計算結(jié)果，從圖中可以看出，本文計算結(jié)果與實驗結(jié)果吻合良好，最大誤差控制在1℃以內(nèi)。

3.2 浸入邊界法驗證

利用本文程序?qū)Ψ角粌?nèi)圓管加熱的空氣自然對流進(jìn)行計算，通過將計算得到的圓管表面穩(wěn)態(tài)傳熱量與文獻(xiàn)結(jié)果[23,31]進(jìn)行對比來驗證本文程序的正確性，結(jié)果如表3 所示。其中，r和l分別為方腔內(nèi)圓管的半徑和方腔特征尺寸，其他參數(shù)見文獻(xiàn)[23,31]。從表中可以看出，本文結(jié)果與文獻(xiàn)結(jié)果吻合良好。

表3 方腔內(nèi)圓管恒壁溫加熱的穩(wěn)定熱流量Table 3 Steady heat flux of heating coil in a square cavity heated with constant temperature

3.3 網(wǎng)格無關(guān)解和時間步長無關(guān)解驗證

圖3 本文計算的量油孔縱向平均油溫與實測結(jié)果的比較Fig.3 Comparison of simulated and measured longitudinal average oil temperature at gauge hatch

計算之前，需要確定本研究中采用的網(wǎng)格數(shù)和時間步長。本節(jié)將對三組網(wǎng)格和三組時間步長下，1000 m3雙盤式浮頂油罐內(nèi)5 根盤管組加熱的含蠟原油融化過程（算例2）進(jìn)行計算，研究網(wǎng)格數(shù)和時間步長變化對罐內(nèi)含蠟原油融化過程的影響。計算之前，首先需要確定油罐融化初場。本節(jié)中，融化之前，罐內(nèi)核心區(qū)油溫約為24.55℃，罐底部約有0.485 m 厚的凝油層。圖4 給出了三組計算網(wǎng)格和三組計算時間步長下，罐內(nèi)含蠟原油平均溫升隨時間的變化曲線。從圖4(a)可以看出，當(dāng)網(wǎng)格總數(shù)分別為55312 個（含蠟原油區(qū)網(wǎng)格為254×178）和37224個（含蠟原油區(qū)網(wǎng)格為170×130）時，各時刻罐內(nèi)平均油溫近似相等，且二者與網(wǎng)格總數(shù)為28704個（含蠟原油區(qū)網(wǎng)格為142×110）的計算結(jié)果相差約3℃，因此，本文計算網(wǎng)格取37224 個（含蠟原油區(qū)網(wǎng)格為170×130）。同理，當(dāng)時間步長分別為0.01 s 和0.005 s 時，罐內(nèi)平均油溫相差不大，且二者與時間步長為0.02 s 的計算結(jié)果差距較大，如圖4(b)所示，因此本文計算中時間步長選取0.01 s。

圖4 罐內(nèi)平均油溫隨時間的變化曲線Fig.4 Variation of average oil temperature inside tank along with time

通過以上浮頂油罐傳熱模型和浸入邊界法驗證可以確定本文模型的準(zhǔn)確性。通過網(wǎng)格無關(guān)解和時間步長無關(guān)解驗證可以確認(rèn)本文空間步長和時間步長選擇的合理性。下文中將利用已經(jīng)驗證的模型以及篩選的空間步長和時間步長對三種盤管組布置條件下罐內(nèi)含蠟原油的湍流融化規(guī)律進(jìn)行研究。

4 結(jié)果及分析

本節(jié)主要研究不同盤管組布置方式對含蠟原油融化規(guī)律的影響，并對盤管組傳入浮頂油罐的熱量以及其分配情況進(jìn)行分析。在研究之前，本文對融化初始溫度場、速度場、析蠟情況等進(jìn)行計算。

4.1 融化初場確定

計算之前，首先需要確定融化初場。根據(jù)1.3.2節(jié)融化初場的設(shè)定規(guī)則，以化工區(qū)儲存一定時間的浮頂油罐溫度場、速度場等作為融化初場，如圖5所示。

圖5 浮頂油罐融化初場設(shè)定Fig.5 Initial fields for melting process of floating roof tank

分析圖5可知，融化開始之前，罐內(nèi)的流動主要發(fā)生在油罐中上部和罐頂空氣層，其中罐內(nèi)含蠟原油和罐頂空氣最大流速分別約為0.025 和0.015 m/s。融化開始后，在盤管組內(nèi)通入溫度恒定為80℃且快速流動的熱水來加熱融化罐內(nèi)含蠟原油，罐內(nèi)含蠟原油溫度場和流場分布如圖5(a)所示。在融化之前，罐內(nèi)膠凝含蠟原油主要分布在罐底部，最大析蠟量超過9.5%。隨著液高的增加，析蠟量逐漸減小。距罐底約0.4 m 范圍內(nèi)，蠟晶多孔介質(zhì)已經(jīng)形成，含蠟原油已經(jīng)膠凝（根據(jù)本文選取的含蠟原油，析蠟量超過2.67%時，含蠟原油將會膠凝），如圖5(b)所示。

4.2 管組加熱的浮頂油罐內(nèi)含蠟原油融化規(guī)律研究

圖6～圖8 分別給出了三組算例條件下1000 m3浮頂油罐內(nèi)含蠟原油湍流融化溫度場和速度矢量場的演化過程。分析可知，盤管組加熱的浮頂油罐內(nèi)含蠟原油融化過程可以分為兩個階段：融化階段和溫升階段。

（1）0～300 s，罐內(nèi)含蠟原油融化階段。在這一階段，盤管組傳入浮頂油罐的熱量主要用于融化油罐中上部懸浮的固態(tài)蠟晶和罐底部蠟晶多孔介質(zhì)，并加熱罐內(nèi)含蠟原油。在融化階段，罐內(nèi)含蠟原油的湍流自然對流過程依次經(jīng)歷了形成、發(fā)展并最終穩(wěn)定。具體分析可知，當(dāng)融化時間為0～30 s，熱量從盤管組傳入罐內(nèi)含蠟原油，加熱融化盤管組附近的含蠟原油和蠟晶。與此同時，被加熱的含蠟原油在浮升力作用下向上流動，湍流自然對流開始形成。當(dāng)融化時間為30～150 s，熱含蠟原油到達(dá)罐頂，并逐漸擴(kuò)展到整個油罐區(qū)域。在這一階段，罐內(nèi)湍流自然對流區(qū)域逐漸擴(kuò)散，強度逐漸增大，最大流速將達(dá)到0.7 m/s。此外，這一階段，在湍流自然對流的加熱作用下，不僅油罐中上部懸浮的固態(tài)蠟晶逐漸融化，而且罐底靠近罐中心的蠟晶多孔介質(zhì)在熱含蠟原油的沖刷下，逐漸解體并融化。當(dāng)融化時間為150～300 s 時，罐內(nèi)湍流自然對流過程逐漸穩(wěn)定，形成了從上至下順時針和逆時針的兩個主旋渦，如圖6(d)、(e)，圖7(d)、(e)和圖8(d)、(e)所示。而且，在湍流自然對流攪動下，罐內(nèi)油溫近似均勻，且溫度整體升高。而且，在這一階段，罐底靠近罐壁處的蠟晶多孔介質(zhì)經(jīng)過熱油的沖刷，最終完全融化。至此，罐內(nèi)含蠟原油融化過程完全結(jié)束。

（2）300 s 以后，罐內(nèi)含蠟原油溫升階段。這一階段，罐內(nèi)含蠟原油已經(jīng)完全融化，湍流自然對流已經(jīng)擴(kuò)展到整個油罐區(qū)域，且已經(jīng)形成了上下兩個主要的旋渦。此時，從盤管組傳入浮頂油罐的熱量主要用于加熱罐內(nèi)含蠟原油，并逐漸向罐頂空氣層、罐底土壤層以及外界大氣傳遞。而且，隨著罐內(nèi)油溫的升高，盤管組與含蠟原油之間的溫差逐漸減小，湍流自然對流強度逐漸減弱，表現(xiàn)為罐內(nèi)含蠟原油最大流速降低到0.35 m/s。

值得指出的是，在浮頂油罐加熱融化過程中，熱含蠟原油不斷沖擊罐頂鋼板層和罐底鋼板層，導(dǎo)致罐頂空氣層形成了一系列正逆相間的旋渦，如圖6(e)、(f)，圖7(e)、(f)和圖8(e)、(f)所示。而對于罐底土壤層，由于土壤區(qū)域較大、熱導(dǎo)率較小，而且加熱時間較短（1200 s），傳熱方式為導(dǎo)熱，因此罐底土壤層溫度場未出現(xiàn)明顯的變化。

圖6 算例1條件下浮頂油罐內(nèi)融化溫度場和速度矢量場分布Fig.6 Temperature and velocity fields in tank during melting process of case 1

從圖6(b)、圖7(b)、圖8(b)可以看出，在融化初期，當(dāng)盤管組傾角θ逐漸增大時，被加熱的含蠟原油上浮過程將傾向于盤管組正上方。分析其原因可知，盤管組左右兩側(cè)分別為罐中心和罐壁（圖1），且罐中心體積明顯小于罐壁處，因此，在上升熱油的加熱下，罐中心附近含蠟原油的溫升速率將明顯快于罐壁處。當(dāng)盤管組非水平布置時（θ=-5.71°或5.71°），被盤管組加熱的含蠟原油將沿盤管組向上流動。在以上兩方面的作用下，當(dāng)θ較小時，上升熱油將向左偏移[圖6(b)、圖7(b)]，當(dāng)θ較大時，熱油垂直上升[圖8(b)]。而且，當(dāng)熱油到達(dá)罐頂，并在罐頂堆積時，造成上升熱油溫差減小，阻力增大，因此在靠近罐壁的盤管上方將形成上升油流[圖6(b)、圖8(b)]，融化過程加速。

圖9給出了罐內(nèi)平均油溫及融化液相分?jǐn)?shù)隨時間的變化曲線。分析可知，罐內(nèi)平均油溫和融化液相分?jǐn)?shù)隨時間持續(xù)增大，且增長率逐漸減小。結(jié)合前文的分析可知，在融化初期（0～70 s），算例1 條件下，融化區(qū)域主要位于罐中心處，而算例2 和算例3對應(yīng)的融化區(qū)域位于盤管組上方，因此3 個算例表現(xiàn)出相似的融化規(guī)律，即三者液相分?jǐn)?shù)和平均油溫相差不大。隨著融化過程的進(jìn)行（70 s 以后），θ=5.71°時（算例3），靠近罐壁處的盤管組上方率先出現(xiàn)了上升油流[圖8(c)]，罐內(nèi)含蠟原油融化效率顯著提高，且該過程一直持續(xù)到融化過程結(jié)束。對于溫升過程而言，當(dāng)融化時間超過300 s 后，罐內(nèi)含蠟原油已經(jīng)完全融化，湍流自然對流過程已經(jīng)穩(wěn)定，此時，算例1 和算例2 對應(yīng)的罐內(nèi)含蠟原油流動過程較為相似，尤其是靠近罐底部的逆時針主旋渦[圖6(e)、圖7(e)]，因此二者表現(xiàn)出相似的溫升規(guī)律，且均小于算例3條件下的結(jié)果。

圖7 算例2條件下浮頂油罐內(nèi)融化溫度場和速度矢量場分布Fig.7 Temperature and velocity fields in tank during melting process of case 2

綜上所述，在本文計算條件下，與算例1和算例2相比，算例3表現(xiàn)出較高的融化效率和較快的溫升速率。

4.3 盤管組熱流量及其分配情況分析

盤管組傳入浮頂油罐的熱量是驅(qū)動罐內(nèi)含蠟原油融化和溫升的動力，本節(jié)對盤管組的熱流量及其分配情況進(jìn)行研究。式（18）給出了盤管組總熱流量的計算式。

式中，qtot、qht、qtr和qpc分別為盤管組總熱流量，從罐頂、罐底和罐壁鋼板層傳出的熱流量，單位時間內(nèi)含蠟原油因為溫升和融化而吸收的熱量，MW。盤管組總傳熱量Qtot是從融化開始到第n個時間步內(nèi)盤管組傳入油罐的熱量，由式（19）計算，單位為MJ。

圖10 給出了盤管組的總熱流量及其分配。圖10(a)為算例2 條件下，罐內(nèi)盤管組熱流量及其分配情況。從圖中可以看出，時層熱流量（qtot）經(jīng)歷了先快速減小，然后小幅度反彈，最后緩慢波動式減小的過程，這與浮頂油罐內(nèi)層流融化類似[13]。分析其原因，融化初始時刻，盤管組與周圍含蠟原油的溫差高達(dá)55.8℃，熱量以導(dǎo)熱方式快速傳入罐內(nèi)，最大熱流量達(dá)到100 MW。隨著傳熱過程的進(jìn)行，盤管組與周圍含蠟原油的溫差快速減小，造成熱流量快速減小到約30 MW。隨著融化過程的進(jìn)行，罐內(nèi)湍流自然對流過程逐漸形成，熱量以對流的方式傳遞到整個油罐區(qū)域，熱流量出現(xiàn)了波動式的小幅反彈，最大熱流量約40 MW。隨著罐內(nèi)油溫的升高，盤管組與含蠟原油的溫差逐漸減小，自然對流強度逐漸減小，熱流量緩慢波動式減小。就熱流量分配而言，當(dāng)融化時間在0～200 s 時，熱流量主要用于加熱和融化罐內(nèi)含蠟原油。當(dāng)融化時間超過200 s后，熱流量基本全部用于加熱含蠟原油。此外，從含蠟原油傳入罐頂空氣層、罐底土壤層和罐壁保溫層的平均熱流量約400 W。

圖8 算例3條件下浮頂油罐內(nèi)融化溫度場和速度矢量場分布Fig.8 Temperature and velocity fields in tank during melting process of case 3

圖9 浮頂油罐內(nèi)含蠟原油融化特性Fig.9 Melting characteristic of waxy crude oil in floating roof tank

圖10 盤管組熱流量隨時間的變化Fig.10 Variations of heat flux of coils along with time

如前文所述，盤管組布置方式將對罐內(nèi)含蠟原油融化過程產(chǎn)生較大的影響。融化熱流量表現(xiàn)出“山峰狀”變化的特點，即先快速增大，然后波動式減小。這是因為自然對流形成后，熱油上浮到達(dá)罐頂之前，融化區(qū)域快速擴(kuò)大，熱流量快速增大。當(dāng)熱油上浮到達(dá)罐頂之后，自然對流在罐內(nèi)擴(kuò)展，罐內(nèi)流動過程變得非常復(fù)雜[圖6(d)、圖8(d)]，此時罐內(nèi)中上部油溫快速升高，最終導(dǎo)致熱流量波動式減小。罐內(nèi)湍流自然對流的發(fā)展對時層熱流量和總熱量具有決定性的影響。結(jié)合前文分析可知，由于融化初期和溫升階段罐內(nèi)含蠟原油的流動過程類似，因此在這兩個階段時層熱流量變化過程類似，而在融化中期，不同盤管組布置方式下，罐內(nèi)含蠟原油自然對流過程非常復(fù)雜，時層熱流量出現(xiàn)了一些不同。就總熱量而言，由于融化過程中，罐內(nèi)流動過程存在差異，導(dǎo)致不同盤管組布置方式下，盤管組總熱量出現(xiàn)了差異，并在整個溫升過程中維持不變（總熱量差異約1000 MJ）。此外，從圖10 中還可以看出，各時層盤管組的熱流量處于劇烈波動之中，這是因為融化和溫升過程中，罐內(nèi)含蠟原油始終處于強湍流自然對流過程的緣故。

綜上所述，在算例3對應(yīng)的盤管組布置方式下，相同時刻盤管組總熱量較大，傳熱效率較高。因此，為了提高浮頂油罐融化和加熱效率，宜選用算例3盤管組布置方式。

5 結(jié) 論

利用在盤管組中通入熱水來加熱融化罐內(nèi)含蠟原油是最常用的油罐加熱方式之一，而盤管組布置方式對于罐內(nèi)含蠟原油加熱融化效率等具有重要的影響?；诖?，本研究在考慮含蠟原油融化過程中形態(tài)和物性變化的基礎(chǔ)上，建立了盤管組加熱的浮頂油罐內(nèi)含蠟原油融化物理數(shù)學(xué)模型。為提高計算效率，本文采用浸入邊界法在結(jié)構(gòu)化網(wǎng)格上對罐內(nèi)規(guī)則區(qū)域和非規(guī)則區(qū)域進(jìn)行統(tǒng)一處理，實現(xiàn)了油罐融化過程的一體化耦合求解。以實際1000 m3雙盤式浮頂油罐為例，研究了不同盤管組布置方式下，罐內(nèi)含蠟原油的融化和溫升規(guī)律，分析了盤管組布置方式對浮頂油罐融化溫升過程的影響，得到以下結(jié)論。

（1）盤管組加熱的浮頂油罐內(nèi)含蠟原油融化過程分為兩個階段：融化階段和溫升階段。在融化階段，罐內(nèi)湍流自然對流依次經(jīng)歷形成、發(fā)展和穩(wěn)定過程。而且，在形成期和發(fā)展期，熱油在罐內(nèi)上浮、擴(kuò)展，流動過程相對復(fù)雜且混亂。此時，融化主要發(fā)生在油罐中上部和罐底靠罐中心處。在穩(wěn)定期，罐內(nèi)含蠟原油流動過程趨于穩(wěn)定，盤管組傳入的熱量主要用于加熱罐底靠近罐壁處的蠟晶多孔介質(zhì)。在溫升階段，罐內(nèi)形成了從上至下順時針和逆時針的兩個穩(wěn)定的旋渦。而且，在這兩個旋渦的攪動下，罐內(nèi)油溫近似均勻，且溫度整體上升。

（2）時層熱流量表現(xiàn)出先快速減小，然后小幅度反彈，最后緩慢波動式減小的特點。融化吸熱流量隨時間經(jīng)歷“山峰狀”的變化過程，最大融化吸熱流量達(dá)到約10 MW。三種盤管組布置方式下，盤管組總熱量隨時間的變化過程相似。在融化后期，算例3的總熱量比算例1高約1000 MJ。

（3）本文算例中，算例3 表現(xiàn)出較快的溫升速率、較大的傳熱量以及較高的融化效率。

符 號 說 明

cp——比定壓熱容，J/(kg·℃)

d——計算位置到最近邊界的距離，m

——拉格朗日節(jié)點的分布力向量，kg/(m3·s2)

f——歐拉節(jié)點的分布力向量，kg/(m2·s2)

fi——由浸入邊界法引起的分布力密度，kg/(m2·s2)

gi——重力加速度向量，m/s2

gs——析蠟分?jǐn)?shù)

ΔH——相變潛熱，J/kg

K——稠度系數(shù)

Kd——蠟晶多孔介質(zhì)滲透率，m2

L——融化總潛熱，J/kg

Ls——亞格子混合長度，m

n——流動特性指數(shù)

Prt——湍流Prandtl數(shù)

q——由浸入邊界法引起的分布熱密度，kg·℃/(m2·s)

qht,qpc,qtr,qtot——分別為向外界傳出的熱流量、融化吸熱量、溫升吸熱量、盤管組總熱流量，MW

qj——大渦模擬引起的熱通量，kg·℃/(m2·s)

Ra——Rayleigh數(shù)

rp1,rp2,rp3,rp4,rp5——從左向右盤管中心徑向坐標(biāo)，m

r0——盤管半徑，m

δr——相鄰盤管間距，m

ST——能量方程源項，J/(m3·s)

Sui——動量方程源項，kg/(m2·s2)

Δsi——盤管表面表面第i段弧長，m

Ta,Tt,Tw——分別為含蠟原油反常點、顯觸點和析蠟點，℃

T0——盤管壁面溫度，℃

t——時間，s

xp1,xp2,xp3,xp4,xp5——從左向右盤管中心軸向坐標(biāo)，m

β——體積膨脹系數(shù)，℃-1

Δ——過濾尺寸，m

θ——盤管組所在平面與r軸正方向夾角，(°)

κ——馮卡門常數(shù)，1

λ——熱導(dǎo)率，W/(m·℃)

μ,μt,μa——分別為動力黏度、亞格子黏度和表觀黏度，Pa·s

ρ——密度，kg/m3

τij——剪應(yīng)力，kg/(m·s2)