王敏，李敬法，張欣雨，宇波

1 油氣管道輸送安全國家工程實(shí)驗(yàn)室/城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室/中國石油大學(xué)(北京)，北京 102249 2 中國石油化工集團(tuán)國際石油勘探開發(fā)公司，北京 100029 3 機(jī)械工程學(xué)院，北京石油化工學(xué)院，北京 102617

單盤式浮頂油罐內(nèi)含蠟原油溫降規(guī)律數(shù)值計(jì)算研究

王敏1，李敬法1，張欣雨2，宇波3*

1 油氣管道輸送安全國家工程實(shí)驗(yàn)室/城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室/中國石油大學(xué)(北京)，北京 102249 2 中國石油化工集團(tuán)國際石油勘探開發(fā)公司，北京 100029 3 機(jī)械工程學(xué)院，北京石油化工學(xué)院，北京 102617

建立了考慮含蠟原油非牛頓性和析蠟相變過程的單盤式浮頂油罐內(nèi)含蠟原油溫降物理數(shù)學(xué)模型，并發(fā)展了一體化耦合求解方法。模型中，采用冪律方程描述罐內(nèi)含蠟原油的非牛頓特性；采用湍流大渦模擬方法計(jì)算罐內(nèi)含蠟原油的湍流流動過程；采用焓多孔介質(zhì)理論計(jì)算罐內(nèi)含蠟原油的析蠟相變過程，并跟蹤罐內(nèi)含蠟原油的相變界面。利用文獻(xiàn)中的結(jié)果對本文模型進(jìn)行驗(yàn)證。基于驗(yàn)證的程序，對實(shí)際單盤式浮頂油罐內(nèi)含蠟原油溫降過程進(jìn)行計(jì)算，研究罐內(nèi)含蠟原油溫降過程的演化規(guī)律，并對罐頂、罐底凝油層的增長以及罐頂、罐底和罐壁處傳熱量的變化規(guī)律進(jìn)行研究。研究結(jié)果可為單盤式浮頂油罐的設(shè)計(jì)和生產(chǎn)方案的制定提供一定的指導(dǎo)。

單盤式浮頂油罐；溫降規(guī)律；非牛頓性；析蠟相變；大渦模擬

0 背景

中國所產(chǎn)原油大多為含蠟原油。含蠟原油具有易凝、高黏、流變性復(fù)雜等特點(diǎn)，在實(shí)際生產(chǎn)儲存中通常需要加熱處理。單盤式浮頂油罐是常用的含蠟原油儲罐之一[1]。含蠟原油在罐內(nèi)儲存時會向外界環(huán)境放熱而被冷卻。當(dāng)油溫降低到含蠟原油析蠟點(diǎn)時，含蠟原油開始析蠟；當(dāng)析蠟分?jǐn)?shù)達(dá)到原油總質(zhì)量的2%～3%時，含蠟原油將會膠凝[2]，嚴(yán)重影響實(shí)際生產(chǎn)，甚至造成生產(chǎn)事故。研究單盤式浮頂油罐內(nèi)含蠟原油的溫降過程，準(zhǔn)確把握罐內(nèi)含蠟原油的溫降規(guī)律，對科學(xué)設(shè)計(jì)浮頂油罐的保溫結(jié)構(gòu)、合理制定浮頂油罐的周轉(zhuǎn)方案具有重要意義。

目前，研究油罐內(nèi)含蠟原油溫降規(guī)律的方法主要有兩種：實(shí)驗(yàn)研究[3-5]和數(shù)值模擬研究[6-9]，其中數(shù)值模擬因?yàn)橥顿Y小、研究周期短以及研究結(jié)果全面等優(yōu)點(diǎn)成為近年來最主要的研究手段。在數(shù)值模擬研究方面，Lin等[6]采用數(shù)值計(jì)算方法對高徑比為1/3～1，Pr數(shù)為1～1 000、Ra數(shù)為6×106～6×1010的油罐內(nèi)溫降規(guī)律進(jìn)行研究。在此基礎(chǔ)上，采用尺度分析方法擬合了無量綱溫度與時間、儲罐高徑比、Pa數(shù)和Ra數(shù)的關(guān)系式。Oliveski等[7]采用實(shí)驗(yàn)方法和數(shù)值計(jì)算方法對Ra數(shù)分別為3.1×109和4.3×109的實(shí)驗(yàn)油罐內(nèi)原油溫降15 h的原油溫度場和速度場進(jìn)行計(jì)算，給出了不同溫降時間下，油罐豎直中心線上的油溫分布以及罐頂和罐壁處原油的速度矢量場分布。Oliveski[8]采用有限容積法對不同容積的實(shí)驗(yàn)罐內(nèi)原油的溫降過程展開研究，并采用量綱分析法回歸了油罐的平均Nu數(shù)與儲罐無量綱高徑比、無量綱總傳熱系數(shù)以及Ra數(shù)和Pr數(shù)的換熱關(guān)聯(lián)式。李旺[9]采用大渦模擬方法對10萬方雙盤式浮頂油罐內(nèi)原油溫降過程進(jìn)行計(jì)算，并分析了浮盤類型、保溫層厚度、液位、油品物性、太陽輻射以及風(fēng)速對罐內(nèi)原油溫降過程的影響。研究發(fā)現(xiàn)，浮盤類型、保溫層厚度、液位和太陽輻射對油罐溫降過程具有較大影響，而風(fēng)速和油品物性對油罐溫降過程影響較小。

