李昌均 蔡幸福 肖明清 王 飛 曹曉巖

1（火箭軍工程大學(xué) 西安 710025）

2（空軍工程大學(xué) 西安 710051）

氚是氫的一種放射性同位素，通常以氣體氚和氚水的形式存在，具有極強(qiáng)的滲透和擴(kuò)散能力，所以氚大多被貯存在高壓儲(chǔ)存容器內(nèi)。在意外條件作用下，高壓容器可能會(huì)出現(xiàn)破損，導(dǎo)致氚泄漏，瞬間釋放并擴(kuò)散到環(huán)境中，對(duì)人員和環(huán)境造成危害。因此，氚泄漏是典型的核事故場(chǎng)景，開展氚泄漏數(shù)值模擬研究對(duì)核事故應(yīng)急決策與處置具有重要的理論意義和現(xiàn)實(shí)價(jià)值。

目前，在高壓氣體泄漏研究及氣體狀態(tài)方程應(yīng)用方面，趙芳等［1］利用理想氣體模型和真實(shí)氣體模型開展了高壓氣體定容積定流量放氣特性研究，結(jié)果表明：高壓放氣下，理想氣體模型存在偏差，真實(shí)氣體模型更接近于實(shí)驗(yàn)值；王振華等［2］利用Abel-Noble氣體狀態(tài)方程模擬了高壓氫氣泄漏過(guò)程，并建立了氫氣泄漏噴射火預(yù)測(cè)模型，通過(guò)實(shí)驗(yàn)數(shù)據(jù)比對(duì)，驗(yàn)證Abel-Noble狀態(tài)方程模擬高壓氣體泄漏的準(zhǔn)確性；王建吉等［3］建立了基于Abel-Noble 的狀態(tài)方程模型，用于計(jì)算和預(yù)測(cè)變截面渦旋壓縮機(jī)徑向泄漏量的大小，拓展了Abel-Noble 狀態(tài)方程模型的應(yīng)用領(lǐng)域；魯寨軍等［4-6］利用UDF（User defined functions）分別將Span-Wagner 狀態(tài)方程和P-R 狀態(tài)方程嵌入Fluent求解器中，模擬計(jì)算了高壓CO2氣體泄漏瞬態(tài)行為，實(shí)現(xiàn)了氣體狀態(tài)方程與Fluent耦合求解；在氚泄漏研究方面，Iwai等［7-9］利用流體動(dòng)力學(xué)模型對(duì)室內(nèi)環(huán)境氚泄漏行為進(jìn)行了研究，得到的氚氣濃度模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。劉貴議等［10］針對(duì)產(chǎn)氚模塊氚提取系統(tǒng)（Tritium Extraction System of Tritium Breeder Module，TBM TES）氚泄漏事故，模擬計(jì)算了氚濃度在不同換氣情況下達(dá)到安全值的條件；Li等［11］對(duì)低壓儲(chǔ)氚容器在密閉、開放空間的氚泄漏事故，研究了氚氣的擴(kuò)散行為。

盡管前人對(duì)氚泄漏事故已經(jīng)開展了一些研究，但研究對(duì)象大多是低壓氣體連續(xù)泄漏，且事故場(chǎng)景各不相同，對(duì)于高壓儲(chǔ)氚容器在特定場(chǎng)景空間發(fā)生泄漏問(wèn)題，缺少對(duì)高壓氚泄漏的瞬態(tài)行為分析，以及泄漏后氚氣在不同工況下的時(shí)空分布規(guī)律。

綜上所述，針對(duì)高壓儲(chǔ)氚容器氚泄漏問(wèn)題，提出基于AN-EOS（Abel-Noble Equation of State）模型的高壓氚泄漏數(shù)值模擬方法，建立基于AN-EOS 的氚泄漏瞬態(tài)求解模型，定義UDF 自定義函數(shù)，實(shí)現(xiàn)AN-EOS與Fluent求解器的耦合計(jì)算，利用Fluent軟件［12］計(jì)算氚氣流動(dòng)，得到氚氣濃度的時(shí)空分布規(guī)律。

1 氚泄漏模擬方法

1.1 AN-EOS模型

