儲罐中心區填充多孔材料對LNG流動特征影響的數值模擬

邢志祥，張淑淑，汪李金，張 瑩

(常州大學環境與安全工程學院，江蘇 常州 213164)

多孔介質是一種多相物質共存、由固相和氣相或液相組成的物質，固相作為固體骨架，氣相或液相充滿了孔隙空間，孔隙之間必須相互連通[1]。Nield等[2]最先開始了多孔介質的研究，認為多孔材料具有比表面積大、質量輕、體積小、導熱性能好的特點。近年來，針對多孔介質的試驗和模擬研究很多，研究主要集中在填充多孔材料后液體流動的情況以及多孔材料內部火焰傳播的情況，而研究對象通常是管道、密閉容器等，但關于儲罐內填充多孔材料的研究則相對較少。邢志祥等[3]模擬了圓柱形儲罐內填充聚氨酯多孔材料后可燃氣體火焰傳播的情況，并與13 L圓柱形密閉儲罐現場試驗的結果進行了對比，結果表明：模擬結果與試驗結果相吻合，填充物平均孔徑越小，阻火抑爆性能越好；張健中等[4]分析了在加油站埋地油罐中是否有必要使用網狀鋁合金阻隔防爆材料等問題;田宏等[5-6]對填充在液化石油氣儲罐中的多孔金屬材料做了系統的介紹，但并未通過試驗來舉證。

液化天然氣(Liquified Natural Gas,LNG)儲罐是儲存液化天然氣的重要設備，其運行處于超低溫狀態時容易引發各類安全事故，而LNG翻滾是其中一種事故類型。儲罐內由于LNG密度差產生分層，上層液體與下層液體之間存在液液分界面，而漏熱導致分界面被破壞進一步引發翻滾現象。分界面被破壞的原因有兩種[7]：一種是儲罐內壁的邊界層穿透，導致分界面被破壞；另一種是上層液體在中心處向下的射流和下層液體在中心處向下的卷攜，導致中心處的分界面被破壞，進一步造成了整個分界面被破壞。考慮到多孔材料的優點，本文將其應用到儲罐中研究其是否能起到抑制LNG翻滾的作用，并選取儲罐中心區填充多孔材料，對儲罐內LNG的流動特征進行了數值模擬。

1 填充多孔材料的儲罐內液體流動計算模型的建立

本文選取儲罐中心區填充多孔材料，并分別設置多孔材料的填充厚度為1 m、1.5 m、2 m，通過建立儲罐中心區填充多孔材料的物理模型和數學模型，對儲罐內LNG的流動情況進行了數值模擬計算。

1. 1 物理模型

本文選取6 000 m3的圓柱形儲罐，內罐直徑為24 m，高度為17.4 m，設計液位高15.76 m，并取儲罐的軸截面建立多孔介質中湍流流動的二維模型。坐標原點位于左側壁面和底部壁面的交匯處，X軸向右為正，Y軸向上為正；模擬過程不考慮氣相空間，只針對儲罐內的液相部分來模擬LNG的流動情況，液體充裝高度為14 m；多孔材料填充在儲罐的中心部位；對儲罐結構進行了簡化，不考慮壁厚，并忽略儲罐內部的管線、循環泵等部件，由此建立的儲罐中心區填充多孔材料簡化后的物理模型如圖1所示。

圖1 儲罐中心區填充多孔材料簡化后的物理模型Fig.1 Simplified physical model of the filling porous material of the storage tank central area

對儲罐中心區填充多孔材料時儲罐內LNG流動情況的物理模型(中間標明H的區域為多孔材料區域，左右兩側區域為純流體區域)進行了如下的簡化：①儲罐內部流體分成兩個高度相等的分層，分層高度為7 m，上分層流體密度為424 kg/m3，下分層流體密度為425 kg/m3，初始密度差為1 kg/m3，兩個分層內的密度均勻一致；②儲罐內流體密度符合Boussinesq假設；③選取金屬鋁多孔材料，且視多孔材料為均勻各向同性；④多孔材料的孔隙為球形空洞，孔隙之間相互連通，即為開孔結構；⑤多孔材料的物性參數為常數；⑥多孔材料固體骨架與流體處于局部熱平衡，且無任何化學反應；⑦儲罐內流體流動為湍流流動；⑧忽略流體流動過程中的黏性耗散；⑨不考慮輻射換熱。

1. 2 數學模型

采用非穩態模型、VOF模型和標準k-ε湍流模型，并基于Boussinesq假設，計算過程中多孔區域動量方程采用Darcy-Brinkman-Forchheimer模型，能量方程采用局部熱平衡模型，得到的控制方程如下：

