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

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

0 引言

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

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

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

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

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

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

1 模型描述

1.1 實驗裝置

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

表1 傳感器的具體位置坐標(biāo)Table 1 Specific position coordinates of sensors

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

1.2 守恒方程

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

連續(xù)方程和動量方程分別為：

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

其中，有效剪應(yīng)力的表達式為：

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

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

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

其中，比總焓的定義為：

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

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

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

組分i的質(zhì)量擴散通量矢量的定義為：

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

組分i的質(zhì)量分?jǐn)?shù)需滿足如下守恒方程：

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

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

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

1.3 幾何模型及網(wǎng)格劃分

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

圖2 幾何模型Fig.2 Geometric model

根據(jù)幾何模型的結(jié)構(gòu)特點，本文采用結(jié)構(gòu)化網(wǎng)格，基本單元為六面體，單元最小長度為0.5 mm，網(wǎng)格總數(shù)為213600個。模型泄漏口直徑為0.2 cm，為提高數(shù)值模擬的準(zhǔn)確性和計算穩(wěn)定性，在泄漏口附近的X坐標(biāo)方向和Y坐標(biāo)方向?qū)W(wǎng)格進行了加密。幾何模型的計算域網(wǎng)格劃分如圖3所示。

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

1.4 模型設(shè)置

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

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

2 實驗結(jié)果分析

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

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

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

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

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

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

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

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

圖6 監(jiān)測點處氫氣體積分?jǐn)?shù)隨時間的變化關(guān)系Fig.6 Variation of hydrogen volume fraction at monitoring points with time

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

通過與實驗結(jié)果的對比可以看出，本文所采用的MPPBuoyantEpplFoam 求解器在模擬氫氣泄漏后的流動和分層方面表現(xiàn)出了較高的準(zhǔn)確性。氫氣泄漏停止后，數(shù)值模擬結(jié)果中氫氣的擴散速度比實驗結(jié)果得到的慢，因此需要更長時間達到均勻混合?？傮w來說，MPPBuoyantEpplFoam 求解器的模擬準(zhǔn)確性得到了驗證，后續(xù)可以通過優(yōu)化氫氣擴散模型來提升模擬準(zhǔn)確度。

3 結(jié)論

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

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

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

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

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