目前，描述氣體泄漏的熱力學(xué)模型主要有兩類：理想氣體模型和實(shí)際氣體狀態(tài)方程模型［13］，理想氣體模型假定氣體壓縮因子為1，實(shí)際氣體狀態(tài)方程是經(jīng)過(guò)大量實(shí)驗(yàn)對(duì)理想氣體狀態(tài)方程修正后的半經(jīng)驗(yàn)?zāi)Ｐ?。AN-EOS［14］模型是在理想氣體模型基礎(chǔ)上，考慮了氣體分子的體積帶來(lái)的影響，對(duì)氣體壓縮因子修正后得到的氣體狀態(tài)方程。在高壓儲(chǔ)存條件下，氣體分子行為偏離理想氣體行為，AN-EOS模型比理想氣體模型更能準(zhǔn)確地描述氣體泄漏的物理特性［15］。因此，本文選用AN-EOS 模型開展高壓氣體泄漏的瞬態(tài)行為研究，AN-EOS具體形式為：

式中：p為壓強(qiáng)；v為比體積；b為氣體分子比積體修正 項(xiàng)；b=7.691×10-3m3·kg-1；Rg為 氣 體 常 數(shù)；T為溫度。

假定氣體的泄漏過(guò)程為等熵過(guò)程，由AN-EOS模型可得：

式中：k為比熱容比；constant為常量。

定容積儲(chǔ)氣罐高壓氣體泄漏過(guò)程主要分為兩個(gè)階段：音速流動(dòng)階段和亞音速流動(dòng)階段。容器內(nèi)壓強(qiáng)與環(huán)境壓強(qiáng)之比大于臨界壓強(qiáng)比時(shí)，泄漏口處氣流速度達(dá)到當(dāng)?shù)芈曀?，即?/p>

其中：

式中：vcr為臨界壓強(qiáng)比；cf為當(dāng)?shù)芈曀?；v2為泄漏口處比體積；T2為泄漏口處溫度；p2為泄漏口壓強(qiáng)；pi為儲(chǔ)氣罐內(nèi)壓強(qiáng)。

泄漏口處壓強(qiáng)p2=pi?vcr，此時(shí)，泄漏口處的質(zhì)量流率為：

式中：Qm為泄漏口質(zhì)量流率；A為泄漏口面積。

容器內(nèi)壓強(qiáng)與環(huán)境壓強(qiáng)（pamb）之比小于或等于臨界壓強(qiáng)比時(shí)，泄漏口氣流為亞音速流動(dòng)，泄漏口壓強(qiáng)減小到環(huán)境壓強(qiáng)，即p2=pamb，此后保持不變，泄漏口處氣流速度為：

此時(shí)，泄漏口處的質(zhì)量流率為：

由于儲(chǔ)氣罐容積不變，則有d（V）=0，即：

對(duì)式（2）求微分后聯(lián)立式（7）可得：

經(jīng)計(jì)算可知，該微分方程難以求得理論解，因此采用迭代法求解方程，則：

式中：Q′m為亞音速階段泄漏口質(zhì)量流率；V為儲(chǔ)氣罐體積；mi為罐內(nèi)氣體質(zhì)量；vi為罐內(nèi)比體積；Ti為罐內(nèi)溫度；j為迭代步數(shù)；Δt為時(shí)間步長(zhǎng)。

1.2 模型構(gòu)建

1.2.1模型假設(shè)

本文研究的對(duì)象為高壓儲(chǔ)氚容器，內(nèi)部壓強(qiáng)為34.5 MPa，溫度為300 K，忽略儲(chǔ)氚容器內(nèi)氚氣質(zhì)量分?jǐn)?shù)因衰變發(fā)生的變化，對(duì)模型作出以下假設(shè)：

1）氚泄漏事故場(chǎng)景，假設(shè)氚泄漏事故發(fā)生在一個(gè)長(zhǎng)方體空間，長(zhǎng)寬高分別為5.0 m、3.0 m和2.5 m，在正面墻體上方設(shè)置一個(gè)進(jìn)風(fēng)口，右側(cè)墻體上方設(shè)置一個(gè)排風(fēng)口，兩者均為0.5 m×0.5 m 的正方形（圖1）。