雖然國內(nèi)外學(xué)者已經(jīng)對單盤式浮頂油罐內(nèi)原油的溫降規(guī)律開展了一定的研究，但這些研究通常將原油簡化為牛頓流體，沒有考慮含蠟原油的析蠟相變過程。而且研究中往往只給出油罐內(nèi)某些特定位置處油溫隨溫降時間的變化。就作者查閱的文獻(xiàn)而言，目前還沒有見到全面研究罐內(nèi)含蠟原油溫降演化規(guī)律以及罐內(nèi)含蠟原油湍流流動與傳熱特性的研究成果。基于此，本文建立了考慮含蠟原油非牛頓性和析蠟相變過程的單盤式浮頂油罐內(nèi)含蠟原油湍流溫降物理數(shù)學(xué)模型，并發(fā)展了一體化耦合求解技術(shù)。通過對實(shí)際單盤式浮頂油罐內(nèi)原油溫降過程進(jìn)行計(jì)算，研究罐內(nèi)含蠟原油溫度場及流場的演化規(guī)律，并對凝油層增長以及傳熱量變化進(jìn)行分析。

1 物理數(shù)學(xué)模型及定解條件

1.1 物理模型

單盤式浮頂油罐包括：鋼板層、保溫層和含蠟原油，其傳熱過程包括液固(含蠟原油與鋼板層)，固固(鋼板層和保溫層)耦合傳熱以及罐體與外界大氣之間強(qiáng)制對流換熱和油罐與罐底土壤層之間的導(dǎo)熱。因此，單盤式浮頂油罐的傳熱過程非常復(fù)雜，計(jì)算難度較大。

對含蠟原油溫降過程分析發(fā)現(xiàn)，初始時刻含蠟原油溫度較高，含蠟原油以純液態(tài)形式存在，且表現(xiàn)出牛頓流體特性；當(dāng)油溫開始降低到含蠟原油析蠟點(diǎn)Tw(含蠟原油中開始有蠟晶析出時的溫度)以下時，含蠟原油中開始有蠟晶析出，固液分散體系開始形成。此時含蠟原油仍然表現(xiàn)出牛頓流體特性；當(dāng)油溫降低到含蠟原油反常點(diǎn)Ta(牛頓流體與非牛頓流體的轉(zhuǎn)換溫度)以下時，含蠟原油開始表現(xiàn)出非牛頓流體特性，此時含蠟原油仍以固液分散體系形式存在；當(dāng)油溫繼續(xù)降低到含蠟原油顯觸點(diǎn)Tt(含蠟原油顯現(xiàn)觸變性的最高溫度)以下時，含蠟原油中析出的蠟晶開始發(fā)生交聯(lián)，形成蠟晶多孔介質(zhì)結(jié)構(gòu)[2]。表1給出了溫降過程中含蠟原油的狀態(tài)以及流變性的變化過程。

由于實(shí)際單盤式浮頂油罐的尺寸較大，罐內(nèi)含蠟原油傳熱方式為湍流自然對流(實(shí)際1 000方單盤式浮頂油罐半徑為6 m，當(dāng)溫差為30 ℃時，罐內(nèi)含蠟原油Ra數(shù)高達(dá)2.0×1013)，因此要對這樣的三維實(shí)際油罐進(jìn)行數(shù)值計(jì)算，其計(jì)算難度和計(jì)算量都將非常大，計(jì)算非常耗時?？紤]到單盤式浮頂油罐的對稱性，本文將三維油罐簡化為二維油罐，如圖1所示。采用湍流大渦模擬方法對考慮含蠟原油的非牛頓性和析蠟相變過程的單盤式浮頂油罐內(nèi)含蠟原油湍流溫降過程進(jìn)行計(jì)算，研究罐內(nèi)含蠟原油的溫降規(guī)律。

表1 溫降過程中含蠟原油的狀態(tài)和流變性變化Table 1 The state and rheological behavior of waxy crude oil in temperature drop process

圖1 二維單盤式浮頂油罐物理模型Fig. 1 Sketch map of single-plate fl oating roof oil tank

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

對浮頂油罐內(nèi)含蠟原油溫降過程分析發(fā)現(xiàn)，其本質(zhì)是一個非穩(wěn)態(tài)自然對流過程，其中涉及相變、固液相作用、非牛頓性、流固耦合以及湍流。結(jié)合1.1節(jié)建立的單盤式浮頂油罐溫降物理模型，建立考慮含蠟原油非牛頓性[2,10-11]和析蠟相變過程[12-14]的單盤式浮頂油罐內(nèi)含蠟原油湍流溫降數(shù)學(xué)模型?？紤]到物理模型中涉及到流固耦合，為了實(shí)現(xiàn)一體化耦合求解，將固體(鋼板層和保溫層)區(qū)域設(shè)定為黏度為無窮大的流體，其余參數(shù)按照實(shí)際參數(shù)取值。此外，本文采用大渦模擬方法[15]計(jì)算罐內(nèi)含蠟原油的湍流過程。二維圓柱坐標(biāo)系下單盤式浮頂油罐溫降控制方程如式(1)～(4)所示：

其中，ρ、β和cp分別為密度、體積膨脹系數(shù)和定壓熱容。和分別為采用大渦模擬方法過濾后x方向和r方向的速度分量、壓力和溫度。sux和sur為蠟晶多孔介質(zhì)對含蠟原油的束縛作用，本文稱之為達(dá)西源項(xiàng)。sh為含蠟原油析蠟潛熱項(xiàng)。含蠟原油的析蠟相變處理方法將在下文中詳細(xì)介紹。μeff和λeff分別為有效黏度和有效導(dǎo)熱系數(shù)，本文中將分別由式(5)和式(6)計(jì)算：

