肖開喜 侯磊 趙智 李延豪 黃亞楠 張蕊

1中國石油大學(xué)（北京）機(jī)械與儲運(yùn)工程學(xué)院·油氣管道輸送安全國家工程實(shí)驗(yàn)室

2玉門油田分公司工程技術(shù)研究院

石油及其產(chǎn)品在儲存和收發(fā)作業(yè)過程中，由于受到工藝技術(shù)和設(shè)備的限制，不可避免地會有部分輕質(zhì)液態(tài)烴蒸發(fā)到罐內(nèi)氣體空間，通過呼吸閥逸入大氣，造成無法預(yù)估的油品損失和環(huán)境污染[1-3]。早在20 世紀(jì)40 年代，國外相關(guān)機(jī)構(gòu)和學(xué)者就開始對油品蒸發(fā)損耗理論進(jìn)行研究，試圖從理論上闡述油品蒸發(fā)的機(jī)理，探明影響蒸發(fā)損耗的因素和相關(guān)參數(shù)的變化規(guī)律[4-6]。與油品蒸發(fā)性質(zhì)最為相關(guān)的參數(shù)是溫度，對罐內(nèi)溫度分布及其變化規(guī)律的研究是油品蒸發(fā)損耗研究工作的重要組成部分。

楊宏偉等[7-10]實(shí)地測量了罐外環(huán)境溫度、罐內(nèi)油品溫度和油罐氣體空間不同位置的溫度，對實(shí)測數(shù)據(jù)進(jìn)行分析和處理，得到了氣體空間溫度分布規(guī)律。梁穎[11]在對某固定頂油罐溫度分布實(shí)測數(shù)據(jù)進(jìn)行分析的基礎(chǔ)上，利用ANSYS 軟件對溫度場進(jìn)行模擬，根據(jù)模擬數(shù)據(jù)擬合出氣體溫度分布公式，確定“小呼吸”起始與終了時(shí)的溫度。VISSER 等[12]對液化氣儲罐內(nèi)非穩(wěn)態(tài)傳熱進(jìn)行模擬，得到儲罐內(nèi)溫度分布和到達(dá)液化氣的熱流密度，用溫度分布和熱流密度來預(yù)測儲罐的蒸發(fā)損失，通過最小化蒸發(fā)損失優(yōu)化儲氣罐的傳熱設(shè)計(jì)。張麗娜等[13]為研究呼吸閥凍結(jié)問題，建立了油罐氣體空間溫度分布二維穩(wěn)態(tài)導(dǎo)熱模型，采用有限差分方法將其差分格式化，運(yùn)用列主元消去法求解方程，得到了氣體空間溫度分布規(guī)律。馬曉宇等[14]根據(jù)實(shí)測數(shù)據(jù)建立了某固定頂油罐氣體空間傳熱模型，對計(jì)算溫度值進(jìn)行分段線性擬合分析，為深入研究蒸發(fā)損耗提供了理論依據(jù)。

對油罐氣體空間溫度分布及其變化規(guī)律的研究主要有現(xiàn)場溫度實(shí)測、專業(yè)軟件模擬和數(shù)學(xué)模型計(jì)算三種方式，目前研究主要針對某一時(shí)刻的溫度分布，對其變化規(guī)律的研究不夠深入，對溫度分布的影響因素處理較為簡單。本文考慮地理環(huán)境、大氣溫度、太陽輻射和油品溫度等影響因素，建立油罐氣體空間二維傳熱數(shù)學(xué)模型；采用有限容積法對數(shù)學(xué)模型和邊界條件進(jìn)行離散處理，引入當(dāng)量溫度概念轉(zhuǎn)化邊界條件，簡化模型求解過程；結(jié)合實(shí)際案例編寫計(jì)算程序，對計(jì)算結(jié)果進(jìn)行關(guān)聯(lián)分析和誤差分析，驗(yàn)證了模型具有很高的準(zhǔn)確度。

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

1.1 模型建立

油罐氣體空間溫度分布主要受到環(huán)境溫度、太陽輻射和油面溫度等因素的影響，罐內(nèi)氣體與環(huán)境和油品的熱量交換是一個(gè)復(fù)雜的氣-液-固三態(tài)傳熱過程，傳熱包括熱傳導(dǎo)、對流換熱和輻射三種形式[15-16]。為方便建立數(shù)學(xué)模型，作以下假設(shè)：油罐為密閉空間不存在自然通風(fēng)；油氣內(nèi)部只存在導(dǎo)熱形式的傳熱；油氣溫度分布只是高度和半徑的函數(shù)；罐頂和罐壁受到的太陽輻射分布均勻。氣體空間溫度變化的熱量來源主要包括經(jīng)罐頂、罐壁和油面?zhèn)魅氲臒崃髁?，?jīng)罐頂、罐壁傳入的熱流量包括大氣傳入的熱流量和油罐接收的輻射熱流量，油面與氣體的換熱是以對流換熱形式完成的。罐內(nèi)氣體空間換熱結(jié)構(gòu)原理如圖1 所示。