2）高壓儲(chǔ)氚容器簡(jiǎn)化為一個(gè)立方體，邊長(zhǎng)為0.3 m，放置在場(chǎng)景地面中心位置（圖1）；泄漏孔簡(jiǎn)化為一個(gè)圓形孔口，直徑為10 mm，且泄漏中孔口直徑不變。

圖1 氚泄漏場(chǎng)景物理模型Fig.1 Physical model of tritium leakage scene

3）監(jiān)測(cè)點(diǎn)位設(shè)置，為了準(zhǔn)確描述氚氣在泄漏中的濃度變化，在密閉空間設(shè)置8 個(gè)監(jiān)測(cè)點(diǎn)位P1~P8，假設(shè)操作人員平均身高1.7 m，P1~P5主要監(jiān)測(cè)氚氣在空間內(nèi)z=1.7 m 的濃度分布情況，P1 位于泄漏源上方，P2~P5位于4個(gè)角落，其中P6、P7為進(jìn)風(fēng)口、排風(fēng)口處監(jiān)測(cè)點(diǎn)，P8 為空間頂部監(jiān)測(cè)點(diǎn)，各監(jiān)測(cè)點(diǎn)坐標(biāo)如表1所示。

表1 監(jiān)視點(diǎn)坐標(biāo)Table 1 Monitor point coordinates

1.2.2網(wǎng)格劃分

整個(gè)物理模型采用六面體非結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分，為精細(xì)化描述泄漏口和邊界區(qū)域的氣體流動(dòng)情況，對(duì)泄漏口進(jìn)行加密網(wǎng)格處理，對(duì)壁面添加邊界層，流體區(qū)域根據(jù)模型尺寸自動(dòng)設(shè)置網(wǎng)格大小，模型總網(wǎng)格數(shù)313 303個(gè)，節(jié)點(diǎn)數(shù)1 275 844個(gè)（圖2）。網(wǎng)格檢查結(jié)果顯示：模型中最小尺寸為4.24×10-11m3，最大尺寸為2.8×10-4m3，無(wú)負(fù)網(wǎng)格情況；網(wǎng)格質(zhì)量檢查顯示：最小正交質(zhì)量為0.58，滿足網(wǎng)格質(zhì)量要求，可作為Fluent求解的模型輸入。

圖2 劃分的網(wǎng)格模型 (a)整體網(wǎng)格（含x軸方向切割面），(b)局部網(wǎng)格Fig.2 Partition of the grid model (a)Overall grid(including an x-axis cutting surface),(b)Partial grid

1.2.3Fluent求解器設(shè)置

在模型求解設(shè)置中，氚氣在空氣中的流動(dòng)屬于多組分傳輸問(wèn)題，打開組分傳輸模型；高壓儲(chǔ)氚容器內(nèi)的壓強(qiáng)高達(dá)數(shù)百個(gè)大氣壓，當(dāng)發(fā)生氚泄漏時(shí)，氚氣初始速度較大，在密閉空間中的流動(dòng)包括射流和擴(kuò)散，因此，采用能適用于復(fù)雜流動(dòng)類型的Realizablek-ε模型［16］。

在邊界類型設(shè)置中，將氚氣泄漏源設(shè)置為質(zhì)量入口，通過(guò)UDF 自定義函數(shù)定義入口質(zhì)量流率、泄漏口壓強(qiáng)和溫度隨時(shí)間的變化關(guān)系。將進(jìn)風(fēng)口定義為速度入口（Velocity inlet），換氣系統(tǒng)對(duì)密閉空間換氣頻率為每小時(shí)4 次，即平均風(fēng)速為1.5 m?s-1，將排風(fēng)口定義為速度入口，平均風(fēng)速為-1.5 m?s-1；密閉條件下將進(jìn)風(fēng)口和排風(fēng)口設(shè)置為墻面。在材料設(shè)置中，添加氚氣材料［17-18］，設(shè)置相關(guān)參數(shù)（表2），將密閉空間壁面材料設(shè)置為碳酸鈣，將儲(chǔ)氚容器材料設(shè)置為鋼鐵，求解流場(chǎng)的數(shù)值方法選擇基于壓力耦合方程組的半隱式方法——SIMPLE（Semi-Implicit Method for Pressure Linked Equations）算法。