其中，μt和αt分別為亞格子渦黏系數(shù)和亞格子渦擴(kuò)散系數(shù)。亞格子渦黏系數(shù)的構(gòu)造是大渦模擬方法的關(guān)鍵，本文中將采用Smagorinshy-Lilly模型[16]構(gòu)造，具體表達(dá)式如式(7)所示：

式中，Ls為亞格子混合長度，其定義式為Ls=min(κd, CsΔ)[16]?；旌祥L度的計(jì)算中考慮了邊界效應(yīng)，對近壁區(qū)做了壁面修正；κ為馮卡門常數(shù)，本文中取0.41；d為計(jì)算位置距最近邊界的距離；Cs為模型常數(shù)，本文中取0.1；Δ為過濾尺寸。

對于亞格子渦擴(kuò)散系數(shù)，本文中采用式(8)計(jì)算。其中，Prt取0.87[16]。

1.2.1 非牛頓性

當(dāng)油溫降低到含蠟原油的反常點(diǎn)以下時，含蠟原油為假塑性流體，其剪應(yīng)力與剪切率的關(guān)系可由冪律方程很好地描述。此外，當(dāng)油溫繼續(xù)降低到含蠟原油顯觸點(diǎn)以下一定范圍內(nèi)時，冪律方程仍能很好的描述觸變性含蠟原油的剪應(yīng)力與剪切率的關(guān)系[2]，因此選擇冪律方程描述含蠟原油的非牛頓性。式(9)給出了含蠟原油的表觀黏度μa與剪切率˙之間的關(guān)系式。

其中，K為稠度系數(shù)，單位為Pa·sn；n為流動特性指數(shù)。K和n與油溫T的關(guān)系由實(shí)驗(yàn)數(shù)據(jù)擬合得到。為罐內(nèi)含蠟原油的剪切率，由式(10)計(jì)算[17]。

1.2.2 析蠟相變過程

當(dāng)油溫降低到含蠟原油析蠟點(diǎn)以下時，含蠟原油開始析蠟。此時，析出的蠟晶漂浮在含蠟原油中，隨著含蠟原油一起流動。當(dāng)油溫降低到含蠟原油顯觸點(diǎn)以下時，析出的蠟晶開始相互交聯(lián)，形成蠟晶多孔介質(zhì)結(jié)構(gòu)，并對包裹在其孔隙中的含蠟原油產(chǎn)生束縛。隨著油溫的繼續(xù)降低，蠟晶多孔介質(zhì)孔隙度逐漸減小，束縛作用逐漸增強(qiáng)。本文采用達(dá)西定律描述這種束縛作用，達(dá)西源項(xiàng)的表達(dá)式如式(11)和式(12)所示。

其中，K0為滲透率常數(shù)，本文中取10-7[14,18]。gs為析蠟分?jǐn)?shù)，即某一溫度下，含蠟原油中析出的蠟晶質(zhì)量占原油總質(zhì)量的比例，由DSC實(shí)驗(yàn)確定。

含蠟原油在析蠟相變過程中還將釋放相變潛熱，相變潛熱表現(xiàn)為能量方程中的源項(xiàng)sh，如式(14)所示：

式中，ΔH為相變潛熱。本文中ΔH=(1-gs)L，L為含蠟原油的析蠟相變總潛熱。

1.3 邊界條件和初始條件

1.3.1 邊界條件

單盤式浮頂油罐罐頂、罐壁以及浮頂油罐右側(cè)土壤表層均為第三類邊界條件，其對流換熱系數(shù)為25.63 W/(m2·℃)[2]。油罐中心軸為對稱邊界，土壤層右邊界為絕熱邊界，下邊界為恒溫邊界，溫度取10 ℃。此外，在油罐邊界條件設(shè)定中，包含了外界氣溫，本文以中國某地區(qū)11月到12月的實(shí)測氣溫為基礎(chǔ)數(shù)據(jù)來計(jì)算。

1.3.2 初始條件

本文中規(guī)定初始時刻罐內(nèi)含蠟原油、空氣層、鋼板層和保溫層的溫度均為30 ℃。地表溫度為初始時刻的外界氣溫，為1 ℃。不同深度處的土壤溫度由土壤恒溫層和地表溫度線性插值得到。

2 數(shù)值計(jì)算方法及模型驗(yàn)證

2.1 數(shù)值計(jì)算方法

本文采用有限容積法對控制方程進(jìn)行離散，其中對流項(xiàng)采用GAMMA格式[19]離散，擴(kuò)散項(xiàng)采用中心差分格式離散，時間項(xiàng)采用一階全隱格式離散。采用基于非均分同位網(wǎng)格的SIMPLE算法耦合壓力和速度，求解器選用強(qiáng)隱方法[20]。為提高計(jì)算速度，本文采用多重網(wǎng)格方法求解罐底土壤溫度場，采用一體化求解方法對油罐部分(包括罐頂空氣層和鋼板層、罐底鋼板層和罐壁保溫層)溫度場和流場等進(jìn)行計(jì)算，再與罐底土壤溫度場進(jìn)行耦合。如此反復(fù)，直至兩部分均達(dá)到收斂標(biāo)準(zhǔn)為止。本文中時間步長取5 s。此時收斂標(biāo)準(zhǔn)為：油罐部分質(zhì)量不平衡量小于10-4，罐底土壤層導(dǎo)熱方程余量小于10-6。

另一個關(guān)鍵問題是相變潛熱ΔH的更新。在數(shù)值計(jì)算中，每個時層內(nèi)相變潛熱根據(jù)各時層內(nèi)油溫進(jìn)行更新，相變潛熱能否正確更新會直接影響計(jì)算過程的穩(wěn)定性和收斂性。本文采用式(15)更新相變潛熱[12]。