圖1 油罐熱流圖Fig.1 Heat flow diagram of oil tank

可將罐內(nèi)氣體空間視為圓柱形，忽略向陽或背陽對溫度分布的影響，根據(jù)其對稱性建立柱坐標(biāo)下的二維穩(wěn)態(tài)導(dǎo)熱數(shù)學(xué)模型[17]。

式中：λ為罐內(nèi)氣體導(dǎo)熱熱阻，w/（m·k）；qw為對稱中心處傳入的熱流量，W/m2；qe為經(jīng)罐壁傳入的熱流量，W/m2；Qe為罐壁輻射熱流量，W/m2；Tf為環(huán)境溫度，K；Te,b為罐壁處氣體溫度，K；re1、re2和re3分別為大氣到油罐外壁傳熱熱阻、油罐鋼板導(dǎo)熱熱阻和油罐內(nèi)壁到罐內(nèi)氣體傳熱熱阻，m2·K/W；qs為經(jīng)油面?zhèn)魅氲臒崃髁?，W/m2；To為油品表面溫度，K；Ts,b為油品表面處氣體溫度，K；rs1為油品表面到氣體空間的傳熱熱阻，m2·K/W；qn為經(jīng)罐頂傳入的熱流量，W/m2；Qn為罐頂輻射熱流量，W/m2；Tn,b為罐壁處氣體溫度，K；rn1、rn2和rn3分別為大氣到罐頂外壁傳熱熱阻、罐頂鋼板導(dǎo)熱熱阻和罐頂內(nèi)壁到罐內(nèi)氣體傳熱熱阻，m2·k/W。

1.2 方程離散

有限容積法是求解溫度場微分方程廣泛使用的離散化方法之一，其基本思想是將計(jì)算域分成若干互不重疊的有限容積單元，使每一個(gè)網(wǎng)格結(jié)點(diǎn)都由一個(gè)有限容積單元所包圍，這個(gè)有限容積單元在單位時(shí)間內(nèi)接收或傳出熱量僅與其相鄰的單元有關(guān)，并且等于相鄰單元通過界面?zhèn)魅牖騻鞒鰺崃康目偤蚚18]。柱坐標(biāo)下導(dǎo)熱微分方程有限容積法離散如圖2 所示。

圖2 柱坐標(biāo)下有限容積示意圖Fig.2 Schematic diagram of finite volume in cylindrical coordinates

采用有限容積法對二維穩(wěn)態(tài)導(dǎo)熱微分方程進(jìn)行離散處理，離散后的形式如下

與某一邊界相鄰的內(nèi)點(diǎn)，其對應(yīng)系數(shù)值會發(fā)生變化，為了便于分析計(jì)算，將邊界處內(nèi)點(diǎn)系數(shù)值和邊界條件同時(shí)考慮和處理。計(jì)算區(qū)域的邊界為第三類邊界條件，邊界節(jié)點(diǎn)的溫度未知，為了在Tp中不出現(xiàn)未知溫度，需利用己知邊界條件消去。

針對右邊界處，式（2）可表示為

其中：ae(Te-Tp)=qerpΔz，qe為流進(jìn)控制體的熱流量。

式中：re4為右邊界控制體邊緣到中心的導(dǎo)熱熱阻，m2·K/W；為同時(shí)考慮大氣傳入熱流量和罐壁輻射熱流量的當(dāng)量溫度，K。

通過引入當(dāng)量溫度概念將復(fù)雜的第三類邊界條件轉(zhuǎn)換成簡單的第一類邊界條件來求解微分方程。對左邊界、下邊界和上邊界采取相同的處理方式，結(jié)果分別為

式中：rs2為下邊界控制體邊緣到中心的導(dǎo)熱熱阻，m2·K/W；rn4為上邊界控制體邊緣到中心的導(dǎo)熱熱阻，m2·K/W；為同時(shí)考慮大氣傳入熱流量和罐頂輻射熱流量的當(dāng)量溫度，K。

2 實(shí)例驗(yàn)證