連續性方程:

?u?x+?v?y=0

(1)

動量方程:

ρφ2?u?t+u?u?x+v?u?y=-?p?x+

μφ?2u?x2+?2u?y2-μKu+ρCFKu2+v2u

(2)

ρφ2?v?t+u?v?x+v?v?y=-?p?y+

μφ?2v?x2+?2v?y2-μKv+ρCFKu2+v2v-ρg

(3)

能量方程:

ρCp?T?t+u?T?x+v?T?y=λeff?2T?x2+?2T?y2

(4)

k方程：

??t(ρk)+??xi(ρkui)=??xjμ+μtσk?k?xj〗+

Gk+Gb-ρε-YM+Sk

(5)

ε方程:

??t(ρε)+??xi(ρεui)=??xjμ+μtσε?ε?xj〗+

C1εεk(Gk+C3εGb)-C2ερε2k+Sε

(6)

上式中：u為X軸方向的速度分量(m/s)；v為Y軸方向的速度分量(m/s)；ρ為液體密度(kg/m3)；t為時間(s)；p為壓力(Pa)；μ為動力黏度(Pa·s)；g為重力加速度(9.81 m2/s)；φ為孔隙率；K為滲透率(m2)；CF為Forchheimer系數；Cp為定壓比熱容[J/(kg·K)]；T為溫度(K)；λeff為多孔區域的有效導熱系數[W/(m·K)],由流體的導熱系數λf和固體的導熱系數λs的體積平均值計算得到，即λeff=(1-φ)λs+φλf；k為湍動能(J)；ε為湍流擴散率；σk和σε分別為k方程和ε方程的普朗特數，σk=1.0，σε=1.3；μt為湍流黏性系數，μt=Cμρk2/ε，其中Cμ=0.09；Gk為平均速度梯度引起的湍動能k的產生項；Gb為浮力引起的湍動能k的產生項；YM為可壓縮湍流流動中脈動擴張的貢獻；C1ε、C2ε、C3ε為經驗常數，C1ε=1.44，C2ε=1.92，C3ε=1；Sk和Sε為用戶自定義源項。

當φ=1時，k→∞，λeff=λf，動量方程的源項趨近于0，能量方程也變為標準能量方程，則上述控制方程變為純流體區域的控制方程。

1. 3 邊界條件和初始條件

儲罐側壁面和底壁面取為無滑移壁面邊界條件，并采用定熱流的方式加熱，熱流密度恒定為30 W/m2；氣液交界面取為壓力出口；多孔區域采用金屬鋁多孔材料，孔隙率設定為0.95，并需設定黏性阻力系數和慣性阻力系數；多孔區域與純流體區域的交界面設為內部邊界。

初始時刻：初始表壓為15 000 Pa，X和Y軸方向的速度分量都為0，分層區初始溫度為111 K，主流區初始溫度為111.5 K。

在計算過程中，多孔區域和純流體區域通過φ的取值進行區分，并統一進行求解。

1. 4 相關參數介紹

(1) 孔隙率φ和孔徑dp[8]：孔隙率φ是指孔隙體積占多孔材料總體積的比值；孔密度是指多孔材料每英寸上孔的數目(1 inch=0.025 4 m)，孔徑dp=0.025 4/孔密度(m)。

(2) 滲透率K和Forchheimer系數CF[9]：滲透率K代表了多孔介質中孔隙的表面積和彎曲程度，通常由Vafai總結的經驗公式確定，即

當多孔材料中的流體流動為湍流流動時，達西定律便不能描述流體流動的情況，此時可采用其修正模型進行描述，即在動量方程中加入源項，源項由兩項組成，一項為黏性阻力項，一項為慣性阻力項，慣性阻力項的系數用CF表示，稱為Forchheimer系數，可表示為

CF=1.75150φ3/2

(3) 壓差Δp：壓差指儲罐內上下分層之間的壓力差(Pa)，Δp=p2-p1，其中p2為底壁面平均壓力(Pa)，p1為壓力出口平均壓力(Pa)。

2 儲罐內LNG流動特征的數值模擬與分析

本文選取儲罐的軸截面建立二維模型，應用Fluent 6.3軟件對儲罐中心區無多孔材料填充和有多孔材料填充時儲罐內LNG的流動進行了數值模擬，模擬斷面均取自儲罐的同一軸截面，并選取速度云圖和流線圖的模擬結果進行了分析。速度云圖是描述流體速度大小分布情況的圖；流線圖是描述某一時刻流體運動趨勢的圖，可以通過流線圖中的流線方向和流線疏密，判斷出某一時刻流體的運動方向和區分流動強弱區域，流線圖中箭頭表示流體的運動方向，流線的疏密反映流速大小，流線越密，流體流速越大，反之亦然。