其中，α為亞松弛因子，本文中取0.2；m和m-1分別為當(dāng)前迭代步和上一迭代步；為相變潛熱的反函數(shù)[12]，即含蠟原油溫度。根據(jù)前面的分析可知，當(dāng)油溫低于含蠟原油析蠟點(diǎn)時，ΔH=(1.0-gs)L，析蠟分?jǐn)?shù)gs與油溫之間的關(guān)系可以由DSC實(shí)驗(yàn)測得，由此，可以得到的計(jì)算式。

根據(jù)建立的單盤式浮頂油罐溫降物理數(shù)學(xué)模型，本文采用Fortran語言編寫了計(jì)算程序，相應(yīng)的程序設(shè)計(jì)流程圖如圖2所示。

以實(shí)際1 000方的單盤式浮頂油罐為例，對其溫降規(guī)律展開研究。其中1 000方單盤式浮頂油罐的尺寸x1、x2、x3、x4、r1、r2和r3分別為-5.0 m、0.01 m、8.01 m、8.02 m、6.0 m、6.11 m和11.11 m[1]。罐內(nèi)含蠟原油、空氣層、鋼板層、土壤層和保溫層的物性如表2。

考慮非牛頓性和析蠟相變過程的含蠟原油黏度由式(16)計(jì)算。

含蠟原油的表觀黏度μa由式(9)計(jì)算得到。其中，稠度系數(shù)K和流動特性指數(shù)n由實(shí)驗(yàn)數(shù)據(jù)擬合得到，本文中的計(jì)算式見式(17)和(18)。

對含蠟原油進(jìn)行DSC實(shí)驗(yàn)，擬合得到析蠟分?jǐn)?shù)gs與油溫T的關(guān)系見式(19)。

根據(jù)析蠟潛熱ΔH、析蠟分?jǐn)?shù)gs及油溫的關(guān)系，可以推導(dǎo)出f-1(ΔH)的表達(dá)式如式(20)。

其中，含蠟原油的析蠟點(diǎn)Tw、反常點(diǎn)Ta和顯觸點(diǎn)Tt分別為27 ℃、25 ℃和23 ℃。

圖2 程序設(shè)計(jì)流程圖Fig. 2 Flow chart of the calculation program

2.2 大渦模擬方法驗(yàn)證

為了驗(yàn)證大渦模擬方法的準(zhǔn)確性，本文在102×102的非均分網(wǎng)格上對直角坐標(biāo)系方腔內(nèi)Pr數(shù)和Ra數(shù)分別為0.71和1.58×109的湍流自然對流過程進(jìn)行計(jì)算，將計(jì)算結(jié)果與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果[21]進(jìn)行比較。圖3給出了由本文大渦模擬方法計(jì)算的方腔水平中心線上無量綱溫度Θ和無量綱速度Uy與文獻(xiàn)結(jié)果的比較，從圖中可以看出，本文計(jì)算結(jié)果與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果吻合良好。

2.3 非牛頓流體模型驗(yàn)證

為了驗(yàn)證非牛頓流體模型的準(zhǔn)確性，在102×102的均分網(wǎng)格上，對Pr數(shù)為100，Ra數(shù)分別為104和105的方腔非牛頓流體穩(wěn)態(tài)自然對流過程進(jìn)行計(jì)算，并將計(jì)算得到的熱邊界平均Nu數(shù)與文獻(xiàn)中的結(jié)果[10,11]進(jìn)行比較，結(jié)果如表3所示。從表中可以看出，本文計(jì)算結(jié)果與文獻(xiàn)結(jié)果吻合良好。

2.4 相變模型驗(yàn)證

為了驗(yàn)證相變模型的正確性，本文采用流固耦合方法對半徑為31.9 mm，高為59.8 mm，壁面和底部包裹有厚度分別為2.85 mm的聚碳酸酯和5.97 mm聚丙烯酸固體保溫材料的圓柱形區(qū)域內(nèi)正二十烷的融化過程進(jìn)行計(jì)算。計(jì)算中，正二十烷的初始溫度為23 ℃，相變溫度為36.4 ℃，頂部絕熱，壁面和底部溫度分別為45 ℃和32 ℃，物性等其他數(shù)據(jù)參見Jones等文獻(xiàn)[13]。圖4(a)給出了融化9601.6 s后區(qū)域內(nèi)正二十烷的融化界面位置以及液態(tài)正二十烷中的溫度場和速度矢量場分布。圖4(b)給出了不同融化時間下，區(qū)域內(nèi)正二十烷融化界面位置與實(shí)驗(yàn)結(jié)果的對比，從圖中可以看出，本文計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果[13]吻合良好。

表2 含蠟原油、空氣層、鋼板層、土壤層和保溫層的物性Table 2 Thermo-physical properties of different materials

3 單盤式浮頂油罐內(nèi)含蠟原油溫降規(guī)律研究

3.1 單盤式浮頂油罐溫度場演化

