摘 要：液氫泄漏后產(chǎn)生的易燃易爆氫氣云團存在嚴重安全隱患。通過基于OpenFOAM開發(fā)的求解器，對液氫泄漏及云團擴散行為進行分析。該求解器采用簡化蒸發(fā)模型確定液氫蒸發(fā)速率，通過氣相和固相求解器分別模擬云團的擴散以及地基的傳熱，使用共軛傳熱方法保證氣-固區(qū)域間溫度和熱通量的連續(xù)性，并采用考慮大氣穩(wěn)定度的對數(shù)風廓線作為大氣邊界條件。通過風廓線觀測數(shù)據(jù)和液氫泄漏實驗數(shù)據(jù)驗證提出的求解器。研究結果表明，穩(wěn)定大氣條件、高溫以及大泄漏量均導致可燃氫氣云團橫向擴散范圍增大，采用滯留坑方案可將蒸發(fā)速率降低99%，并使可燃云團的擴散距離及高度分別降低46%和76%，有效縮小了危險范圍。

關鍵詞：氫；液體泄漏；大氣邊界層；擴散；數(shù)值模擬

中圖分類號：TK91 """"""""""""""""""""""""""""""""""文獻標志碼：A

0 引 言

氫由于低污染和高能量密度被廣泛認為是一種具有巨大潛力的能源載體。低溫液態(tài)儲氫是氫氣儲存和運輸?shù)姆绞街唬湓诖笠?guī)模以及長距離運輸方面存在一定優(yōu)勢，未來有望與高壓氣態(tài)儲氫互補發(fā)展，為氫能產(chǎn)業(yè)提供重要支持。然而，一旦液氫儲罐破裂，將會產(chǎn)生嚴重安全隱患。大規(guī)模的液氫泄漏將形成液氫池［1］，液氫池吸收周圍環(huán)境如地面、空氣和太陽輻射的熱量，劇烈沸騰并蒸發(fā)。蒸發(fā)的氫氣會以氣態(tài)形式擴散到大氣環(huán)境中，形成云狀結構。這種擴散會形成易燃和易爆的云團，對公眾構成嚴重的安全威脅。因此，準確預測液氫泄漏及云團的擴散行為至關重要。

目前液氫泄漏實驗較少，且主要針對于小規(guī)模、封閉或半封閉的泄漏情形［2-5］，鮮有研究者對危害范圍更大的大氣環(huán)境下大規(guī)模液氫泄露情形開展實驗研究。文獻［6］在兩個建筑物之間進行了一系列液氫泄露實驗，發(fā)現(xiàn)建筑物和風向顯著影響氫云的擴散行為；文獻［7］開展了稍小規(guī)模的液氫泄漏實驗，分析了液氫泄漏高度和方向對后續(xù)擴散的影響；文獻［8］開展了目前最具代表性的大氣環(huán)境液氫泄漏實驗研究，其被廣泛用于數(shù)值模擬研究的模型驗證［9-10］。然而由于大氣穩(wěn)定度、風向等一系列復雜因素的影響，文獻［8］開展的7組實驗中僅第6組成功監(jiān)測到局部區(qū)域的氫氣濃度以及可視氫氣云團的最遠擴散距離、最大高度和消失時間，此外該組實驗未能捕獲除檢測塔以外云團擴散的過程。

鑒于實驗成本、監(jiān)測難度和安全隱患等諸多限制，眾多學者采用CFD方法對該問題進行研究。針對液氫泄漏云團擴散行為，研究者提出考慮多種因素的數(shù)值模型，包括相變［11-12］、濕度［13］和風向波動［14］。唐鑫等［15］采用考慮多組分相變的三維瞬態(tài)模型模擬了不同環(huán)境條件下液氫蒸發(fā)和擴散行為，并分析了氫氣云團的形態(tài)變化、持續(xù)時間和危險范圍；Hansen［16］和蒲亮等［11］發(fā)現(xiàn)低溫氫氣云團在近場區(qū)域與稠密氣體行為相似，而在遠場區(qū)域迅速上升；Giannissi等［12］的計算結果表明空氣中氮氣、氧氣的相變增加了氫氣云團的浮力；文獻［10，13］使用Ansys Fluent研究風速、風溫、空氣濕度和地面溫度對氫云濃度分布的影響；邵翔宇等［17］研究了不同天氣條件下氫氣云團擴散行為，發(fā)現(xiàn)風速的增加或大氣壓的降低促進了氫氣云團與周圍空氣的混合。

然而，當前的數(shù)值研究缺乏對大氣穩(wěn)定度和泄漏源等因素的研究。實際上，大氣穩(wěn)定度對風廓線有著顯著影響［18-19］，并將進一步影響云團的擴散。此外泄漏源對應的入口條件同樣至關重要［20］，擴展后的液氫池尺寸和蒸發(fā)速率有著重要影響［19，21］。

因此，本文基于開源的CFD平臺OpenFOAM開發(fā)的求解器對液氫泄漏和氫氣云團在大氣中的擴散進行數(shù)值模擬。建立簡化的液氫蒸發(fā)模型，以便快速確定入口條件；采用RANS方法對云團的擴散進行模擬，并考慮了大氣穩(wěn)定度的影響以及地面?zhèn)鳠岬挠绊?。通過文獻中的兩組實驗［8，22］對開發(fā)的求解器進行驗證，并分析大氣穩(wěn)定度、環(huán)境溫度、泄漏量以及滯留坑對危險范圍的影響。

1 數(shù)學模型

1.1 氣相求解器

氫氣云團擴散過程涉及氣液兩相流動問題，并伴隨著相變現(xiàn)象［12］。鑒于相變通常發(fā)生在近場區(qū)域，對遠場擴散的影響較小，同時實驗結果表明氫氣與空氣混合引起的升溫在云團的浮力效應中占主導地位［8］。因此，本文目前忽略了相變和濕度因素，這些因素對浮力的影響將留待未來工作進行研究。采用Realizable k-ε模型進行湍流閉合?？紤]浮力效應的三維可壓縮N-S方程組：

[?ρ?t+▽·ρu=0] （1）

[?ρu?t+▽·ρuu=-Δp+▽·τeff+ρg] （2）

[?ρYi?t+▽·ρuYi=▽·μeffYi] （3）

[?ρhs?t+▽·ρuhs=▽·αeff▽hs+i=1n（▽·hsiρDi-α▽Yi）] （4）

式中：[ρ]——密度，kg/m3；[t]——時間，s；[u]——速度，m/s；[p]——壓力，Pa；[hs]——顯焓，J/kg；[D]——組分擴散系數(shù)，m2/s；[α]——熱擴散系數(shù)，m2/s；[μ]——動力黏度，kg/（m·s）；[τ]——應力張量，m2·s2；[g]——重力系數(shù)，N/kg；下標[i]代表第[i]個組分，eff代表有效值，具有標識符號“—”的變量為空間濾波變量。

1.2 固相求解器

氫氣云團溫度較低，來自地面的熱量將增加其浮力效應，因此對地面?zhèn)鳠徇M行求解。采用恒定的熱力屬性，同時忽略地面植被對傳熱的影響。三維控制方程為：

[ρsCp，s?T?t=▽·λs?T] （5）

式中：[λs]——固體導熱系數(shù)，W/（m·K）；[ρs]——固體密度，kg/m3；[Cp，s]——固體定壓比熱，J/（kg·K）。

1.3 共軛傳熱

通過氣固交界面處的邊界條件，保證交界面處的溫度和熱通量連續(xù)性。交界面邊界條件為：

[Tf，int=Ts，int] （6）

[Qcon=ks·?Ts/?y] （7）

式中：[Qcon]——對流傳熱熱流密度，W/m2；[y]——交界面法向量上的距離；[hc]——對流傳熱系數(shù)，W/（m2·K），采用適用于2～20 m/s風速的經(jīng)驗公式進行計算：

[hc=10.45-u+10u0.5] （8）

1.4 蒸發(fā)速率

液氫泄漏分為連續(xù)泄漏以及瞬時泄漏，本文針對危害性更大的瞬時泄漏進行研究?；谝韵录僭O，計算地面?zhèn)鬟f至液氫的熱流密度。1）地面的初始溫度為均勻分布；2）地面外表面溫度突然下降；3）距離地面內表面無限遠處位置的溫度不會發(fā)生變化。因此隨時間變化的熱流密度為：

[Hct=λsT0-Tb/as·π·tsp/2，tlt;tspλsT0-Tb/as·π·t，t≥tsp] （9）

式中：[Hct]——熱流密度，W/m2；[T0]——地面初始溫度，K；[Tb]——液氫沸點溫度，K；as——固體地面熱擴散系數(shù)，m2/s；[tsp]——延遲時間，s。為避免初始階段的預測蒸發(fā)速率過高，同時考慮到液氫并非瞬間實現(xiàn)與地面的完全接觸，應考慮液氫池的蔓延擴展過程，因此需設定延遲時間，其物理意義為液氫池擴展至最遠位置所需的時間，取為20 s［23］。

由于液氫吸收的熱量主要來自于地面的熱傳導，而輻射及空氣對流的熱量僅為地面導熱熱量的2.4%［1］，因此僅考慮地面的熱傳導作用，液氫蒸發(fā)速率為：

[qv=Hct/Lv] （10）

式中：[qv]——液氫蒸發(fā)速率，kg/（m·s）；[Lv]——蒸發(fā)潛熱，J/kg。

1.5 大氣邊界條件

氫氣云團的擴散過程發(fā)生在大氣邊界層內，因此高度以及大氣穩(wěn)定度對風廓線存在一定影響。采用對數(shù)型的風廓線模擬風場，并考慮大氣穩(wěn)定度對風速以及湍流的作用：

[U=u*lnzz0-φm/Ck] （11）

式中：[U]——平均風速，m/s；[Ck]——Karman常數(shù)，取0.4［24］；[z0]——地面粗糙度，m；[u*]——摩擦速度，m/s，[u*=U0Ck/[lnzref/z0-φm，ref]]，其中[U0]為參考高度處的風速，m/s；[φm]、[φm，ref]——穩(wěn)定度修正項和參考高度處的穩(wěn)定度修正項。

[φm=2ln1+r2+ln1+r22-2arctanr+π"不穩(wěn)定0，中立-171-exp-0.29zL，穩(wěn)定] （12）

式中：[L]——莫寧-奧布霍夫長度，m，[L=Ls/（lnz0/zs）]，其中[Ls]、[zs]為與大氣穩(wěn)定度相關的常數(shù)［24］；[r]——高度和[L]相關的函數(shù)，計算[r=（1-16z/L）0.25]。

不穩(wěn)定條件下湍流動能[ku]和湍流耗散率[εu]分別為：

[ku=0.36w2*+0.85u2*1-3zL23，z≤0.1habl0.36+0.9zhabl231-0.8zhabl2w2*， "zgt;0.1habl] （13）

[εu=u3*kz1+0.5zL2332，z≤0.1hablw3*habl0.8-0.3zhabl，zgt;0.1habl] （14）

式中：[ku]——不穩(wěn)定條件下的湍流動能，m2/s2；[εu]——不穩(wěn)定條件下的湍流耗散率，m2/s3；[w*]——地面熱流密度，W/m2；[habl]——邊界層高度［24］，m。

中性和穩(wěn)定條件下湍流動能和湍流耗散率分別為：

[kn，s=6u2*，z≤0.1habl6u2*1-zhabl74，zgt;0.1habl] （15）

[εn，s=u3*kz1.24+4.3zL，z≤0.1hablu3*kz1.24+4.3zL1-0.85zhabl32，zgt;0.1habl ] （16）

式中：[kn，s]——中性和穩(wěn)定條件下的湍流動能，m2/s2；εn，s——中性和穩(wěn)定條件下的湍流耗散率，m2/s3。

大氣邊界層高度［14］為：

[habl=1500，不穩(wěn)定0.3u*/fc，中立0.4（u*L/fc）1"穩(wěn)定] （17）

式中：[fc]——科里奧利參數(shù)。

1.6 數(shù)值方法

基于開源CFD平臺OpenFOAM，根據(jù)上述模型方法搭建內部求解器“chtAtmDisFOAM”。該求解器中采用PIMPLE算法處理壓力和速度之間的相互耦合。時間項采用一階精度的Euler離散格式，動量方程以及標量方程中的對流項分別采用二階精度的LUST格式以及二階精度的TVD格式，擴散項采用二階精度的高斯離散格式。

2 問題描述

通過中國科學院大氣物理研究所監(jiān)測的不同大氣穩(wěn)定度下平均風廓線監(jiān)測數(shù)據(jù)［22］以及NASA第6次實驗數(shù)據(jù)［8］對開發(fā)的求解器進行驗證。前者測量了錫林郭勒草原中心位置的不同大氣穩(wěn)定度風廓線特征，其地面條件為平坦草原，10 m參考高度的平均風速約為6 m/s，詳細條件見文獻［22］。后者對5.7 m3的液氫進行釋放，形成的液氫池直徑約為5 m，環(huán)境溫度約為288 K，環(huán)境壓力為標準大氣壓。液氫溫度為20.35 K，儲存于690 kPa的杜瓦瓶中。10 m高處的參考風速為2.2 m/s，大氣穩(wěn)定度接近于穩(wěn)定條件，地面為濕潤的沙土，導熱系數(shù)［23］為0.25 W/（m·K），比熱容［23］為800 J/（kg·K），密度［23］為1800 kg/m3。該實驗屬于儲存在高于大氣壓的泄漏情形，發(fā)生閃蒸的液氫占比較低，因此忽略閃蒸產(chǎn)生的影響。

對液氫泄漏源進行一定簡化：假設液氫瞬間釋放，在地面形成的液池迅速達到最小厚度以及最大直徑。液氫池最小厚度取為地表的表面粗糙度［23］。因此可通過泄漏液氫的總體積及液池厚度確定液池面積，并使用提出的簡化蒸發(fā)模型確定液池蒸發(fā)速率。形成的液氫池通常接近圓形，為便于網(wǎng)格劃分假設形成等面積的方形液池。

計算區(qū)域如圖1所示，包括氣相開放空間以及固體地面區(qū)域，長寬高分別為200 m/100 m/100 m和200 m/100 m/10 m。液氫池設置于風入口邊界下游40 m處。計算網(wǎng)格通過前處理工具blockMesh生成。氫氣入口設置16×16個網(wǎng)格，并在其附近以及近場區(qū)域進行局部加密：氣相區(qū)域高度方向設置80層網(wǎng)格，由下往上的網(wǎng)格膨脹率為1."長度方向總共設置170層網(wǎng)格，由氫氣入口至下游邊界的網(wǎng)格膨脹率為1."寬度方向總共設置60個網(wǎng)格，氫氣入口至兩側的網(wǎng)格膨脹率為1.2；固相區(qū)域高度方向設置40層網(wǎng)格，由上往下的網(wǎng)格膨脹率為1.2。氣相區(qū)域和固體區(qū)域網(wǎng)格總數(shù)分別為816000和408000。

3 驗 證

由于計算區(qū)域較大，需采用合理的網(wǎng)格尺寸以同時滿足計算效率以及計算分辨率要求。對上文所設置的計算網(wǎng)格進行整體加密及加粗（保持網(wǎng)格膨脹率不變，增大/減小每個方向上的網(wǎng)格數(shù)量）。氫氣入口處的網(wǎng)格數(shù)量分別為24×24個以及10×10個。加密及加粗后的氣相區(qū)域網(wǎng)格數(shù)量約為120萬和40萬，固體區(qū)域網(wǎng)格數(shù)量約為60萬和20萬。圖2為3種網(wǎng)格分辨率下預測可視云團擴散距離，其中中等網(wǎng)格和細網(wǎng)格預測結果的相對誤差不超過4%。因此后續(xù)計算采用中等網(wǎng)格進行計算。

通過中國科學院大氣物理研究所監(jiān)測的不同大氣穩(wěn)定度下平均風廓線監(jiān)測數(shù)據(jù)［22］對3種大氣穩(wěn)定度下的風廓線進行驗證。圖3為不同穩(wěn)定度條件下預測風廓線與測量數(shù)據(jù)［22］的對比。隨著高度的增加，不同穩(wěn)定度對應的歸一化風速差異逐漸增大，穩(wěn)定條件下歸一化風速明顯高于其他兩者。

在3種穩(wěn)定度條件下，預測結果在100 m以下的下層表現(xiàn)良好，上層風速略微高估，總體可認為使用的大氣邊界條件滿足計算要求。

數(shù)值模擬捕捉20.94 s的氫氣濃度分布云圖如圖4所示。從圖4可發(fā)現(xiàn)，此時氫氣云團處于重氣擴散階段，并向下游位置繼續(xù)擴散，且尚未完全脫離地面。圖4中預測云團濃度分布、運動趨勢與實驗結果［8］大致相符，且云團傾角相似。值得注意的是，預測結果中可燃氫氣云團前端與地面發(fā)生分離的距離與實驗存在一定差距，該誤差原因可能是此時云團仍處于近場，忽略相變所導致。

采樣點處的預測及實驗瞬時氫氣濃度［8］對比如圖5所示。采樣點位為監(jiān)測塔2高度為1 m處以及監(jiān)測塔5高度為9.4 m處。圖5中，實驗監(jiān)測氫氣濃度存在明顯的波動，預測結果同樣捕獲到了這些波動，而劉元亮等［10］的預測結果中，氫氣濃度變化較為平穩(wěn)，未捕捉到氫氣濃度的波動。對比結果表明，本文提出的求解器能較為準確的預測近場氫氣濃度的瞬時變化情況。

圖6為監(jiān)測塔7、8、9處的預測瞬時氫氣濃度與實驗監(jiān)測最大氫氣濃度［8］對比。實驗中高度為1、9.4和18.4 m處的最大氫氣濃度分別為0%、18.7%以及19%，而對應的預測最大值分別為0.34%、19.6%以及21.5%，與前者十分接近。這表明開發(fā)的求解器能準確預估近場局部位置的氫氣最大濃度，可用于判斷氫氣濃度是否滿足安全要求。

圖7為預測可視云團（氫氣濃度8%）擴散最大高度、最遠距離以及消失時間與實驗觀測值［8］的對比。預測可視氫氣云團的最大高度和最遠距離均與實驗數(shù)據(jù)接近，相對誤差分別為3.5%和5.6%。預測消失時間為82.5 s，接近實驗值90 s。整體上看，所提數(shù)值模型可用于預測液氫泄漏后大氣環(huán)境下的云團擴散行為。

4 影響因素研究

氫氣云團的擴散行為受眾多因素影響。本文研究大氣穩(wěn)定度、環(huán)境溫度以及泄漏量等因素對氫氣云團擴散行為的影響，并分析泄漏點處設置滯留坑方案的可行性。若無明確說明，算例參考風速均為2.2 m/s，其他計算條件如表1所示。

當儲罐整體破裂時，大量液氫會在瞬間以極高的速率蒸發(fā)。圖8為不同條件下液氫蒸發(fā)速率。圖8中，泄漏量為1 t時，隨著環(huán)境溫度的升高，蒸發(fā)速率有所增加，而蒸發(fā)所需時間減少。當泄漏量分別達到5 t和10 t時，蒸發(fā)速率達200和401 kg/s，并均在30 s內完成蒸發(fā)。無滯留坑時，10 t液氫形成的液氫池尺寸達78.81 m。而當滯留坑尺寸取為8 m時，滯留坑的存在將阻礙液氫池的蔓延擴展，將其限制在極小區(qū)域內，減小了液氫與地面的接觸面積，進而減小液氫接收到的熱量，導致蒸發(fā)速率減小達99%。此外，蒸發(fā)所需時間急劇增加，如圖8中局部放大圖所示。

圖9為不同大氣穩(wěn)定度條件下預測氫氣可燃云團（氫氣濃度4%）的擴散范圍。由于初始時刻處于重氣擴散階段，云團尚未脫離地面，不同穩(wěn)定度條件下預測可燃云團高度幾乎一致。而此時擴散距離受風速影響，穩(wěn)定條件下可燃云團向下游擴散的速度明顯高于其他兩種情況。隨著時間的推移，云團脫離地面，進入浮升擴散階段，由于高度不斷上升3種大氣穩(wěn)定度對應的風速差異逐漸增大，此時穩(wěn)定條件下的最大擴散距離明顯大于其他兩者。與之相反的是，由于高風速對云團在高度方向的擴散存在抑制作用，穩(wěn)定條件下的云團最大擴散高度最低。此外，由于強不穩(wěn)定條件下的高湍流作用，其可燃云團最早消失。

表2為不同環(huán)境溫度、不同液氫泄漏量以及不同地面邊界條件對可燃云團擴散范圍的影響。預測結果中，隨著環(huán)境溫度升高，浮力作用增強，可燃云團的擴散高度以及擴散距離均隨之增大：環(huán)境溫度增加40 K后，可燃云團的擴散距離和高度分別增大6.2%及3.8%。相比于環(huán)境溫度，泄漏量對可燃云團的影響則更為顯著。大規(guī)模泄漏產(chǎn)生的液氫池尺寸巨大，并伴隨著極大的蒸發(fā)速率，導致隨著大量的可燃氫氣聚集，因此所需的稀釋時間以及擴散范圍均隨之增大，泄漏量增大10倍后，可燃云團的擴散距離和高度分別增大93.9%及216.8%。

此外，表2列出了氣固交界面處使用絕熱條件（AC）以及共軛傳熱方法（CHT）所對應的可視云團擴散范圍。計算中考慮了風速較低情形和風速較高情形兩種風速條件。在兩種參考風速下，考慮地面?zhèn)鳠岷徒^熱地面條件獲得的預測結果幾乎一致，僅在高風速時，氫氣云團由于在垂直方向被抑制而更接近于地面，此情況下地面?zhèn)鳠釋U散距離和高度產(chǎn)生了略微影響。

圖10為有滯留坑和無滯留坑情況下可燃云團預測擴散范圍。需注意的是，計算中滯留坑幾何形狀的影響被忽略。有滯留坑情況下，地面?zhèn)飨蛞簹涑氐臒崃看罅拷档?，導致氫氣云團浮力減小，最大擴散高度降低了76%，同時最大擴散距離也減小45.9%。由于蒸發(fā)速率的差異，無滯留坑情形下可燃云團在短時間內（147 s）消失殆盡，而在有滯留坑情形下可燃云團長時間維持在一小范圍區(qū)域。

5 結 論

本文基于開源CFD平臺OpenFOAM開發(fā)一內部求解器，用于模擬大規(guī)模液氫泄漏以及形成的氫氣云團擴散行為。該求解器通過簡化模型確定液氫的蒸發(fā)速率，并結合考慮大氣穩(wěn)定度的大氣邊界層條件以及地面?zhèn)鳠帷Ｍㄟ^中國科學院大氣物理研究所測量的不同穩(wěn)定度下平均風廓線以及NASA的大規(guī)模液氫泄漏實驗數(shù)據(jù)對開發(fā)的求解器進行了驗證。3種穩(wěn)定度下風廓線的預測結果與監(jiān)測數(shù)據(jù)非常吻合，液氫泄漏后的云團形態(tài)、近場處的氫氣瞬時濃度、最大濃度和遠場處的可視氫氣云團最大擴散距離、高度及消失時間均與實驗結果較為符合，證明了求解器的準確性和有效性。

研究發(fā)現(xiàn)，液氫泄漏后的云團擴散行為受大氣穩(wěn)定度、環(huán)境溫度和泄漏量影響：大氣邊界層越穩(wěn)定，氫氣云團的橫向移動范圍越大；強不穩(wěn)定條件下氫氣云團更快被稀釋；環(huán)境溫度對可燃云團安全范圍的影響較小；相比之下泄漏量產(chǎn)生的影響更為顯著；此外采用絕熱地面邊界條件和求解地面共軛傳熱對可視云團擴散范圍幾乎無影響；滯留坑在降低液氫蒸發(fā)速率方面有顯著效果，且使可燃云團擴散范圍更小但更持久。若無滯留坑，可燃云團將擴散至更遠更高的區(qū)域，但將很快消失。

［參考文獻］

［1］ VERFONDERN K， DIENHART B. Experimental and theoretical investigation of liquid hydrogen pool spreading and vaporization［J］. International journal of hydrogen energy， 1997， 22（7）： 649-660.

［2］ HECHT E S， PANDA P P. Mixing and warming of cryogenic hydrogen releases［J］. International journal of hydrogen energy， 2019， 44（17）： 8960-8970.

［3］ KOBAYASHI H， NARUO Y， MARU Y， et al. Experiment of cryo-compressed （90-MPa） hydrogen leakage diffusion［J］. International journal of hydrogen energy， 2018， 43（37）： 17928-17937.

［4］ KOBAYASHI H， DAIMON Y， UMEMURA Y， et al. Temperature measurement and flow visualization of cryo-compressed hydrogen released into the atmosphere［J］. International journal of hydrogen energy， 2018， 43（37）： 17938-17953.

［5］ SHU Z Y， LEI G， LIANG W Q， et al. Experimental investigation of hydrogen dispersion characteristics with liquid helium spills in moist air［J］. Process safety and environmental protection， 202""162： 923-931.

［6］ STATHARAS J C， VENETSANOS A G， BARTZIS J G， et al. Analysis of data from spilling experiments performed with liquid hydrogen［J］. Journal of hazardous materials， 2000， 77（1-3）： 57-75.

［7］ HEDLEY D， HAWKSWORTH S J， RATTIGAN W， et al. Large scale passive ventilation trials of hydrogen［J］. International journal of hydrogen energy， 2014， 39（35）： 20325-20330.

［8］ WITCOFSKI R D， CHIRIVELLA J E. Experimental and analytical analyses of the mechanisms governing the dispersion of flammable clouds formed by liquid hydrogen spills［J］. International journal of hydrogen energy， 1984， 9（5）： 425-435.

［9］ SUN R F， PU L， HE Y C， et al. Phase change modeling of air at the liquid hydrogen release［J］. International journal of hydrogen energy， 2024， 50： 717-731.

［10］ LIU Y L， WEI J J， LEI G， et al. Dilution of hazardous vapor cloud in liquid hydrogen spill process under different source conditions［J］. International journal of hydrogen energy， 2018， 43（15）： 7643-7651.

［11］ PU L， TANG X， SHAO X Y， et al. Numerical investigation on the difference of dispersion behavior between ""cryogenic ""liquid ""hydrogen ""and "methane［J］. International journal of hydrogen energy， 2019， 44（39）： 22368-22379.

［12］ GIANNISSI S G， VENETSANOS A G. Study of key parameters in modeling liquid hydrogen release and dispersion in open environment［J］. International journal of hydrogen energy， 2018， 43（1）： 455-467.

［13］ LIU Y L， WEI J J， LEI G， et al. Modeling the development of hydrogen vapor cloud considering the presence "of "air "humidity［J］. International "journal "of hydrogen energy， 2019， 44（3）： 2059-2068.

［14］ GIANNISSI S G， VENETSANOS A G， MARKATOS N， et al. CFD modeling of hydrogen dispersion under cryogenic release conditions［J］. International journal of hydrogen energy， 2014， 39（28）： 15851-15863.

［15］ TANG X， PU L， SHAO X Y， et al. Dispersion behavior and safety study of liquid hydrogen leakage under different application situations［J］. International journal of hydrogen energy， 2020， 45（55）： 31278-31288.

［16］ HANSEN O R. Liquid hydrogen releases show dense gas behavior［J］. International journal of hydrogen energy， 2020， 45（2）： 1343-1358.

［17］ SHAO X Y， PU L， TANG X， et al. Parametric influence study of cryogenic hydrogen dispersion on theoretical aspect［J］. International journal of hydrogen energy， 2020， 45（38）： 20153-20162.

［18］ 孫壯， 李樹民， 朱蓉， 等. 考慮大氣穩(wěn)定度的典型復雜地形風場CFD模式［J］. 太陽能學報， 2024， 45（4）： 197-205.

SUN Z， LI S M， ZHU R， et al. CFD model of wind field in typical complex terrain considering atmospheric stability［J］. Acta energiae solaris sinica， 2024， 45（4）： 197-205.

［19］ 楊祥生， 趙寧， 田琳琳， 等. 基于Park-Gauss模型考慮大氣穩(wěn)定性對布局優(yōu)化的影響［J］. 太陽能學報， 202""42（3）： 43-47.

YANG X S， ZHAO N， TIAN L L， et al. Influence of atmospheric stability on layout optimization based on park-Gauss model［J］. Acta energiae solaris sinica， 202""42（3）： 43-47.

［20］ MIDDHA P， HANSEN O R， STORVIK I E. Validation of CFD-model for hydrogen dispersion［J］. Journal of loss prevention in the process industries， 2009， 22（6）： 1034-1038.

［21］ VERFONDERN K， DIENHART B. Pool spreading and vaporization of liquid hydrogen［J］. International journal of hydrogen energy， 2007， 32（2）： 256-267.

［22］ 賀園園， 程雪玲， 朱蓉. 適用于平坦草原的近地層以上風廓線推算方法［J］. 太陽能學報， 2024， 45（6）： 451-460.

HE Y Y， CHENG X L， ZHU R. Method of wind speed profile extrapolating from surface observations over flat grassland［J］. Acta energiae solaris sinica， 2024， 45（6）： 451-460.

［23］ BOSCH C， WETERINGS R. Methods for the calculation of physical effects［M］. Hague： Min. VROM， 1997.

［24］ HAN J， ARYA S P， SHAOHUA S， et al. An estimation of turbulent kinetic energy and energy dissipation rate based on atmospheric boundary layer similarity theory［R］. NAS 1.26：210298， 2000.

NUMERICAL SIMULATION OF LIQUID HYDROGEN LEAKAGE AND DISPERSION CONSIDERING ATMOSPHERIC STABILITY

Cheng Chonglyu， Shan Conghui， Zhang Mengfan， Xu Baopeng

（School of Energy and Power， Dalian University of Technology， Dalian 116024， China）

Abstract：The flammable and explosive cloud resulting from a liquid hydrogen leak poses a serious safety hazard. Using a solver developed with OpenFOAM， an in-depth analysis of liquid hydrogen leakage and subsequent dispersion behavior was conducted. The solver employed a simplified evaporation model to determine the evaporation rate. Hydrogen cloud dispersion and heat transfer from the ground were simulated using separate gas-phase and solid-phase solvers. The conjugate heat transfer method was used to ensure continuous temperature and heat flux between the gas and solid regions. A logarithmic wind profile incorporating atmospheric stability was also utilized. The solver was validated using wind profile measurements and liquid hydrogen leakage experiments. The results showed stable atmospheric conditions， high temperatures， and large leakage volumes increased lateral dispersion of flammable hydrogen clouds. Additionally， a retention pit scheme reduced the evaporation rate by 99%， flammable cloud distance by 46%， and cloud height by 76%， effectively decreasing the hazardous range.

Keywords：hydrogen； leakage （fluid）； atmospheric boundary layer； dispersion； numerical simulation