某地區(qū)一座2 000 m3立式拱頂罐，直徑D=14.57 m，罐體總高H=12 m，氣體空間高度Hg=9.7 m。罐內(nèi)氣體導(dǎo)熱系數(shù)λg=0.028 W/（m·K），油面對流換熱系數(shù)h=11 W/（m·K），相關(guān)傳熱熱阻re1=0.08 m2·K/W、re2=0.002 m2·K/W、re3=0.08 m2·K/W，rn1=0.1 m2·K/W、rn2=0.002 m2·K/W、rn3=0.1 m2·K/W。

以初夏某一天為例，查得正午太陽天頂角θ0=35°，測得部分時(shí)間的環(huán)境溫度和油面溫度如表1所示，根據(jù)已知參數(shù)即可求得油罐氣體空間的溫度分布及其變化規(guī)律。

表1 部分環(huán)境溫度和油面溫度測量數(shù)據(jù)Tab.1 Partial measurement data of ambient temperature and oil level temperature

2.1 計(jì)算分析

罐頂或罐壁在不同時(shí)間、不同位置接收到的太陽輻射熱流量不同，考慮大氣對太陽輻射的削弱作用和罐頂、罐壁對輻射的吸收能力，與陽光照射方向垂直的單位面積罐頂、罐壁所吸收的輻射熱流量用式（8）計(jì)算。

以正午太陽輻射強(qiáng)度計(jì)算，得到q0=420.05 W/m2。一天里陽光不能總是垂直照射在罐頂、罐壁上，根據(jù)蘭貝特定律，任一時(shí)刻單位面積罐頂、罐壁所吸收的太陽輻射熱流量可用式（9）計(jì)算。

式中：Qn、Qe分別為罐頂和罐壁單位面積吸收的輻射熱流量，W/m2；圓頻率ω=π/12；日出時(shí)間τrc=6。

計(jì)算得到一天內(nèi)不同時(shí)間單位面積罐頂、罐壁所吸收的輻射熱流量結(jié)果如圖3 所示，根據(jù)環(huán)境溫度和油品溫度的一般變化規(guī)律，擬合得到一天內(nèi)不同時(shí)間的溫度。確定相關(guān)參數(shù)后，迭代計(jì)算罐內(nèi)氣體空間溫度分布，分析其變化規(guī)律。Python 編程計(jì)算流程如圖4 所示。

圖3 罐頂、罐壁輻射熱流量Fig.3 Radiant heat flow of tank top and wall

圖4 溫度分布計(jì)算流程Fig.4 Calculation process of temperature distribution

迭代計(jì)算得到不同時(shí)間氣體空間各個(gè)節(jié)點(diǎn)的溫度值，利用Python 軟件強(qiáng)大的數(shù)據(jù)可視化功能對計(jì)算數(shù)據(jù)進(jìn)行處理，能夠更加直觀地觀測和分析溫度分布情況。圖5 為不同時(shí)刻罐內(nèi)氣體空間平均溫度。由圖5 可知：日出前氣體空間的溫度處于油面溫度和大氣溫度之間；日出后由于太陽輻射熱流量的傳入，氣體空間溫度逐漸增加，很快就超過油面溫度和環(huán)境溫度，在午后2 h 左右達(dá)到最大值，之后溫度又開始降低；日落后溫度又處于油面溫度和大氣溫度之間，這種狀態(tài)一直持續(xù)到第二天日出。圖6 為罐內(nèi)氣體空間不同高度的平均溫度。由圖6可知：日出前大氣溫度低于油面溫度，氣體溫度隨著高度的增加而減?。蝗粘龊蟮饺章淝按髿鉁囟雀哂谟兔鏈囟?，同時(shí)還有太陽輻射的作用，氣體溫度隨著高度的增加而增加；日落時(shí)溫度的遞增趨勢變?yōu)橄仍龊鬁p，之后這種變化趨勢與日出前一樣。

圖5 氣體空間、環(huán)境和油面平均溫度Fig.5 Average temperature of gas space，environment and oil level

圖6 氣體空間不同高度平均溫度Fig.6 Average temperature of gas space at different height

圖7 為不同時(shí)刻油罐氣體空間溫度云圖。日出前沒有太陽輻射，氣體空間溫度分布受大氣溫度和油面溫度影響，除油面附近一薄層外，同一高度氣體空間溫度基本相等，徑向溫差一般不超過1～2 ℃。日出時(shí)，只有罐壁接收到太陽輻射，罐壁附近不僅有大氣傳入的熱流量，還有太陽輻射熱流量，氣體空間高溫主要在罐壁附近，此時(shí)平均徑向溫差達(dá)到2.15 ℃。日出后，罐頂接收到的太陽輻射逐漸增大，罐壁接收到太陽輻射逐漸減小，在12:00 分別達(dá)到最大值和最小值，氣體空間溫度在縱向上出現(xiàn)明顯的分層，正午過后罐壁接收到的太陽輻射逐漸增加，徑向上又開始出現(xiàn)明顯的溫差。日落時(shí)，又只有罐壁接收到太陽輻射，此時(shí)氣體空間出現(xiàn)明顯的溫度差，平均徑向溫差達(dá)到3.02 ℃，在氣體空間中間位置的徑向溫差尤為明顯。日落后，罐頂和罐壁都沒有接收到太陽輻射，氣體空間的溫度只受油面溫度和大氣溫度影響，罐內(nèi)溫度分布與日出前相同。

