受限空間內氫氣泄漏擴散行為的數值模擬

時婷婷，楊福明，常雪倫，孫 晨，于文濤，魯仰輝

(國家電投集團科學技術研究院有限公司，北京 102209)

0 引言

氫能作為21世紀的綠色能源，得到世界各國的普遍關注。隨著氫能的快速發展，氫安全問題正受到越來越多的關注。相比于其他能源形式，氫具有一系列獨特屬性。氫分子是自然界中最小的分子，在高壓系統中很容易通過管道接頭、閥門等位置發生泄漏；氫對很多材料具有氫脆效應，容易導致臨氫環境下的承壓設備失效；氫氣密度低、擴散快，在露天場合不易形成可燃混合物，但是在受限空間內則很容易在局部發生聚集；同時，氫氣具有可燃范圍寬、燃燒熱值高、燃燒速度快、爆炸能量大的特點，聚集的可燃氣云一旦被點燃可能會發生燃爆事故，其后果往往非常嚴重。綜上，在開展氫能及其應用技術的基礎研究過程中，需明確氫安全事故演化行為及后果，探索針對氫安全事故所展開的應急處理措施，為相關法律法規及標準提供可靠的數據支撐，扎實推動氫能技術發展和規模化應用。

泄漏往往是氫安全事故的始發事件，也是后續燃爆行為分析的輸入條件，因此氫氣泄漏擴散行為及演化規律是氫安全研究的基礎性內容。針對受限空間內氫氣泄漏擴散的規律，國內外研究學者利用實驗和數值模擬進行了相應的研究。Denisenko 等[1]和Manins[2]通過實驗研究總結了氫氣在不同形狀、大小和邊界條件的受限空間內的擴散演變機制。實驗結果顯示：在相同泄漏流速的情況下，氫氣在受限空間內的擴散行為可分為2 種，當氫氣的初始泄漏速度小于150 m/s 時，氫氣在受限空間內的擴散方式為填充模式；而當氫氣的初始泄漏速度大于等于150 m/s 時，氫氣在受限空間內的擴散方式變為漸淡模式，從安全角度來講，漸淡模式比填充模式更危險。

de Stefano 等[3]通過實驗分析了受限空間內泄漏位置及泄漏速度對氫氣擴散的影響。實驗結果顯示：氫氣在受限空間內的擴散行為主要取決于泄漏速度而不是泄漏的初始位置，泄漏的初始位置對氫氣擴散的影響隨著泄漏速度的增加而減小。當泄漏速度較小時，氫氣首先聚集在受限空間的底部，然后再上升至天花板，最后下降到其他位置，直至均質化階段。

氦氣的安全性高，可操作性強，且氦氣與氫氣同為低密度氣體，氦氣擴散實驗可作為與在受限空間內使用氫氣相關的安全研究的一部分，在泄漏擴散這一物理過程中，充分了解氦氣這種低密度氣體在空氣中的分散機制至關重要；并且由于氦氣的安全性高，用來模擬氫氣在受限空間內的擴散行為的可操作性相對較強。Cariteau 等[4]通過實驗分析了在1 m×1 m×1 m 的受限空間內3 種不同通風口形狀下氦氣的擴散行為。實驗結果顯示：氦氣的擴散行為受通風口形狀的影響較大；自然通風會顯著改變低密度氣體的擴散行為。

根據上述研究結論，可在氫能基礎設施中設置便于通風的開口。Pitts 等[5-6]、Blais 等[7]通過13 組1∶1 全尺寸實驗，研究了車庫內氫氣泄漏擴散及燃燒的發展過程。實驗結果為氫燃料電池車輛相關標準和規范的開發提供了氫相關安全數據，以減輕氫燃料電池的火災隱患。隨后，Hajji 等[8]通過Fluent 軟件對文獻[5-7]研究的車庫內氫氣泄漏實驗進行了模擬，對比驗證了泄漏位置、釋放速率、泄漏持續時間和氫氣濃度之間的關系，模擬結果與實測結果吻合良好。

準確的數值結果可為氫能設施的傳感器布設、結構優化和通風裝置設計等提供參考，為此，基于OpenFOAM[9]、EPPL 算法[10]和多組分庫Mutation++[11]，本項目研發團隊研發了魯棒、高效、精確的全速域多組分浮力流求解器——MPPBuoyantEpplFoam 求解器。本文采用MPPBuoyantEpplFoam 求解器對氫氣在受限空間內(即室內環境)的泄漏擴散行為進行了數值模擬，對氫氣泄漏后的流動擴散規律開展了研究，并對其多組分擴散功能的準確性進行了驗證。

1 模型描述

1.1 實驗裝置