圖5給出了不同溫降時刻，1000 m3單盤式浮頂油罐內(nèi)含蠟原油溫度場和速度矢量場分布。對圖5進(jìn)行分析可知，初始時刻，罐底部含蠟原油溫度為30 ℃，土壤表層溫度僅為1 ℃，熱量快速從罐內(nèi)含蠟原油傳入罐底土壤層，罐底含蠟原油中形成了層狀分布的溫度場，如圖5(a)。而且，在這一階段，罐頂部含蠟原油將通過灌頂鋼板層向外界大氣傳熱，同時罐頂部含蠟原油損失的熱量將通過浮升力，以自然對流方式從油罐中部獲得補(bǔ)充，湍流自然對流過程開始形成。此時，流動過程僅限于油罐右上部，如圖6(a)。隨著溫降過程的繼續(xù)進(jìn)行，罐內(nèi)含蠟原油湍流流動區(qū)域逐漸擴(kuò)展，溫降100 h時，湍流自然對流區(qū)域逐漸擴(kuò)展整個油罐區(qū)域，罐內(nèi)x方向最大流速達(dá)到0.03 m/s，如圖5(b)，罐底部層狀分布的溫度場消失，罐內(nèi)油溫近似均勻，且整體發(fā)生溫降，如圖5(b)～圖5(c)。當(dāng)罐內(nèi)油溫降低到含蠟原油顯觸點(diǎn)(23 ℃)以下時，罐內(nèi)含蠟原油中析出的蠟晶開始相互交聯(lián)，形成蠟晶多孔介質(zhì)結(jié)構(gòu)，并對包裹在其孔隙中的含蠟原油產(chǎn)生束縛，傳熱方式由固液分散體系自然對流轉(zhuǎn)變?yōu)橄灳Ф嗫捉橘|(zhì)傳熱，類似于導(dǎo)熱，此時，罐內(nèi)含蠟原油流動過程明顯減弱。溫降300 h時，罐內(nèi)x方向最大流速僅約為4.0×10-4m/s，如圖5(d)。而且在罐頂、罐底和罐壁冷量的共同作用下，罐內(nèi)出現(xiàn)了明顯的冷熱油界面，如圖5(d)～圖5(f)。

圖3 大渦模擬方法驗(yàn)證Fig. 3 Validation of LES method

表3 Pr=100時，本文計(jì)算的方腔熱邊界平均Nu數(shù)與文獻(xiàn)結(jié)果的對比Table 3 Comparison of Nu at the high temperature boundary with the numerical results obtained by this paper, Matin[10]and Turan[11]when Pr=100

圖6給出了不同溫降時刻，1 000方單盤式浮頂油罐內(nèi)含蠟原油流函數(shù)的分布。值得指出的是，流函數(shù)的分布反映了罐內(nèi)含蠟原油的流動情況。從圖中可以看出，罐內(nèi)流函數(shù)演化過程可以分為兩個階段：(1)溫降時間小于200 h。在這一階段，罐內(nèi)核心區(qū)油溫高于含蠟原油的顯觸點(diǎn)，核心區(qū)含蠟原油以純液態(tài)或固液分散體系形式存在，流動性較好。而且隨外界氣溫的劇烈波動，罐內(nèi)流函數(shù)處于瞬態(tài)變化之中。此外，在外界大氣冷量以及罐內(nèi)湍流自然對流作用下，罐內(nèi)含蠟原油的流動在整個油罐區(qū)域內(nèi)進(jìn)行。(2)溫降300～600 h。在這一階段，罐內(nèi)核心區(qū)油溫已經(jīng)降低到含蠟原油顯觸點(diǎn)以下，蠟晶多孔介質(zhì)結(jié)構(gòu)已經(jīng)形成，含蠟原油的流動性明顯減弱，表現(xiàn)為流函數(shù)值僅約為溫降200 h的1%。隨著油溫繼續(xù)降低，析蠟量繼續(xù)增大，蠟晶多孔介質(zhì)孔隙度繼續(xù)減小，流動性繼續(xù)減弱。在外界大氣冷量作用下，罐內(nèi)含蠟原油的流動僅限于油罐右上部。

圖4 正二十烷融化過程中不同時刻融化界面位置與實(shí)驗(yàn)數(shù)據(jù)對比Fig. 4 The comparison between calculated values and experimental results

當(dāng)罐內(nèi)含蠟原油還未完全膠凝時，單盤式浮頂油罐內(nèi)油溫近似均勻，且整體發(fā)生溫降，本文稱之為油罐核心區(qū)。本文通過考察罐內(nèi)x=6.16 m，r=4.18 m處的油溫變化來研究核心區(qū)含蠟原油溫降情況，通過考察罐內(nèi)平均油溫來研究油罐整體溫降情況，計(jì)算結(jié)果如圖7所示。其中，平均油溫＜T＞由式(21)計(jì)算。式中，為罐內(nèi)第i、j控制容積的油溫，ΔVij為第i、j控制容積的體積。

圖7給出了單盤式浮頂油罐核心區(qū)油溫和平均油溫隨溫降時間的變化曲線。從圖中可以看出，核心區(qū)溫降曲線存在兩個轉(zhuǎn)折點(diǎn)：含蠟原油的析蠟點(diǎn)(27 ℃)和顯觸點(diǎn)(23 ℃)。具體分析可知，當(dāng)油溫降低到含蠟原油析蠟點(diǎn)以下時，含蠟原油中開始有蠟晶析出，并釋放析蠟潛熱。而析蠟潛熱的釋放會對罐內(nèi)含蠟原油溫降過程起到緩沖作用。當(dāng)油溫繼續(xù)降低到含蠟原油顯觸點(diǎn)以下時，罐內(nèi)已經(jīng)析出的蠟晶開始相互交聯(lián)，形成蠟晶多孔介質(zhì)結(jié)構(gòu)，造成罐內(nèi)含蠟原油的傳熱方式從自然對流轉(zhuǎn)變?yōu)橄灳Ф嗫捉橘|(zhì)傳熱，類似于導(dǎo)熱，從而使得罐內(nèi)含蠟原油的傳熱過程明顯變慢。此外，當(dāng)油溫降低到含蠟原油顯觸點(diǎn)以下時，析蠟潛熱的傳遞，以及油溫降低后從核心區(qū)獲得的熱量補(bǔ)償均變慢，在二者的共同作用下，罐內(nèi)溫降過程出現(xiàn)了階梯狀下降的特點(diǎn)，這與文獻(xiàn)中水冷凝的實(shí)驗(yàn)結(jié)果[22]類似。

