邢志祥,張淑淑,汪李金,張 瑩
(常州大學(xué)環(huán)境與安全工程學(xué)院,江蘇 常州 213164)
多孔介質(zhì)是一種多相物質(zhì)共存、由固相和氣相或液相組成的物質(zhì),固相作為固體骨架,氣相或液相充滿了孔隙空間,孔隙之間必須相互連通[1]。Nield等[2]最先開始了多孔介質(zhì)的研究,認為多孔材料具有比表面積大、質(zhì)量輕、體積小、導(dǎo)熱性能好的特點。近年來,針對多孔介質(zhì)的試驗和模擬研究很多,研究主要集中在填充多孔材料后液體流動的情況以及多孔材料內(nèi)部火焰?zhèn)鞑サ那闆r,而研究對象通常是管道、密閉容器等,但關(guān)于儲罐內(nèi)填充多孔材料的研究則相對較少。邢志祥等[3]模擬了圓柱形儲罐內(nèi)填充聚氨酯多孔材料后可燃氣體火焰?zhèn)鞑サ那闆r,并與13 L圓柱形密閉儲罐現(xiàn)場試驗的結(jié)果進行了對比,結(jié)果表明:模擬結(jié)果與試驗結(jié)果相吻合,填充物平均孔徑越小,阻火抑爆性能越好;張健中等[4]分析了在加油站埋地油罐中是否有必要使用網(wǎng)狀鋁合金阻隔防爆材料等問題;田宏等[5-6]對填充在液化石油氣儲罐中的多孔金屬材料做了系統(tǒng)的介紹,但并未通過試驗來舉證。
液化天然氣(Liquified Natural Gas,LNG)儲罐是儲存液化天然氣的重要設(shè)備,其運行處于超低溫狀態(tài)時容易引發(fā)各類安全事故,而LNG翻滾是其中一種事故類型。儲罐內(nèi)由于LNG密度差產(chǎn)生分層,上層液體與下層液體之間存在液液分界面,而漏熱導(dǎo)致分界面被破壞進一步引發(fā)翻滾現(xiàn)象。分界面被破壞的原因有兩種[7]:一種是儲罐內(nèi)壁的邊界層穿透,導(dǎo)致分界面被破壞;另一種是上層液體在中心處向下的射流和下層液體在中心處向下的卷攜,導(dǎo)致中心處的分界面被破壞,進一步造成了整個分界面被破壞??紤]到多孔材料的優(yōu)點,本文將其應(yīng)用到儲罐中研究其是否能起到抑制LNG翻滾的作用,并選取儲罐中心區(qū)填充多孔材料,對儲罐內(nèi)LNG的流動特征進行了數(shù)值模擬。
本文選取儲罐中心區(qū)填充多孔材料,并分別設(shè)置多孔材料的填充厚度為1 m、1.5 m、2 m,通過建立儲罐中心區(qū)填充多孔材料的物理模型和數(shù)學(xué)模型,對儲罐內(nèi)LNG的流動情況進行了數(shù)值模擬計算。
本文選取6 000 m3的圓柱形儲罐,內(nèi)罐直徑為24 m,高度為17.4 m,設(shè)計液位高15.76 m,并取儲罐的軸截面建立多孔介質(zhì)中湍流流動的二維模型。坐標原點位于左側(cè)壁面和底部壁面的交匯處,X軸向右為正,Y軸向上為正;模擬過程不考慮氣相空間,只針對儲罐內(nèi)的液相部分來模擬LNG的流動情況,液體充裝高度為14 m;多孔材料填充在儲罐的中心部位;對儲罐結(jié)構(gòu)進行了簡化,不考慮壁厚,并忽略儲罐內(nèi)部的管線、循環(huán)泵等部件,由此建立的儲罐中心區(qū)填充多孔材料簡化后的物理模型如圖1所示。
圖1 儲罐中心區(qū)填充多孔材料簡化后的物理模型Fig.1 Simplified physical model of the filling porous material of the storage tank central area
對儲罐中心區(qū)填充多孔材料時儲罐內(nèi)LNG流動情況的物理模型(中間標明H的區(qū)域為多孔材料區(qū)域,左右兩側(cè)區(qū)域為純流體區(qū)域)進行了如下的簡化:①儲罐內(nèi)部流體分成兩個高度相等的分層,分層高度為7 m,上分層流體密度為424 kg/m3,下分層流體密度為425 kg/m3,初始密度差為1 kg/m3,兩個分層內(nèi)的密度均勻一致;②儲罐內(nèi)流體密度符合Boussinesq假設(shè);③選取金屬鋁多孔材料,且視多孔材料為均勻各向同性;④多孔材料的孔隙為球形空洞,孔隙之間相互連通,即為開孔結(jié)構(gòu);⑤多孔材料的物性參數(shù)為常數(shù);⑥多孔材料固體骨架與流體處于局部熱平衡,且無任何化學(xué)反應(yīng);⑦儲罐內(nèi)流體流動為湍流流動;⑧忽略流體流動過程中的黏性耗散;⑨不考慮輻射換熱。
采用非穩(wěn)態(tài)模型、VOF模型和標準k-ε湍流模型,并基于Boussinesq假設(shè),計算過程中多孔區(qū)域動量方程采用Darcy-Brinkman-Forchheimer模型,能量方程采用局部熱平衡模型,得到的控制方程如下:
連續(xù)性方程:
?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系數(shù);Cp為定壓比熱容[J/(kg·K)];T為溫度(K);λeff為多孔區(qū)域的有效導(dǎo)熱系數(shù)[W/(m·K)],由流體的導(dǎo)熱系數(shù)λf和固體的導(dǎo)熱系數(shù)λs的體積平均值計算得到,即λeff=(1-φ)λs+φλf;k為湍動能(J);ε為湍流擴散率;σk和σε分別為k方程和ε方程的普朗特數(shù),σk=1.0,σε=1.3;μt為湍流黏性系數(shù),μt=Cμρk2/ε,其中Cμ=0.09;Gk為平均速度梯度引起的湍動能k的產(chǎn)生項;Gb為浮力引起的湍動能k的產(chǎn)生項;YM為可壓縮湍流流動中脈動擴張的貢獻;C1ε、C2ε、C3ε為經(jīng)驗常數(shù),C1ε=1.44,C2ε=1.92,C3ε=1;Sk和Sε為用戶自定義源項。
當φ=1時,k→∞,λeff=λf,動量方程的源項趨近于0,能量方程也變?yōu)闃藴誓芰糠匠?,則上述控制方程變?yōu)榧兞黧w區(qū)域的控制方程。
儲罐側(cè)壁面和底壁面取為無滑移壁面邊界條件,并采用定熱流的方式加熱,熱流密度恒定為30 W/m2;氣液交界面取為壓力出口;多孔區(qū)域采用金屬鋁多孔材料,孔隙率設(shè)定為0.95,并需設(shè)定黏性阻力系數(shù)和慣性阻力系數(shù);多孔區(qū)域與純流體區(qū)域的交界面設(shè)為內(nèi)部邊界。
初始時刻:初始表壓為15 000 Pa,X和Y軸方向的速度分量都為0,分層區(qū)初始溫度為111 K,主流區(qū)初始溫度為111.5 K。
在計算過程中,多孔區(qū)域和純流體區(qū)域通過φ的取值進行區(qū)分,并統(tǒng)一進行求解。
(1) 孔隙率φ和孔徑dp[8]:孔隙率φ是指孔隙體積占多孔材料總體積的比值;孔密度是指多孔材料每英寸上孔的數(shù)目(1 inch=0.025 4 m),孔徑dp=0.025 4/孔密度(m)。
(2) 滲透率K和Forchheimer系數(shù)CF[9]:滲透率K代表了多孔介質(zhì)中孔隙的表面積和彎曲程度,通常由Vafai總結(jié)的經(jīng)驗公式確定,即
當多孔材料中的流體流動為湍流流動時,達西定律便不能描述流體流動的情況,此時可采用其修正模型進行描述,即在動量方程中加入源項,源項由兩項組成,一項為黏性阻力項,一項為慣性阻力項,慣性阻力項的系數(shù)用CF表示,稱為Forchheimer系數(shù),可表示為
CF=1.75150φ3/2
(3) 壓差Δp:壓差指儲罐內(nèi)上下分層之間的壓力差(Pa),Δp=p2-p1,其中p2為底壁面平均壓力(Pa),p1為壓力出口平均壓力(Pa)。
本文選取儲罐的軸截面建立二維模型,應(yīng)用Fluent 6.3軟件對儲罐中心區(qū)無多孔材料填充和有多孔材料填充時儲罐內(nèi)LNG的流動進行了數(shù)值模擬,模擬斷面均取自儲罐的同一軸截面,并選取速度云圖和流線圖的模擬結(jié)果進行了分析。速度云圖是描述流體速度大小分布情況的圖;流線圖是描述某一時刻流體運動趨勢的圖,可以通過流線圖中的流線方向和流線疏密,判斷出某一時刻流體的運動方向和區(qū)分流動強弱區(qū)域,流線圖中箭頭表示流體的運動方向,流線的疏密反映流速大小,流線越密,流體流速越大,反之亦然。
圖2為儲罐中心區(qū)無多孔材料填充時,t=100 s和t=1 000 s時刻儲罐內(nèi)LNG的流動情況模擬結(jié)果。
圖2 儲罐中心區(qū)無多孔材料填充時t=100 s和t=1 000 s 時刻的儲罐內(nèi)流體的速度云圖和流線圖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時,由于外界環(huán)境的漏熱,儲罐壁面處的流體率先吸收熱量導(dǎo)致密度減小,產(chǎn)生浮力驅(qū)動流,沿著壁面向上移動,形成流動邊界層;由于上下兩層流體存在溫度差和密度差,因此位于分界面處的流體出現(xiàn)自然對流現(xiàn)象,且流體流速較儲罐其他位置的流速大,而分界面處邊界層的流體流速最大[見圖2(a)]。
(2) 在t=1 000 s時,分界面處流體的自然對流發(fā)展到整個儲罐內(nèi)流體的自然對流,且上下兩層流體的自然對流各自獨立,即由液液分界面隔開,被限制在自身區(qū)域內(nèi)。對于上層流體,由于儲罐側(cè)壁面的漏熱,導(dǎo)致邊界層處的流體受熱膨脹密度減小,沿儲罐壁面向上移動,當其到達氣液交界面時,其中一部分流體蒸發(fā)變成氣體釋放到氣相空間,其密度增大導(dǎo)致向下移動,由此形成了上層流體的自然能對流現(xiàn)象;對于下層流體,邊界層處的流體受熱膨脹密度減小,沿著邊界層向上移動,當其到達液液分界面時,因浮力太小而無法穿過分界面進入上層,但會通過分界面向上層流體傳遞熱量,因而溫度降低,密度增大,導(dǎo)致向下移動,由此形成了下層流體的自然對流現(xiàn)象[見圖2(b)]。
(3) 在純自然對流運動中,流體的瑞利數(shù)Ra是判定由浮力產(chǎn)生的對流強度大小的標準[10]。當Ra小于臨界值時,流體之間是熱傳導(dǎo)狀態(tài),不發(fā)生對流運動;當Ra大于臨界值時,才會發(fā)生對流運動,從而在液體中出現(xiàn)宏觀對流花紋,稱其為Benard花紋[11]。儲罐內(nèi)無多孔材料填充時,在t=1 000 s時,上層和下層的自然對流結(jié)構(gòu)明顯不同,但流動都是由類似圓形的滾動圈組成,上層的滾動圈數(shù)為1個,為扁長形,而下層的滾動圈數(shù)為4個,分別為2個胖圓形和2個扁長形的滾動圈;相比于t=100 s時分界面處的流體流速較高,t=1 000 s時分界面處的流體流速反而較小,除此之外,在滾動圈中心處的流速也較小,稱之為滯留區(qū)[見圖2(c)]。
在模擬計算過程中,對二維截面上流體的最大流速進行了監(jiān)測,得到流體的最大流速為0.195 m/s。
圖3為儲罐中心區(qū)多孔材料填充厚度為1 m時,t=100 s和t=1 000 s時刻儲罐內(nèi)LNG的流動情況模擬結(jié)果。
圖3 儲罐中心區(qū)多孔材料填充厚度為1 m時t=100 s和t=1 000 s時刻的儲罐內(nèi)流體的速度云圖和流線圖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時,儲罐內(nèi)流體流動發(fā)展情況相對于無多孔材料填充時比較滯后,分界面處的流體流速較無多孔材料填充時小,除邊界層外,多孔材料填充區(qū)域兩側(cè)自然對流強度較高[見圖3(a)]。
(2) 在t=1 000 s時,多孔材料填充區(qū)域相比相鄰區(qū)域的流體流速小[見圖3(b)];相對于無多孔材料填充時,上層和下層的自然對流結(jié)構(gòu)有明顯不同,上層的自然對流結(jié)構(gòu)由1個扁長形滾動圈變?yōu)閹讉€眼淚形狀的小滾動圈,且由于多孔材料增加了流體流動阻力,在多孔材料填充區(qū)域左側(cè)流體流線豎直向下,穿過多孔材料后,便向氣液交界面移動,而下層的自然對流結(jié)構(gòu)由2個胖圓形和2個扁長形的滾動圈變?yōu)?個胖圓形的滾動圈,且中間的滾動圈比左右兩側(cè)的滾動圈大,滾動圈方向從左到右依次為逆時針、順時針、逆時針,同時發(fā)現(xiàn)在壁面熱邊界條件和多孔材料作用下,罐體兩邊的流動強度較強,中間的流動強度較弱,中間滾動圈的滯留區(qū)面積比兩側(cè)滾動圈大[見圖3(c)]。
儲罐內(nèi)沿中心區(qū)填充多孔材料厚度為1 m時,二維截面上流體的最大流速為0.190 m/s。
圖4為儲罐中心區(qū)多孔材料填充厚度為1.5 m時,t=100 s和t=1 000 s時刻儲罐內(nèi)LNG的流動情況模擬結(jié)果。
圖4 儲罐中心區(qū)多孔材料填充厚度為1.5 m時t=100 s和 t=1 000 s時刻的儲罐內(nèi)流體的速度云圖和流線圖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時,仍是邊界層處流體和分界面處流體的流速較儲罐其他位置的流速大,分界面處出現(xiàn)自然對流現(xiàn)象,但與同時期無多孔材料填充時相比,分界面上層流體的流速較高,在多孔材料填充區(qū)域兩側(cè)各出現(xiàn)了3個“火焰”形的高流速滾動圈,其速度分布同“火焰”的溫度分布一樣,中心速度最大,向四周減弱[見圖4(a)]。
(2) 在t=1 000 s時,多孔材料填充區(qū)域的流體流速最小,且下層的自然對流結(jié)構(gòu)明顯沿儲罐中心線呈對稱分布[見圖4(b)];上層和下層的自然對流結(jié)構(gòu)與無多孔材料填充時有明顯的不同,上層的自然對流結(jié)構(gòu)與多孔材料填充厚度為1 m時一樣,仍為幾個大小不一、眼淚形狀的小滾動圈,而且在多孔材料填充區(qū)域內(nèi)部有一個滾動圈,而下層的自然對流結(jié)構(gòu)沿儲罐中心線呈對稱分布,由2個近似相等的胖圓形和2個扁長形的滾動圈組成,且中間的滾動圈比左右兩側(cè)的滾動圈大,中心線兩側(cè)的滾動圈方向依次為順時針、逆時針、順時針、逆時針,順時針與逆時針交替出現(xiàn),同時發(fā)現(xiàn)與多孔材料填充厚度為1 m時相似,中間滾動圈的滯留區(qū)面積比兩側(cè)滾動圈大,兩側(cè)滾動圈的滯留區(qū)面積幾乎為0[見圖4(c)]。
儲罐內(nèi)沿中心區(qū)填充多孔材料厚度為1.5 m時,二維截面上流體的最大流速為0.371 m/s。
圖5為儲罐中心區(qū)多孔材料填充厚度為2 m時,t=100 s和t=1 000 s時刻儲罐內(nèi)LNG的流動情況模擬結(jié)果。
圖5 儲罐中心區(qū)多孔材料填充厚度為2 m時t=100 s和t=1 000 s時刻的儲罐內(nèi)流體的速度云圖和流線圖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時,儲罐內(nèi)流體流動發(fā)展情況與同時期多孔材料填充厚度為1 m時較相似,分界面處的流體流速較無多孔材料填充時小,除邊界層外,多孔材料填充區(qū)域兩側(cè)自然對流強度較高[見圖5(a)]。
(2) 在t=1 000 s時,相對于無多孔材料填充時,低速區(qū)域范圍變大,高速區(qū)域范圍變小,且速度場與填充厚度為1.5 m時的速度場較相似[見圖5(b)];上層的自然對流結(jié)構(gòu)為2個小的滾動圈,下層的自然對流結(jié)構(gòu)主要為2個大的滾動圈,滾動圈方向依次為順時針、逆時針,此時的流場結(jié)構(gòu)不再沿儲罐中心線呈對稱分布,而是以中心線為分界線,左右兩側(cè)流場具有相似性,說明多孔材料填充厚度為2 m時,極大地增加了流動阻力,好像一堵墻,起到了分割流場的作用[見圖5(c)]。
儲罐內(nèi)沿中心區(qū)填充多孔材料厚度為2 m時,二維截面上流體的最大流速為0.189 m/s。
圖6為多孔材料不同填充厚度(H=0、1 m、1.5 m、2 m)時儲罐內(nèi)二維截面上流體平均流速的變化曲線。
圖6 多孔材料不同填充厚度時儲罐內(nèi)二維截面上流體平均流速的變化曲線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的流體平均流速最高,這也與速度云圖得到的結(jié)果相一致,在600 s之后,流體的自然對流發(fā)展到整個儲罐,有多孔材料填充的儲罐內(nèi)流體的平均流速比無多孔材料填充的要低。
在氣液交界面有蒸發(fā)氣體進入到氣相空間時,氣相空間壓力會增大,當其超過安全閥設(shè)定壓力時,安全閥開啟,釋放蒸發(fā)氣體到大氣空間。當儲罐內(nèi)液體發(fā)生翻滾時,會瞬間產(chǎn)生大量蒸發(fā)氣體,不僅會帶來經(jīng)濟損失還會發(fā)生安全事故。圖7為多孔材料不同填充厚度時儲罐內(nèi)壓力出口的質(zhì)量流量變化曲線。
圖7 多孔材料不同填充厚度時儲罐內(nèi)壓力出口的質(zhì)量流量變化曲線Fig.7 Graph of mass flow rate at the pressure outlet of the storage tank with different thicknesses of porous materials
由圖7可見,不管是哪種情況,儲罐內(nèi)壓力出口的質(zhì)量流量都為負值,但填充多孔材料的儲罐均能不同程度地減少壓力出口的質(zhì)量流量,這樣不僅可減少經(jīng)濟損失,還可減弱事故發(fā)生時的劇烈程度,起到一定的抑制作用。
圖8為多孔材料不同填充厚度時儲罐內(nèi)上下分層之間的壓差變化曲線。
由圖8可見,在儲罐中心區(qū)內(nèi)填充多孔材料可以減小儲罐內(nèi)上下分層之間的壓差,即減小儲罐內(nèi)翻滾事故發(fā)生時釋放的能量,從而在一定程度上抑制了翻滾事故的發(fā)生。
圖8 多孔材料不同填充厚度時儲罐內(nèi)上下分層之間的壓差變化曲線Fig.8 Graph of pressure differentials of layers in the storage tank with different thicknesses of porous materials
通過對儲罐中心區(qū)無多孔材料填充和有多孔材料填充時儲罐內(nèi)液化天然氣的流動情況進行數(shù)值模擬與分析,得出以下結(jié)論:
(1) 在儲罐中心區(qū)內(nèi)填充多孔材料,多孔材料填充部分的流體流速明顯比相鄰區(qū)域小,同時在儲罐壁面熱邊界條件和多孔材料作用下,罐體兩邊的流動強度較強,中間的流動強度較弱,中間滾動圈的滯留區(qū)面積比兩側(cè)滾動圈大,且滾動圈的方向是順時針與逆時針交替出現(xiàn)。
(2) 由于儲罐內(nèi)液體翻滾的發(fā)生需要孕育時間和受計算機的限制,在此前對無多孔材料填充時儲罐內(nèi)液體翻滾進行模擬時需要很長的時間,因此本文在對儲罐中心區(qū)填充多孔材料后儲罐內(nèi)液體的流動情況進行模擬時并未對整個過程進行模擬,只進行了前1 000 s時的模擬計算及分析,發(fā)現(xiàn)儲罐中心區(qū)填充多孔材料后,可減小儲罐內(nèi)流體的平均流速,減少壓力出口的質(zhì)量流量,并減小儲罐內(nèi)上下分層之間的壓差。
[1] 林瑞泰.多孔介質(zhì)傳熱傳質(zhì)引論[M].北京:科學(xué)出版社,1995:39-47.
[2] Nield D A,Bejan A.ConvectioninPorousMedia[M].New York:Springer Velag,1999:23-65.
[3] 邢志祥,杜貞,張成燕,等.密閉儲罐內(nèi)填充非金屬多孔材料后預(yù)混可燃氣體火焰?zhèn)鞑サ臄?shù)值模擬[J].安全與環(huán)境學(xué)報,2014,14(6):91-95.
[4] 張健中,許光,周金廣,等.網(wǎng)狀鋁合金阻隔防爆材料功效及應(yīng)用探討[J].中國安全科學(xué)學(xué)報,2014,24(3):42-46.
[5] 田宏,王旭,高永庭.用于石油液化氣體儲罐填充的多孔金屬材料的防火防爆機理及應(yīng)用[J].消防技術(shù)與產(chǎn)品信息,2000(1):29-30.
[6] 田宏,王旭,高永庭.多孔填充材料的防火防爆機理及應(yīng)用[J].工業(yè)安全與防塵,2000(4):43-46.
[7] 孫福濤,蒲亮,齊迪.大型LNG儲罐分層破壞及翻滾過程的數(shù)值研究[J].低溫工程,2017(1):47-53.
[8] 付全榮.泡沫金屬填充套管換熱器內(nèi)流體流動和傳熱研究[D].太原:太原理工大學(xué),2010.
[9] 于明躍.多孔固體構(gòu)架與氣流對流換熱特性數(shù)值研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2010.
[10] 李泊然.晃蕩條件下LNG液貨艙分層與翻滾現(xiàn)象的數(shù)值模擬[D].哈爾濱:哈爾濱工業(yè)大學(xué),2016.
[11] 秦允毫.熱學(xué)[M].2版.北京:高等教育出版社,2004:154-156.