表2 氚氣物理參數(shù)Table 2 Parameters of tritium

2 結(jié)果與討論

本文首先利用真實(shí)氣體泄漏數(shù)據(jù)［19］對(duì)模型進(jìn)行驗(yàn)證，其次利用UDF自定義函數(shù)定義氚氣泄漏瞬態(tài)過(guò)程，開展密閉和通風(fēng)條件下高壓儲(chǔ)氚容器氚泄漏數(shù)值模擬。由AN-EOS模型計(jì)算可知，當(dāng)t=4.7 s時(shí)，氚氣停止泄漏，如圖3 所示。因此，F(xiàn)luent 求解器設(shè)置模擬時(shí)間應(yīng)大于4.7 s。

圖3 氚氣泄漏質(zhì)量流率變化規(guī)律Fig.3 Variation of mass flow rate of tritium leakage with time

2.1 模型驗(yàn)證

為了驗(yàn)證AN-EOS 高壓氣體泄漏模型的準(zhǔn)確性，將該模型模擬結(jié)果與真實(shí)氣體泄漏數(shù)據(jù)比對(duì)分析，設(shè)定與文獻(xiàn)［19］相同的實(shí)驗(yàn)條件：泄漏氣體壓強(qiáng)為34.5 MPa，氣體溫度300 K，泄漏孔面積為3.17×10-5m3，儲(chǔ)氣罐容積為2.73×10-2m3，環(huán)境壓強(qiáng)為標(biāo)準(zhǔn)大氣壓。將模型參數(shù)按照實(shí)驗(yàn)條件設(shè)置求解，對(duì)比分析模擬結(jié)果與參考數(shù)據(jù)泄漏口壓強(qiáng)和質(zhì)量流率的變化情況（圖4）。

圖4 模型驗(yàn)證 (a)泄漏口壓強(qiáng)變化，(b)質(zhì)量流率變化Fig.4 Model validation (a)Pressure changes at the leakage outlet,(b)Mass flow rate changes

由圖4可以看出，AN-EOS模型與真實(shí)氣體在泄漏口壓強(qiáng)和質(zhì)量流率的變化趨勢(shì)一致。表3給出了模擬結(jié)果與參考數(shù)據(jù)的誤差分析，在氣體泄漏初始時(shí)刻，AN-EOS 模型的泄漏口壓強(qiáng)與質(zhì)量流率均比真實(shí)氣體大，泄漏口壓強(qiáng)最大相對(duì)誤差為17.5%，最小為6.9%，質(zhì)量流率最大相對(duì)誤差為12.1%，最小為1.1%，隨著氣體泄漏的繼續(xù)進(jìn)行，整體誤差逐漸減小，由此可知，本文建立的AN-EOS模型可用于開展高壓氣體泄漏模擬實(shí)驗(yàn)。

表3 誤差分析Table 3 Error analysis

2.2 氚泄漏瞬態(tài)行為

為準(zhǔn)確描述氚泄漏瞬態(tài)行為，對(duì)氚泄漏開展熱力學(xué)模型計(jì)算，得到儲(chǔ)氚容器內(nèi)氚氣質(zhì)量和氚泄漏速度隨時(shí)間變化規(guī)律（圖5），選取通風(fēng)條件下泄漏后0.6 s、2 s和4 s時(shí)刻開展了模擬實(shí)驗(yàn)，得到x=1.5 m處YZ平面氚氣流動(dòng)速度矢量圖（圖6）。

由圖5 可以看出，初始時(shí)刻儲(chǔ)氚容器內(nèi)氚氣質(zhì)量約為1.37 kg，氚氣泄漏速度約為859.23 m·s-1，隨著氚泄漏的持續(xù)進(jìn)行，氚氣質(zhì)量隨時(shí)間非線性變化，氚氣泄漏速度在音速流動(dòng)階段（0~3.6 s）隨時(shí)間非線性變化，在亞音速流動(dòng)階段（3.6~4.7 s）隨時(shí)間近似線性變化；當(dāng)t=2.3 s 時(shí)，氚氣質(zhì)量為0.137 kg，儲(chǔ)氚容器泄漏量約90%，當(dāng)t=4.7 s 時(shí)，氚氣停止泄漏，容器內(nèi)氚氣質(zhì)量約為0.036 kg。由此可知，為保持儲(chǔ)氚容器與外界環(huán)境動(dòng)態(tài)平衡，容器內(nèi)會(huì)保留一定量的氚氣。