此外，從圖7中還可以看出，當(dāng)油溫高于含蠟原油顯觸點(diǎn)時，罐內(nèi)核心區(qū)油溫與平均油溫近似相等，而當(dāng)油溫降低到顯觸點(diǎn)以下時，核心區(qū)油溫與平均油溫的差異逐漸顯現(xiàn)。這是因?yàn)?，?dāng)罐內(nèi)油溫高于含蠟原油顯觸點(diǎn)時，罐內(nèi)含蠟原油傳熱方式為湍流自然對流，在湍流自然對流過程的攪動下，罐內(nèi)油溫近似均勻，核心區(qū)油溫與平均油溫近似相等。當(dāng)油溫降低到含蠟原油顯觸點(diǎn)以下時，罐內(nèi)含蠟原油已經(jīng)膠凝，從罐頂、罐底和罐壁處傳入的冷量開始冷卻罐內(nèi)含蠟原油，罐內(nèi)出現(xiàn)了明顯的溫度梯度，如圖5(d)～圖5(f)，核心區(qū)油溫將明顯高于平均油溫。

3.2 單盤式浮頂油罐罐底和罐頂凝油層的增長

當(dāng)油溫降低到含蠟原油顯觸點(diǎn)以下時，罐內(nèi)含蠟原油中析出的蠟晶開始相互交聯(lián)，逐漸形成蠟晶多孔介質(zhì)結(jié)構(gòu)。該結(jié)構(gòu)一旦形成就會嚴(yán)重束縛含蠟原油，造成含蠟原油流動性顯著變差。因此，含蠟原油的顯觸點(diǎn)可看成凝油層開始形成的臨界溫度。

圖5 不同溫降時刻，罐內(nèi)含蠟原油溫度場和速度矢量場分布Fig. 5 Oil temperature and velocity vector distributions in the tank at different moments

圖8給出了單盤式浮頂油罐內(nèi)r=5.52 m處罐頂和罐底凝油層隨溫降時間的增長。從圖中可以看出，隨著溫降過程的進(jìn)行，罐頂和罐底凝油層均表現(xiàn)出指數(shù)增長的特點(diǎn)。當(dāng)罐內(nèi)含蠟原油完全膠凝后，罐頂和罐底凝油層厚度達(dá)到最大，均為油罐液高。具體分析可知，初始時刻，罐底部將形成層狀的油溫分布，如圖5(a)，凝油層出現(xiàn)。隨著罐內(nèi)湍流自然對流的形成，罐底凝油層將受到來自核心區(qū)含蠟原油的加熱和沖刷，凝油層緩慢融化，溫降45 h后，罐底凝油層消失。隨著溫降過程的繼續(xù)進(jìn)行，罐內(nèi)油溫逐漸降低，罐底和罐頂凝油層又相繼出現(xiàn)，且快速增長。其中，罐底是罐內(nèi)油溫最低的區(qū)域，且傳熱方式以導(dǎo)熱為主，因此罐底部凝油層出現(xiàn)相對較早，且增長緩慢。而罐頂部湍流自然對流過程比較劇烈，含蠟原油損失的熱量能快速從油罐核心區(qū)得到補(bǔ)充。因此，罐頂部凝油層出現(xiàn)相對較晚，且隨著罐內(nèi)油溫的降低，凝油層將快速增長，如圖8所示。

3.3 單盤式浮頂油罐罐頂、罐底和罐壁處的傳熱特性

3.1節(jié)和3.2節(jié)分別研究了單盤式浮頂油罐內(nèi)含蠟原油的溫降規(guī)律和凝油層增長規(guī)律，本節(jié)主要研究單盤式浮頂油罐內(nèi)含蠟原油與罐頂、罐底和罐壁鋼板層之間的傳熱特性，其中，傳熱量由式(22)計(jì)算。

方程中，i代表計(jì)算位置編號；N代表計(jì)算邊界網(wǎng)格總數(shù)；is代表傳熱面積；代表法向溫度梯度。

圖9給出了單盤式浮頂油罐罐頂、罐底和罐壁處傳熱量隨溫降時間的變化曲線。從圖中可以看出，初始時刻，單盤式浮頂油罐內(nèi)含蠟原油與外界環(huán)境之間的熱交換主要發(fā)生在罐底部，其中，最大傳熱量達(dá)到約5.5 kW。隨著罐底溫度梯度的形成，罐底傳熱量快速減小。當(dāng)溫降時間超過220 h后，罐底部出現(xiàn)了凝油層，傳熱方式由自然對流轉(zhuǎn)變?yōu)橄灳Ф嗫捉橘|(zhì)傳熱，類似于導(dǎo)熱，傳熱量小于0.2 kW。對于罐頂部，當(dāng)溫降時間小于220 h時，罐頂部以湍流自然對流形式向外傳熱，且隨著湍流強(qiáng)度的增大，傳熱量緩慢增大，最大傳熱量超過2.0 kW。隨著溫降過程的進(jìn)行，罐頂部開始出現(xiàn)凝油層，傳熱方式轉(zhuǎn)變?yōu)橄灳Ф嗫捉橘|(zhì)傳熱，傳熱量逐漸減小，最終穩(wěn)定在0.8 kW左右。此外，由于罐壁處包裹有0.1 m厚的保溫層，在整個溫降過程中，罐壁處的傳熱量相對穩(wěn)定，不超過0.3 kW。