本文以de Stefano 等[3]的實驗裝置作為模擬對象。實驗場景為某核電廠內易發生氫氣泄漏的隔間，為了便于調試氫氣流速，實驗裝置為該隔間1/15 縮比模型；實驗裝置的形狀為立方體，長×寬×高的尺寸為0.47 m×0.33 m×0.20 m；泄漏口位于裝置頂部。為了探究不同位置氫氣濃度隨時間的變化關系，實驗裝置內設有16個傳感器，分別固定在A、B、C、D 4 根支架上，每根支架上有4個傳感器。實驗裝置的尺寸及布置的傳感器位置示意圖如圖1所示，傳感器的具體位置坐標如表1所示。

表1 傳感器的具體位置坐標Table 1 Specific position coordinates of sensors

圖1 實驗裝置尺寸及布置的傳感器位置示意圖Fig.1 Diagram of size of experimental device and sensor location

1.2 守恒方程

MPPBuoyantEpplFoam 求解器采用了先進的算法模型，具有計算準確、穩定性高、網格適應性好等優點；同時由于其可以采用較大的CFL數進行模擬，整體計算效率較高；此外，該求解器還引進了最新開發的全速域求解算法，可以開展亞聲速到超聲速的工況模擬。該求解器采用多組分瞬態可壓縮流動的控制方程，包括連續方程、動量方程、能量方程和組分方程。

連續方程和動量方程分別為：

式中：ρ為密度，kg/m3；t為時間，s；V為速度，m/s；p為壓強，Pa；g為重力加速度，m/s2；τeff為有效剪應力，Pa。

其中，有效剪應力的表達式為：

式中：μ為混合物的分子黏度，Pa·s；μt為混合物的湍流動力黏度，Pa·s。

對于多組分流，比總焓方程形式如下：

式中：H為比總焓，J/kg；D為壓力隨時間變化的擴散系數，m2/s；q為熱擴散通量矢量，kg/(m2·s)。

其中，比總焓的定義為：

式中：h為比焓，J/kg。

熱擴散通量矢量的定義為：

式中：λ為混合物的平均導熱系數，W/(m·K)；Cv為定容比熱容，J/(kg·K)；Prt為湍流普朗特數；T為氣體溫度，K；Ns為組分的個數；hi為組分i的比焓，J/kg；ji為組分i的質量擴散通量矢量，kg/(m2·s)。

組分i的質量擴散通量矢量的定義為：

式中：Dim為組分i擴散到其他組分的平均擴散系數，計算方法參考文獻[12]；Sct為湍流施密特數；Yi為組分i的質量分數。

組分i的質量分數需滿足如下守恒方程：

此外，假定壓強和密度滿足理想氣體狀態方程，即：

式中：R為摩爾氣體常數，J/(kmol·K)；Wi為組分i的摩爾分子質量，kg/kmol。

以上各式中，混合物的分子黏度、混合物的平均導熱系數和組分i擴散到其他組分的平均擴散系數的取值與混合物的溫度和組成有關，通過多組分庫Mutation++計算；混合物的湍流動力黏度等湍流參數與湍流模型相關，由于考慮的實驗條件受湍流的影響相對較弱，本文暫不予以考慮。

1.3 幾何模型及網格劃分

根據前述實驗裝置的尺寸和傳感器的布置情況，建立幾何模型并劃分網格。幾何模型如圖2所示。

圖2 幾何模型Fig.2 Geometric model

根據幾何模型的結構特點，本文采用結構化網格，基本單元為六面體，單元最小長度為0.5 mm，網格總數為213600個。模型泄漏口直徑為0.2 cm，為提高數值模擬的準確性和計算穩定性，在泄漏口附近的X坐標方向和Y坐標方向對網格進行了加密。幾何模型的計算域網格劃分如圖3所示。

圖3 幾何模型的計算域網格劃分Fig.3 Mesh generation in computational domain of geometric model

1.4 模型設置

設定初始狀態下模型內充滿靜止空氣，初始壓力為1 atm(1 atm=101325 Pa)(絕對壓強)，初始溫度為298 K。氫氣從模型頂部泄漏口豎直向下注入計算域，泄漏速度為2.2 m/s，氫氣注入持續時間為60 s。本模型考慮重力影響，重力方向為豎直向下，大小為9.8 m/s2。模型總計算時間為220 s，時間步長為0.0001 s，結果輸出間隔為1 s。

為了探究不同位置氫氣濃度隨時間的變化關系并與文獻[3]中的實驗數據進行對比，按照表1所示的傳感器具體位置坐標在模型中設置了氫氣體積分數監測點。

2 實驗結果分析

氫氣自模型頂部泄漏后，在計算域內泄漏口處的泄漏速度切片如圖4所示。

圖4 泄漏口處的泄漏速度切片圖Fig.4 Slice diagram of leakage velocity at leakage port