圖5 氚泄漏瞬態(tài)變化 (a)質(zhì)量變化，(b)速度變化Fig.5 Transient change of tritium leakage (a)Mass change,(b)Velocity change

由圖6可以看出，儲(chǔ)氚容器發(fā)生泄漏時(shí)，氚氣流動(dòng)初始速度較大，使得氚氣的流動(dòng)方向主要沿射流中心線，且氚氣在極短時(shí)間內(nèi)到達(dá)空間頂部，由于受到封閉空間的影響，使得氣流在射流中心線兩側(cè)形成渦流，如圖6（a）所示；通風(fēng)條件下，隨著氚泄漏的持續(xù)進(jìn)行，儲(chǔ)氚容器內(nèi)壓強(qiáng)迅速降低，氚泄漏速度逐漸減小，氚氣的流動(dòng)受空氣流動(dòng)影響逐漸偏向排風(fēng)口方向，如圖6（b）和（c）所示。

圖6 x=1.5 m處YZ平面的速度矢量圖 (a)t=0.6 s，(b)t=2 s，(c)t=4 sFig.6 Velocity vector profiles in YZ plane at x=1.5 m (a)t=0.6 s,(b)t=2 s,(c)t=4 s

2.3 氚氣濃度分布

為描述密閉和通風(fēng)環(huán)境下氚氣濃度在事故場(chǎng)景空間的流動(dòng)分布規(guī)律，采用場(chǎng)景空間的8 個(gè)監(jiān)測(cè)點(diǎn)位，分別測(cè)量不同時(shí)刻監(jiān)測(cè)的氚氣濃度（圖7、8），選取通風(fēng)條件下泄漏后21 s 開展模擬仿真實(shí)驗(yàn)，得到x=2.0 m處YZ平面氚氣濃度分布圖（圖9）。

圖7為密閉和通風(fēng)環(huán)境下P2~P5 監(jiān)測(cè)點(diǎn)氚氣濃度隨時(shí)間的變化規(guī)律，當(dāng)儲(chǔ)氚容器發(fā)生泄漏時(shí)，即0~4.7 s 內(nèi)，空間各監(jiān)測(cè)點(diǎn)位氚氣濃度迅速增加，其中，密閉環(huán)境下氚氣濃度穩(wěn)定增加，通風(fēng)環(huán)境下氚氣濃度的增加受空氣流動(dòng)影響出現(xiàn)波動(dòng)。當(dāng)停止泄漏后（t＞4.7 s），密閉條件下監(jiān)測(cè)點(diǎn)氚氣濃度繼續(xù)增加（4.7~35 s），35 s后氚濃度值穩(wěn)定，如圖7（a）所示，說(shuō)明停止泄漏后，氚氣在初始速度的作用下繼續(xù)流動(dòng)至各角落，直至均勻分布于場(chǎng)景空間；通風(fēng)環(huán)境下氚氣濃度先增加后減小，其中P4監(jiān)測(cè)點(diǎn)在17.5 s時(shí)，濃度值驟然減小，且明顯低于其他各監(jiān)測(cè)點(diǎn)，如圖7（b）所示，說(shuō)明停止泄漏后，通風(fēng)環(huán)境加劇了氚氣的流動(dòng)，使得停留在P4監(jiān)測(cè)點(diǎn)的氚氣較少。

圖7 空間監(jiān)測(cè)點(diǎn)位氚氣濃度隨時(shí)間的變化 (a)密閉環(huán)境，(b)通風(fēng)環(huán)境Fig.7 Variation of tritium concentration with time at space monitoring points(a)Herimetic environment,(b)Ventilated environment