2.1 儲罐中心區無多孔材料填充時儲罐內LNG的流動情況

圖2為儲罐中心區無多孔材料填充時,t=100 s和t=1 000 s時刻儲罐內LNG的流動情況模擬結果。

圖2 儲罐中心區無多孔材料填充時t=100 s和t=1 000 s 時刻的儲罐內流體的速度云圖和流線圖Fig.2 Velocity cloud and stream traces of fluid in the storage tank with no filling porous materials at t=100 s and t=1 000 s

由圖2可以看出：

(1) 在t=100 s時，由于外界環境的漏熱，儲罐壁面處的流體率先吸收熱量導致密度減小，產生浮力驅動流，沿著壁面向上移動，形成流動邊界層；由于上下兩層流體存在溫度差和密度差，因此位于分界面處的流體出現自然對流現象，且流體流速較儲罐其他位置的流速大，而分界面處邊界層的流體流速最大[見圖2(a)]。

(2) 在t=1 000 s時，分界面處流體的自然對流發展到整個儲罐內流體的自然對流，且上下兩層流體的自然對流各自獨立，即由液液分界面隔開，被限制在自身區域內。對于上層流體，由于儲罐側壁面的漏熱，導致邊界層處的流體受熱膨脹密度減小，沿儲罐壁面向上移動，當其到達氣液交界面時，其中一部分流體蒸發變成氣體釋放到氣相空間，其密度增大導致向下移動，由此形成了上層流體的自然能對流現象；對于下層流體，邊界層處的流體受熱膨脹密度減小，沿著邊界層向上移動，當其到達液液分界面時，因浮力太小而無法穿過分界面進入上層，但會通過分界面向上層流體傳遞熱量，因而溫度降低，密度增大，導致向下移動，由此形成了下層流體的自然對流現象[見圖2(b)]。

(3) 在純自然對流運動中，流體的瑞利數Ra是判定由浮力產生的對流強度大小的標準[10]。當Ra小于臨界值時，流體之間是熱傳導狀態，不發生對流運動；當Ra大于臨界值時，才會發生對流運動，從而在液體中出現宏觀對流花紋，稱其為Benard花紋[11]。儲罐內無多孔材料填充時，在t=1 000 s時，上層和下層的自然對流結構明顯不同，但流動都是由類似圓形的滾動圈組成，上層的滾動圈數為1個，為扁長形，而下層的滾動圈數為4個，分別為2個胖圓形和2個扁長形的滾動圈；相比于t=100 s時分界面處的流體流速較高，t=1 000 s時分界面處的流體流速反而較小，除此之外，在滾動圈中心處的流速也較小，稱之為滯留區[見圖2(c)]。

在模擬計算過程中，對二維截面上流體的最大流速進行了監測，得到流體的最大流速為0.195 m/s。

2.2 儲罐中心區多孔材料填充厚度為1 m時儲罐內LNG的流動情況

圖3為儲罐中心區多孔材料填充厚度為1 m時，t=100 s和t=1 000 s時刻儲罐內LNG的流動情況模擬結果。

圖3 儲罐中心區多孔材料填充厚度為1 m時t=100 s和t=1 000 s時刻的儲罐內流體的速度云圖和流線圖Fig.3 Velocity cloud and stream traces of fluid in the storage tank with the filling thickness of porous materials being 1 m at t=100 s and t=1 000 s

由圖3可以看出：

(1) 在t=100 s時，儲罐內流體流動發展情況相對于無多孔材料填充時比較滯后，分界面處的流體流速較無多孔材料填充時小，除邊界層外，多孔材料填充區域兩側自然對流強度較高[見圖3(a)]。

(2) 在t=1 000 s時，多孔材料填充區域相比相鄰區域的流體流速小[見圖3(b)]；相對于無多孔材料填充時，上層和下層的自然對流結構有明顯不同，上層的自然對流結構由1個扁長形滾動圈變為幾個眼淚形狀的小滾動圈，且由于多孔材料增加了流體流動阻力，在多孔材料填充區域左側流體流線豎直向下，穿過多孔材料后，便向氣液交界面移動,而下層的自然對流結構由2個胖圓形和2個扁長形的滾動圈變為3個胖圓形的滾動圈，且中間的滾動圈比左右兩側的滾動圈大，滾動圈方向從左到右依次為逆時針、順時針、逆時針，同時發現在壁面熱邊界條件和多孔材料作用下，罐體兩邊的流動強度較強，中間的流動強度較弱，中間滾動圈的滯留區面積比兩側滾動圈大[見圖3(c)]。

