拱頂罐油氣空間溫度分布及其變化規律研究*

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

1中國石油大學（北京）機械與儲運工程學院·油氣管道輸送安全國家工程實驗室

2玉門油田分公司工程技術研究院

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

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

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

1 數學模型

1.1 模型建立

油罐氣體空間溫度分布主要受到環境溫度、太陽輻射和油面溫度等因素的影響，罐內氣體與環境和油品的熱量交換是一個復雜的氣-液-固三態傳熱過程，傳熱包括熱傳導、對流換熱和輻射三種形式[15-16]。為方便建立數學模型，作以下假設：油罐為密閉空間不存在自然通風；油氣內部只存在導熱形式的傳熱；油氣溫度分布只是高度和半徑的函數；罐頂和罐壁受到的太陽輻射分布均勻。氣體空間溫度變化的熱量來源主要包括經罐頂、罐壁和油面傳入的熱流量，經罐頂、罐壁傳入的熱流量包括大氣傳入的熱流量和油罐接收的輻射熱流量，油面與氣體的換熱是以對流換熱形式完成的。罐內氣體空間換熱結構原理如圖1 所示。

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

可將罐內氣體空間視為圓柱形，忽略向陽或背陽對溫度分布的影響，根據其對稱性建立柱坐標下的二維穩態導熱數學模型[17]。

式中：λ為罐內氣體導熱熱阻，w/（m·k）；qw為對稱中心處傳入的熱流量，W/m2；qe為經罐壁傳入的熱流量，W/m2；Qe為罐壁輻射熱流量，W/m2；Tf為環境溫度，K；Te,b為罐壁處氣體溫度，K；re1、re2和re3分別為大氣到油罐外壁傳熱熱阻、油罐鋼板導熱熱阻和油罐內壁到罐內氣體傳熱熱阻，m2·K/W；qs為經油面傳入的熱流量，W/m2；To為油品表面溫度，K；Ts,b為油品表面處氣體溫度，K；rs1為油品表面到氣體空間的傳熱熱阻，m2·K/W；qn為經罐頂傳入的熱流量，W/m2；Qn為罐頂輻射熱流量，W/m2；Tn,b為罐壁處氣體溫度，K；rn1、rn2和rn3分別為大氣到罐頂外壁傳熱熱阻、罐頂鋼板導熱熱阻和罐頂內壁到罐內氣體傳熱熱阻，m2·k/W。

1.2 方程離散

有限容積法是求解溫度場微分方程廣泛使用的離散化方法之一，其基本思想是將計算域分成若干互不重疊的有限容積單元，使每一個網格結點都由一個有限容積單元所包圍，這個有限容積單元在單位時間內接收或傳出熱量僅與其相鄰的單元有關，并且等于相鄰單元通過界面傳入或傳出熱量的總和[18]。柱坐標下導熱微分方程有限容積法離散如圖2 所示。

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

采用有限容積法對二維穩態導熱微分方程進行離散處理，離散后的形式如下

與某一邊界相鄰的內點，其對應系數值會發生變化，為了便于分析計算，將邊界處內點系數值和邊界條件同時考慮和處理。計算區域的邊界為第三類邊界條件，邊界節點的溫度未知，為了在Tp中不出現未知溫度，需利用己知邊界條件消去。

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

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

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

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

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

2 實例驗證

某地區一座2 000 m3立式拱頂罐，直徑D=14.57 m，罐體總高H=12 m，氣體空間高度Hg=9.7 m。罐內氣體導熱系數λg=0.028 W/（m·K），油面對流換熱系數h=11 W/（m·K），相關傳熱熱阻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°，測得部分時間的環境溫度和油面溫度如表1所示，根據已知參數即可求得油罐氣體空間的溫度分布及其變化規律。

表1 部分環境溫度和油面溫度測量數據Tab.1 Partial measurement data of ambient temperature and oil level temperature

2.1 計算分析

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

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

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

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

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

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

迭代計算得到不同時間氣體空間各個節點的溫度值，利用Python 軟件強大的數據可視化功能對計算數據進行處理，能夠更加直觀地觀測和分析溫度分布情況。圖5 為不同時刻罐內氣體空間平均溫度。由圖5 可知：日出前氣體空間的溫度處于油面溫度和大氣溫度之間；日出后由于太陽輻射熱流量的傳入，氣體空間溫度逐漸增加，很快就超過油面溫度和環境溫度，在午后2 h 左右達到最大值，之后溫度又開始降低；日落后溫度又處于油面溫度和大氣溫度之間，這種狀態一直持續到第二天日出。圖6 為罐內氣體空間不同高度的平均溫度。由圖6可知：日出前大氣溫度低于油面溫度，氣體溫度隨著高度的增加而減小；日出后到日落前大氣溫度高于油面溫度，同時還有太陽輻射的作用，氣體溫度隨著高度的增加而增加；日落時溫度的遞增趨勢變為先增后減，之后這種變化趨勢與日出前一樣。

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

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

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

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

2.2 模型驗證

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

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

為了進一步驗證數學模型的可靠性，進行計算數據與實測數據灰色關聯分析和誤差分析。灰色關聯分析為溫度分布曲線趨勢提供了量化的度量，誤差分析用來量化判定計算數據的精確性，關聯系數計算公式為[19-20]

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

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

表2 計算值與實測值的關聯度和誤差Tab.2 Correlation and error between calculated value and measured value

3 結論

（1）針對油罐氣體空間與環境、油品之間的復雜氣液固三態傳熱問題，綜合考慮環境溫度、油面溫度和太陽輻射等影響因素，建立氣體空間二維穩態傳熱數學模型，采用有限容積法對模型進行離散處理，引入當量溫度概念轉化邊界條件，簡化模型求解過程。

（2）根據某地區一座拱頂罐實測參數，利用Python 軟件編程計算得到不同時間溫度分布數值解，計算結果能夠準確反映罐內氣體空間的溫度分布。為了進一步驗證模型的準確度，用計算溫度和實測溫度進行誤差分析和灰色關聯分析，驗證得到模型具有較高的準確性。

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

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