圖6 不同溫降時刻，罐內(nèi)含蠟原油流場分布Fig. 6 Stream function distribution in the tank at different moments

圖7 罐內(nèi)核心區(qū)油溫和平均油溫變化曲線Fig. 7 The core region temperature and average temperature in the tank

綜上所述，對于單盤式浮頂油罐，在罐壁保溫層作用下，罐壁處傳熱過程相對穩(wěn)定，傳熱量較小。相對而言，罐頂和罐底傳熱量較大，且變化過程比較復(fù)雜。在短期儲存時，罐頂和罐底部是主要的傳熱區(qū)域；而在長期儲存時，罐頂是最主要的傳熱區(qū)域。因此，在單盤式浮頂油罐保溫設(shè)計(jì)中，要重點(diǎn)關(guān)注罐頂和罐底部。

圖8 罐頂和罐底凝油層的增長Fig. 8 Growth of gel oil at the bottom and top of the tank

4 結(jié)論

基于大渦模擬方法，本文建立了考慮含蠟原油非牛頓性和析蠟相變過程的單盤式浮頂油罐內(nèi)含蠟原油湍流溫降物理數(shù)學(xué)模型，并發(fā)展了一體化耦合求解技術(shù)。采用該模型，本文對實(shí)際1 000方單盤式浮頂油罐內(nèi)含蠟原油溫降過程進(jìn)行計(jì)算，研究罐內(nèi)含蠟原油溫度場和流場的演化、凝油層增長以及傳熱規(guī)律，得到以下結(jié)論：

(1)當(dāng)罐內(nèi)含蠟原油還沒完全膠凝時，在湍流自然對流作用下，罐內(nèi)油溫近似均勻，且整體發(fā)生溫降。當(dāng)罐內(nèi)含蠟原油膠凝后，罐內(nèi)傳熱過程為蠟晶多孔介質(zhì)傳熱，類似于導(dǎo)熱。此時，從罐頂、罐底和罐壁處傳入油罐的冷量開始向油罐核心區(qū)推進(jìn)，油罐中出現(xiàn)了明顯的冷熱油界面。

(2)就凝油層增長而言，對于本文算例，溫降初期，罐底部出現(xiàn)了凝油層。然而，隨著罐內(nèi)湍流過程的形成，凝油層逐漸被加熱和沖刷，最終完全融化。到溫降后期，隨著罐內(nèi)油溫的進(jìn)一步降低，罐頂和罐底再次出現(xiàn)凝油層，且凝油層表現(xiàn)出指數(shù)的、波動式增長的特點(diǎn)。

(3)通過對罐內(nèi)傳熱過程的研究發(fā)現(xiàn)，在罐壁保溫層作用下，罐壁處傳熱過程相對穩(wěn)定，傳熱量較小，而罐頂和罐底傳熱量相對較大，且變化過程更加復(fù)雜。就短期儲存而言，罐頂和罐底是主要的傳熱區(qū)域；而就長期儲存而言，罐頂是最主要的傳熱區(qū)域。因此，在單盤式浮頂油罐保溫設(shè)計(jì)中，要重點(diǎn)關(guān)注罐頂和罐底部。

[1] 郭光臣, 董文蘭, 張志廉. 油庫設(shè)計(jì)與管理[M]. 東營： 中國石油大學(xué)出版社, 2006. [GUO G C, DONG W L, ZHANG Z L. The design and management of tank farm[M]. Dongying： China University of Petroleum Press, 2006.]

[2] 楊筱蘅. 輸油管道設(shè)計(jì)與管理[M]. 東營： 中國石油大學(xué)出版社, 2006. [YANG X H. Oil pipeline design and management[M]. Dongying： China University of Petroleum Press, 2006.]

[3] 于達(dá). 大型浮頂油罐測溫系統(tǒng)的研發(fā)[J]. 油氣儲運(yùn), 2005, 24(8)： 41-43. [YU D. The development of the oil temperature measured system in large fl oating roof oil tank[J]. Oil and Gas Storage and Transportation, 2005, 24(8)： 41-43.]

[4] 趙志明. 大慶北油庫浮頂儲油罐非穩(wěn)態(tài)傳熱問題的數(shù)值計(jì)算[D]. 大慶： 大慶石油學(xué)院, 2009. [ZHAO Z M. The numerical computation of transient heat transfer problems of the fl oating-roof tank in Daqing N. Tank farm[D]. Daqing： Daqing Petroleum Institute, 2009.]

[5] 李超, 劉人瑋, 李旺等. 大型浮頂儲罐原油溫度場實(shí)驗(yàn)測試研究[J]. 工程熱物理學(xué)報, 2013, 34(12)： 2 332-2 334. [LI C, LIU R W, LI W, et al. Crude oil temperature measurement in a large-scale fl oating roof tank[J]. Journal of Engineering Thermophysics, 2013, 34(12)： 2 332-2 334.]

[6] LIN W X, ARMFIELD S W. Long-term behavior of cooling fluid in a vertical cylinder[J]. International Journal of Heat and Mass Transfer, 2005, 48： 53-66.

[7] OLIVESKI R D C. Correlation for the cooling process of vertical storage tanks under natural convection for high Prandtl number[J]. International Journal of Heat and Mass Transfer, 2013, 57(1)： 292-298.

[8] OLIVESKI R D C, MACAGNAN M H, COPETTI J B, et al. Natural convection in a tank of oil： Experimental validation of a numerical code with prescribed boundary condition[J]. Experimental Thermal and Fluid Science, 2005, 29(6)： 671-680.