儲罐內沿中心區填充多孔材料厚度為1 m時，二維截面上流體的最大流速為0.190 m/s。

2.3 儲罐中心區多孔材料填充厚度為1.5 m時儲罐內LNG的流動情況

圖4為儲罐中心區多孔材料填充厚度為1.5 m時，t=100 s和t=1 000 s時刻儲罐內LNG的流動情況模擬結果。

圖4 儲罐中心區多孔材料填充厚度為1.5 m時t=100 s和 t=1 000 s時刻的儲罐內流體的速度云圖和流線圖Fig.4 Velocity cloud and stream traces of fluid in the storage tank with the filling thickness of porous materials being 1.5 m at t=100 s and t=1 000 s

由圖4可以看出：

(1) 在t=100 s時，仍是邊界層處流體和分界面處流體的流速較儲罐其他位置的流速大，分界面處出現自然對流現象，但與同時期無多孔材料填充時相比，分界面上層流體的流速較高，在多孔材料填充區域兩側各出現了3個“火焰”形的高流速滾動圈，其速度分布同“火焰”的溫度分布一樣，中心速度最大，向四周減弱[見圖4(a)]。

(2) 在t=1 000 s時，多孔材料填充區域的流體流速最小，且下層的自然對流結構明顯沿儲罐中心線呈對稱分布[見圖4(b)]；上層和下層的自然對流結構與無多孔材料填充時有明顯的不同，上層的自然對流結構與多孔材料填充厚度為1 m時一樣，仍為幾個大小不一、眼淚形狀的小滾動圈，而且在多孔材料填充區域內部有一個滾動圈，而下層的自然對流結構沿儲罐中心線呈對稱分布，由2個近似相等的胖圓形和2個扁長形的滾動圈組成，且中間的滾動圈比左右兩側的滾動圈大，中心線兩側的滾動圈方向依次為順時針、逆時針、順時針、逆時針，順時針與逆時針交替出現，同時發現與多孔材料填充厚度為1 m時相似，中間滾動圈的滯留區面積比兩側滾動圈大，兩側滾動圈的滯留區面積幾乎為0[見圖4(c)]。

儲罐內沿中心區填充多孔材料厚度為1.5 m時，二維截面上流體的最大流速為0.371 m/s。

2.4 儲罐中心區多孔材料填充厚度為2 m時儲罐內LNG的流動情況

圖5為儲罐中心區多孔材料填充厚度為2 m時，t=100 s和t=1 000 s時刻儲罐內LNG的流動情況模擬結果。

圖5 儲罐中心區多孔材料填充厚度為2 m時t=100 s和t=1 000 s時刻的儲罐內流體的速度云圖和流線圖Fig.5 Velocity cloud and stream traces of fluid in the storage tank with the filling thickness of porous materials being 2 m at t=100 s and t=1 000 s

由圖5可以看出：

(1) 在t=100 s時，儲罐內流體流動發展情況與同時期多孔材料填充厚度為1 m時較相似，分界面處的流體流速較無多孔材料填充時小，除邊界層外，多孔材料填充區域兩側自然對流強度較高[見圖5(a)]。

(2) 在t=1 000 s時，相對于無多孔材料填充時，低速區域范圍變大，高速區域范圍變小，且速度場與填充厚度為1.5 m時的速度場較相似[見圖5(b)]；上層的自然對流結構為2個小的滾動圈，下層的自然對流結構主要為2個大的滾動圈，滾動圈方向依次為順時針、逆時針，此時的流場結構不再沿儲罐中心線呈對稱分布，而是以中心線為分界線，左右兩側流場具有相似性，說明多孔材料填充厚度為2 m時，極大地增加了流動阻力，好像一堵墻，起到了分割流場的作用[見圖5(c)]。

儲罐內沿中心區填充多孔材料厚度為2 m時，二維截面上流體的最大流速為0.189 m/s。

2.5 儲罐中心區多孔材料填充厚度對儲罐內LNG流動情況的影響分析

圖6為多孔材料不同填充厚度(H=0、1 m、1.5 m、2 m)時儲罐內二維截面上流體平均流速的變化曲線。

圖6 多孔材料不同填充厚度時儲罐內二維截面上流體平均流速的變化曲線Fig.6 Graph of average flow velocity of the fluid at the two-dimensional cross-section with different thicknesses of porous materials