圖7 不同時(shí)刻氣體空間溫度分布Fig.7 Gas space temperature distribution at different time

2.2 模型驗(yàn)證

通過對數(shù)值計(jì)算結(jié)果進(jìn)行分析表明，建立的數(shù)學(xué)模型以及處理方式能夠較準(zhǔn)確地反映出油罐內(nèi)部氣體空間溫度分布及其變化規(guī)律。圖8 為氣體空間不同高度的計(jì)算溫度和測量溫度對比圖。由圖8 可知，計(jì)算結(jié)果變化趨勢與實(shí)測值一致，誤差也較小。

圖8 不同高度的平均溫度Fig.8 Average temperature at different height

為了進(jìn)一步驗(yàn)證數(shù)學(xué)模型的可靠性，進(jìn)行計(jì)算數(shù)據(jù)與實(shí)測數(shù)據(jù)灰色關(guān)聯(lián)分析和誤差分析。灰色關(guān)聯(lián)分析為溫度分布曲線趨勢提供了量化的度量，誤差分析用來量化判定計(jì)算數(shù)據(jù)的精確性，關(guān)聯(lián)系數(shù)計(jì)算公式為[19-20]

式中：x0(k)、xi(k)分別為不同高度位置對應(yīng)的測量溫度值和計(jì)算溫度值；α為分辨系數(shù)，一般取0.5。

數(shù)值計(jì)算溫度與實(shí)測溫度的平均絕對關(guān)聯(lián)度和平均絕對誤差值如表2 所示。由表2 可知，所選取7 個(gè)時(shí)間段計(jì)算出的平均關(guān)聯(lián)度均大于0.6，平均溫度誤差均小于2 ℃，計(jì)算溫度與實(shí)測溫度具有良好的關(guān)聯(lián)性和較小的誤差，能夠充分反映油罐氣體空間的溫度分布情況。

表2 計(jì)算值與實(shí)測值的關(guān)聯(lián)度和誤差Tab.2 Correlation and error between calculated value and measured value

3 結(jié)論

（1）針對油罐氣體空間與環(huán)境、油品之間的復(fù)雜氣液固三態(tài)傳熱問題，綜合考慮環(huán)境溫度、油面溫度和太陽輻射等影響因素，建立氣體空間二維穩(wěn)態(tài)傳熱數(shù)學(xué)模型，采用有限容積法對模型進(jìn)行離散處理，引入當(dāng)量溫度概念轉(zhuǎn)化邊界條件，簡化模型求解過程。

（2）根據(jù)某地區(qū)一座拱頂罐實(shí)測參數(shù)，利用Python 軟件編程計(jì)算得到不同時(shí)間溫度分布數(shù)值解，計(jì)算結(jié)果能夠準(zhǔn)確反映罐內(nèi)氣體空間的溫度分布。為了進(jìn)一步驗(yàn)證模型的準(zhǔn)確度，用計(jì)算溫度和實(shí)測溫度進(jìn)行誤差分析和灰色關(guān)聯(lián)分析，驗(yàn)證得到模型具有較高的準(zhǔn)確性。

（3）從日出到日落，油罐氣體空間為蓄熱過程，熱量經(jīng)罐頂和罐壁傳至氣體空間，再傳到油品使油品溫度升高；從日落到第二天日出，氣體空間為放熱過程，熱量經(jīng)罐頂和罐壁傳遞給大氣，同時(shí)油品溫度也下降；罐內(nèi)氣體空間溫度分布受太陽輻射影響最為明顯，油氣空間晝夜溫差可比環(huán)境溫差大10～20 ℃。

（4）溫度分布直接影響到罐內(nèi)壓力和油氣濃度，罐內(nèi)壓力決定呼吸開始和終了時(shí)間，油氣濃度決定呼出氣體中油氣的含量。根據(jù)溫度分布能夠研究罐內(nèi)壓力變化和濃度分布，準(zhǔn)確計(jì)算不同工況下油品蒸發(fā)損耗量，可為采取降低蒸發(fā)損耗措施提供指導(dǎo)依據(jù)。