[9] 李旺. 大型浮頂油罐溫度場數(shù)值模擬方法及規(guī)律研究[D]. 北京： 中國石油大學(xué)(北京), 2013. [LI W. Study on numerical simulation methods and regularities for temperature fi eld of large fl oating roof oil tank[D]. Beijing： China University of Petroleum (Beijing), 2013.]

[10] MATIN M H, POP I, KHANCHEZAR S. Natural convection of power-law fl uid between two-square eccentric duct annuli[J]. Journal of Non-Newtonian Fluid, 2013, 197： 11-23.

[11] TURAN O, SACHDEVA A, CHAKRABORTY N, et al. Laminar natural convection of power-law fl uids in a square enclosure with differentially heated side walls subjected to constant temperatures[J]. Journal of Non-Newtonian Fluid, 2011, 166： 1 049-1 063.

[12] VOLLER V R, PRAKASH C. A fi xed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[J]. International Journal of Heat and Mass Transfer. 1987, 30(8)： 1 709-1 719.

[13] JONES B J, SUN D, KRISHNAN S, et al. Experimental and numerical study of melting in a cylinder[J]. International Journal of Heat and Mass Transfer. 2006, 49(15-16)： 2 724-2 738.

[14] OMARI K E, KOUSKSOU T, GUER Y L. Impact of shape of container on natural convection and melting inside enclosures used for passive cooling of electronic devices[J]. Applied Thermal Engineering. 2011, 31(14-15)： 3 022-3 035.

[15] 張兆順, 崔貴香, 徐春曉. 湍流大渦數(shù)值模擬的理論與應(yīng)用[M]. 北京： 清華大學(xué)出版社, 2008. [ZHANG Z S, CUI G X, XU C X. Theory and application of Large Eddy Simulation[M]. Beijing： Tsinghua University Press, 2008.]

[16] ANSYS FLUENT Theory Guide[M]. ANSYS Inc., 2011.

[17] 袁世偉. 冪律非牛頓流體流動的數(shù)值計(jì)算與實(shí)驗(yàn)研究[D]. 上海： 華東理工大學(xué), 2014. [YUAN S W. Numerical simulation andexperimental study of power-law fl uid[D]. Shanghai： East China University of Science and Technology, 2014.]

[18] WANG M, YU G J, ZHANG X Y, et al. Numerical investigation of melting of the waxy crude oil in an oil tank[J]. Applied Thermal Engineering. 2017, 115： 81-90.

[19] JASAK H, WELLER H G, GOSMAN A D. High resolution NVD differencing scheme for arbitrarily unstructured meshes[J]. International Journal for Numerical Methods in Fluids. 1998, 31(2)： 431-449.

[20] 陶文銓. 計(jì)算傳熱學(xué)的近代進(jìn)展[M]. 北京： 科學(xué)出版社, 2001. [TAO W Q. Recent advances in computational heat transfer[M]. Beijing： Science Press, 2001.]

[21] AMPOFO F, KARAYIANNIS T G. Experimental benchmark data for turbulent natural convection in an air filled square cavity[J]. International Journal of Heat and Mass Transfer. 2003, 46(19)： 3 551-3 572.

[22] KOWALCZYK W, HARTMANN C, DELGADO A. Modeling and numerical simulation of convection driven high pressure induced phase change[J]. International Journal of Heat and Mass Transfer. 2004, 47(5)： 1 079-1 089.

Numerical study on the temperature drop characteristics of waxy crude oil in a single-plate fl oating roof oil tank

WANG Min1, LI Jingfa1, ZHANG Xinyu2, YU Bo3
1 National Engineering Laboratory for Pipeline Safety/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology/ China University of Petroleum-Beijing, Beijing 102249, China 2 Sinopec International Petroleum Exploration and Production Corporation, Beijing, 100029 China 3 School of Mechanical Engineering, Beijing Institute of Petrochemical Technology, Beijing 102617, China

Physical and mathematical models for the temperature drop process of the waxy crude oil in the single-plate fl oating roof oil tank are established with the consideration of non-Newtonian behavior and wax precipitation, and a coupled integrative numerical procedure is developed. In this research, the non-Newtonian behavior is described by the power-law equation. The turbulent fl ow is calculated by Large Eddy Simulation (LES) method. Wax precipitation of the waxy crude oil in the tank is simulated by the enthalpy-porous media theory, in which the gel oil interface is tracked. The numerical program is validated by the results obtained from the literatures. Based on the verif i ed program, temperature drop process in a single-plate fl oating roof oil tank is calculated, and the evolution characteristics of the waxy crude oil in the tank are studied. Furthermore, gel oil growth on the tank top and tank bottom and heat fl ux on tank top, tank bottom and tank wall are also investigated in this research. The research results provide guidance for the scientif i c designing of the single-plate fl oating roof oil tank and making schedule of waxy crude oil turnover reasonably.

single-plate fl oating roof oil tank; temperature drop characteristics; non-newtonian behavior; phase change; LES

10.3969/j.issn.2096-1693.2017.02.025

(編輯 馬桂霞)

王敏, 李敬法, 張欣雨, 宇波. 單盤式浮頂油罐內(nèi)含蠟原油溫降規(guī)律數(shù)值計(jì)算研究. 石油科學(xué)通報, 2017, 02： 267-278

WANG Min, LI Jingfa, ZHANG Xinyu, YU Bo. Numerical study on the temperature drop characteristics of waxy crude oil in a single-plate fl oating roof oil tank. Petroleum Science Bulletin, 2017, 02： 267-278. doi： 10.3969/j.issn.2096-1693.2017.02.025

*通信作者, yubobox@vip.163.com

2017-03-15

國家自然科學(xué)基金(No. 51325603)資助