由圖6可見，600 s是一個分界點，在600 s之前，填充多孔材料厚度為1.5 m的流體平均流速最高，這也與速度云圖得到的結果相一致，在600 s之后，流體的自然對流發展到整個儲罐，有多孔材料填充的儲罐內流體的平均流速比無多孔材料填充的要低。

在氣液交界面有蒸發氣體進入到氣相空間時，氣相空間壓力會增大，當其超過安全閥設定壓力時，安全閥開啟，釋放蒸發氣體到大氣空間。當儲罐內液體發生翻滾時，會瞬間產生大量蒸發氣體，不僅會帶來經濟損失還會發生安全事故。圖7為多孔材料不同填充厚度時儲罐內壓力出口的質量流量變化曲線。

圖7 多孔材料不同填充厚度時儲罐內壓力出口的質量流量變化曲線Fig.7 Graph of mass flow rate at the pressure outlet of the storage tank with different thicknesses of porous materials

由圖7可見，不管是哪種情況，儲罐內壓力出口的質量流量都為負值，但填充多孔材料的儲罐均能不同程度地減少壓力出口的質量流量，這樣不僅可減少經濟損失，還可減弱事故發生時的劇烈程度，起到一定的抑制作用。

圖8為多孔材料不同填充厚度時儲罐內上下分層之間的壓差變化曲線。

由圖8可見，在儲罐中心區內填充多孔材料可以減小儲罐內上下分層之間的壓差，即減小儲罐內翻滾事故發生時釋放的能量，從而在一定程度上抑制了翻滾事故的發生。

圖8 多孔材料不同填充厚度時儲罐內上下分層之間的壓差變化曲線Fig.8 Graph of pressure differentials of layers in the storage tank with different thicknesses of porous materials

3 結 論

通過對儲罐中心區無多孔材料填充和有多孔材料填充時儲罐內液化天然氣的流動情況進行數值模擬與分析，得出以下結論：

(1) 在儲罐中心區內填充多孔材料，多孔材料填充部分的流體流速明顯比相鄰區域小，同時在儲罐壁面熱邊界條件和多孔材料作用下，罐體兩邊的流動強度較強，中間的流動強度較弱，中間滾動圈的滯留區面積比兩側滾動圈大，且滾動圈的方向是順時針與逆時針交替出現。

(2) 由于儲罐內液體翻滾的發生需要孕育時間和受計算機的限制，在此前對無多孔材料填充時儲罐內液體翻滾進行模擬時需要很長的時間，因此本文在對儲罐中心區填充多孔材料后儲罐內液體的流動情況進行模擬時并未對整個過程進行模擬，只進行了前1 000 s時的模擬計算及分析，發現儲罐中心區填充多孔材料后，可減小儲罐內流體的平均流速，減少壓力出口的質量流量，并減小儲罐內上下分層之間的壓差。

[1] 林瑞泰.多孔介質傳熱傳質引論[M].北京:科學出版社,1995:39-47.

[2] Nield D A,Bejan A.ConvectioninPorousMedia[M].New York:Springer Velag,1999:23-65.

[3] 邢志祥,杜貞,張成燕,等.密閉儲罐內填充非金屬多孔材料后預混可燃氣體火焰傳播的數值模擬[J].安全與環境學報,2014,14(6):91-95.

[4] 張健中,許光,周金廣,等.網狀鋁合金阻隔防爆材料功效及應用探討[J].中國安全科學學報,2014,24(3):42-46.

[5] 田宏,王旭,高永庭.用于石油液化氣體儲罐填充的多孔金屬材料的防火防爆機理及應用[J].消防技術與產品信息,2000(1):29-30.

[6] 田宏,王旭,高永庭.多孔填充材料的防火防爆機理及應用[J].工業安全與防塵,2000(4):43-46.

[7] 孫福濤,蒲亮,齊迪.大型LNG儲罐分層破壞及翻滾過程的數值研究[J].低溫工程,2017(1):47-53.

[8] 付全榮.泡沫金屬填充套管換熱器內流體流動和傳熱研究[D].太原：太原理工大學,2010.

[9] 于明躍.多孔固體構架與氣流對流換熱特性數值研究[D].哈爾濱：哈爾濱工業大學,2010.

[10] 李泊然.晃蕩條件下LNG液貨艙分層與翻滾現象的數值模擬[D].哈爾濱：哈爾濱工業大學,2016.

[11] 秦允毫.熱學[M].2版.北京:高等教育出版社,2004:154-156.