圖8為通風(fēng)環(huán)境泄漏源上方、進(jìn)風(fēng)口、排風(fēng)口監(jiān)測(cè)點(diǎn)氚氣濃度隨時(shí)間變化規(guī)律，當(dāng)儲(chǔ)氚容器發(fā)生泄漏時(shí)，氚氣濃度在氣體泄漏方向（容器正上端）迅速增加，隨著容器內(nèi)壓強(qiáng)急劇減小，氣體泄漏速度逐漸減小，氚氣流動(dòng)偏離泄漏方向，沿容器正上端氚氣濃度迅速減小，如圖8（a）所示；場(chǎng)景空間中進(jìn)、排風(fēng)口受通風(fēng)條件影響，在容器發(fā)生泄漏階段，氚氣濃度出現(xiàn)波動(dòng)，泄漏停止后，排風(fēng)口氚氣濃度逐漸減小，而進(jìn)風(fēng)口氚氣濃度呈斷崖式降低，如圖8（b）所示，說(shuō)明當(dāng)場(chǎng)景空間內(nèi)發(fā)生氚泄漏時(shí)，通風(fēng)能有效達(dá)到除氚目的，且空間內(nèi)存在濃度相對(duì)較低區(qū)域。

圖8 通風(fēng)環(huán)境空間監(jiān)測(cè)點(diǎn)氚氣濃度隨時(shí)間的變化 (a)泄漏源上方，(b)進(jìn)、排風(fēng)口Fig.8 Variation of tritium concentration with time at space monitoring points under ventilation environment(a)Above the leakage source,(b)Inlet and exhaust vents

綜上所述，通風(fēng)條件下，當(dāng)氚氣停止泄漏后，進(jìn)風(fēng)口處氚氣濃度值小于1×10-6kg·m-3，位于進(jìn)風(fēng)口下端P4 監(jiān)測(cè)點(diǎn)處氚氣濃度小于1×10-2kg·m-3，兩個(gè)區(qū)域?qū)儆谙鄬?duì)安全區(qū)域，如圖9所示。

圖9 x=2 m處YZ平面的氚氣濃度分布Fig.9 Tritium concentration distribution in YZ plane at x=2 m

3 結(jié)語(yǔ)

本文針對(duì)高壓儲(chǔ)氚容器氚泄漏問(wèn)題，開展了基于AN-EOS 模型的數(shù)值模擬研究，描述了氚泄漏瞬態(tài)行為，構(gòu)建了典型場(chǎng)景模型，模擬了氚氣泄漏過(guò)程，分析了密閉和通風(fēng)環(huán)境下氚氣濃度的時(shí)空分布規(guī)律，得出如下結(jié)論：

1）高壓儲(chǔ)氚容器發(fā)生泄漏時(shí)，氚氣流動(dòng)經(jīng)歷了兩個(gè)階段：音速和亞音速流動(dòng)階段，音速流動(dòng)階段氚氣質(zhì)量和泄漏速度隨時(shí)間非線性變化，亞音速流動(dòng)階段氚氣泄漏速度隨時(shí)間近似線性變化。

2）發(fā)生氚泄漏時(shí)，密閉環(huán)境中氚氣流動(dòng)主要受泄漏時(shí)初始速度影響。通風(fēng)環(huán)境中氚氣流動(dòng)主要受空氣流動(dòng)影響；且通風(fēng)條件下空間內(nèi)監(jiān)測(cè)點(diǎn)氚氣濃度值明顯低于密閉條件，通風(fēng)能有效達(dá)到除氚目的；進(jìn)風(fēng)口區(qū)域氚氣濃度值明顯低于其他區(qū)域，屬于相對(duì)安全域。

氚泄漏事故下，影響其濃度分布的因素復(fù)雜，本文對(duì)儲(chǔ)氚容器變壓和不同工況下的氚泄漏問(wèn)題進(jìn)行了研究，得出了一些有意義的結(jié)論，但進(jìn)/排風(fēng)口的位置、排風(fēng)速度、儲(chǔ)氚容器體積和壓強(qiáng)等因素必然會(huì)對(duì)氚泄漏后的分布情況帶來(lái)影響，是有待進(jìn)一步深入研究的重點(diǎn)問(wèn)題。

作者貢獻(xiàn)聲明李昌均：模擬計(jì)算，起草文章；蔡幸福：論文整體設(shè)計(jì)，數(shù)據(jù)分析，研究經(jīng)費(fèi)支持；肖明清：選題指導(dǎo)；王飛：理論指導(dǎo)，文章修改；曹曉巖：仿真計(jì)算指導(dǎo)，文章修改。