從圖4可知：氫氣從泄漏口豎直向下噴出，泄漏速度為2.2 m/s。噴出后擴散速度迅速減小，這主要受2個方面的影響：一方面，氫氣射流受到周圍靜止空氣的阻滯作用，發生動量交換；另一方面，氫氣密度是空氣密度的1/14，離開噴射口后的氫氣受到較強的向上浮升力作用，從而使向下的擴散速度迅速降低。

氫氣開始泄漏后，不同時間下體積分數為4%時的氫氣濃度云圖如圖5所示，監測點處氫氣體積分數隨時間的變化關系如圖6所示。圖中：Cmax為最大氫氣體積分數；th為中間的間隔時間；Δa為受限空間內上部氫氣與下部氫氣的濃度差，泄漏剛停止時的濃度差為12.00%，泄漏停止20 s 后的濃度差為8.03%。

從圖5和圖6可以看出：氫氣泄漏后的整個過程可以分為注入階段、擴散階段和均勻混合階段3個階段。在注入階段，氫氣在初始動量的作用下快速在計算域中蔓延，氫氣濃度迅速增加，并且在浮力作用下呈現分層狀態；氫氣注入結束，模型上部氫氣濃度達到最大值，此時模型上下部氫氣濃度差也達到最大值；此后進入擴散階段，上下部氫氣濃度差在分子擴散作用下逐漸減少；最終系統達到氫氣和空氣均勻混合的穩定狀態，即均勻混合階段。

圖5 氫氣開始泄漏后，不同時間下體積分數為4%時的氫氣濃度云圖Fig.5 Cloud diagram of hydrogen concentration when volume fraction is 4% at different time after hydrogen leakage

文獻[3]中，傳感器測得的實驗裝置中氫氣體積分數隨時間的變化關系如圖7所示。圖中：1、2、3 分別代表注入階段、擴散階段和均勻混合階段；Δc為受限空間內上部氫氣與下部氫氣的濃度差。

對比圖6和圖7可以看出：圖6的數值模擬結果與圖7的實驗結果曲線趨勢基本一致。氫氣泄漏后，由于氫氣密度較低，其濃度分布在空間中表現出明顯的分層效應。泄漏停止后，氫氣在擴散作用下分層逐漸減弱，實驗裝置上部的氫氣濃度水平降低，而下部的濃度水平升高，最后達到均勻混合。從定量上來看，數值模擬獲得的最大氫氣體積分數為12.77%，實驗測得的最大氫氣體積分數為13.75%，兩者僅相差0.98%。泄漏停止時，數值模擬結果和實驗結果均顯示實驗裝置上部和下部的氫氣濃度差為12%；20 s 后，數值模擬結果和實驗結果均顯示實驗裝置上部和下部的氫氣濃度差為8%。實驗結果顯示，在停止泄漏后約130 s，實驗裝置內的氫氣與空氣達到均勻混合；而數值模擬結果在130 s 后，實驗裝置內的氫氣和空氣還未達到均勻混合，這是由于MPPBuoyantEpplFoam 求解器完全使用層流模型來進行擴散的計算分析，而在實際實驗中很有可能受到了湍流的影響，因此氫氣和空氣混合的速度較快。

圖6 監測點處氫氣體積分數隨時間的變化關系Fig.6 Variation of hydrogen volume fraction at monitoring points with time

圖7 傳感器測得的實驗裝置中氫氣體積分數隨時間的變化關系[3]Fig.7 Variation of hydrogen volume fraction in the experimental device with time measured by sensor[3]

通過與實驗結果的對比可以看出，本文所采用的MPPBuoyantEpplFoam 求解器在模擬氫氣泄漏后的流動和分層方面表現出了較高的準確性。氫氣泄漏停止后，數值模擬結果中氫氣的擴散速度比實驗結果得到的慢，因此需要更長時間達到均勻混合。總體來說，MPPBuoyantEpplFoam 求解器的模擬準確性得到了驗證，后續可以通過優化氫氣擴散模型來提升模擬準確度。

3 結論

本文采用自主開發的MPPBuoyantEpplFoam求解器對受限空間內氫氣泄漏擴散行為進行了數值模擬。結果顯示，該求解器模擬得到的氫氣濃度分布及變化情況與文獻中的實驗結果吻合較好，初步驗證了多組分擴散功能的準確性。在本文所研究的受限空間內，氫氣泄漏后表現出如下流動擴散規律：

1)氫氣從模型頂部豎直向下噴入計算域后，在周圍靜止空氣阻滯作用和浮升力作用下，擴散速度迅速減小，流動的動量控制區變短。

2)在氫氣注入階段，計算域氫氣濃度分布會在浮力作用下呈現明顯的分層效應，上部的氫氣濃度高，下部的氫氣濃度低。

3)氫氣停止注入后，在分子擴散作用下，計算域上下部氫氣濃度差逐漸變小，并最終達到均勻混合。

致謝：感謝哈爾濱工程大學的張慧杰博士在本文工作上給予的幫助和指導！