高压氚气泄漏事故模拟与应急通风优化方法

蔡幸福， 李昌均,2， 薛冬林， 霍勇刚， 王浩伟

(1.火箭军工程大学 核工程学院，陕西 西安 710025； 2.空军工程大学 研究生院，陕西 西安 710051； 3.火箭军装备部装备项目管理中心，陕西 西安 710051)

氚是氢的放射性同位素[1]，其化学性质活泼，易发生氧化和置换反应，渗透能力强，被广泛应用于军事和民用领域，在氚气的生产、贮存和使用过程中，由于自然灾害或人为失误等原因，氚气包容物可能因意外力的作用产生破损，导致氚泄漏事故发生，一旦发生氚泄漏事故，将给人员和环境带来放射性危害[2]。因此，氚安全是人类有效利用氚的核心内容，开展事故情况下的氚气泄漏、扩散行为及应急处置方法研究，对核安全管理、核事故应急处置具有重要现实意义和理论价值，是当前氚安全研究的热点问题[3]。

国内外主要从试验和模拟2方面开展了氚泄漏行为和应急处置研究。在试验研究方面，主要以美国洛斯阿拉莫斯国家实验室的氚系统测试组件(tritium system test assembly，TSTA)和日本原子能院的氚安全Caisson组件(caisson assembly for tritium safety，CATS)为主，TSTA项目始于1978年，试验装置体积约为3 000 m3，旨在开发、设计和演示聚变反应堆氚处理系统的技术和安全运行状况，开展了多次氚气有意释放试验，得到了不同组分、不同氚量、不同泄漏位置、不同通风条件下的氚气泄漏、扩散及除氚效率数据，2003年退役[4-8]；CATS是美日合作项目，建造于1998年，试验装置体积为12.07 m3，旨在开展大量的聚变反应堆安全试验研究，如氚限制性能演示验证、动态模拟代码开发、早期氚行为数据积累以及新安全设备功能测试等[9-18]。在模拟研究方面，主要基于TSTA和CATS试验数据，开展了计算机数值模拟研究，先后开发和使用的软件程序包括：CFX-F3D[5]、FLOW-3D[6]、TBEHAVIOR、GASFLOW-II和Fluent[13]，验证的氚气扩散模型包括：RNG模型[14]、Launder-Sharma模型[15]、Launder-Spalding模型[16]、k-ε湍流模型[17]等，氚气早期释放行为均采用停留时间函数进行描述[18]；国内也针对中国聚变工程试验堆开展了氚泄漏数值模拟研究，LI Wei等[19]基于Fluent软件的开展了氚缓慢泄漏与扩散的数值模拟研究，获取了开放空间、密闭空间2种不同场景的亚音速氚射流、速度剖面、氚扩散距离以及沿射流轨迹的氚浓度；WEI Shiping等[20-21]研究了TEP氚处理管线双端断裂叠加手套箱破口的氚泄漏事故，评估了通风、温湿度、同位素交换率及壁面材料对除氚效率的影响；南华大学建立了氚提取系统手套箱和操作间的实际尺寸三维几何模型，使用Fluent软件模拟了逆流回热器和扩散器的氚泄漏和扩散，分析了通风对氚泄漏、扩散和浓度分布的影响，并评估了操作间的辐射危险[22-23]；哈尔滨工程大学设计开发了不确定性分析平台CUSA，也可为氚泄漏与扩散精确数值模拟提供有意义的借鉴[24]。

目前，绝大多数研究成果集中在聚变堆氚气的缓慢泄漏与扩散行为研究，对高压容器氚泄漏的研究较少，解放军火箭军工程大学在该方面作出了初步研究成果[25-27]，仍需从应急处置角度出发，给出特定场景下最佳通风位置及所需的最短通风时间。本文针对高压贮氚容器氚泄漏问题，建立了基于AN-EOS的氚气泄漏行为模型，编写用户定义函数(user defined function，UDF)接口代码，实现与Fluent求解器耦合计算，并开展氚泄漏事故数值模拟和应急通风优化研究，给出最佳通风位置下的最短应急通风时间。

1 高压氚气泄漏行为模型

1.1 AN-EOS模型

描述高压气体泄漏的热力学模型主要有2类：理想气体模型和实际气体状态方程模型[28]，理想气体模型假定气体压缩因子为1，实际气体状态方程是经过大量实验对理想气体状态方程修正后的半经验模型。AN-EOS模型是在理想气体模型基础上，考虑了气体分子体积的影响，对气体压缩因子修正后得到的气体状态方程。在高压储存条件下，气体分子行为偏离理想气体行为，AN-EOS模型比理想气体模型更能准确地描述气体泄漏的物理特性[29]，因此，本文选用AN-EOS模型开展高压氚气泄漏的瞬态行为研究，AN-EOS模型形式为：

p(v-b)=RgT

(1)

式中：p为压强；v为比体积；b为气体分子比积体修正项，b=7.691×10-3m3/kg；Rg为气体常数；T为温度。

假定气体的泄漏过程为等熵过程，由AN-EOS模型可得：

p(v-b)k=c

(2)

式中：k为比热容比；c为常量。

定容积贮存容器高压气体泄漏过程主要分为2个阶段：音速流动阶段和亚音速流动阶段。容器内压强与环境压强之比大于临界压强比时，泄漏口处气流速度达到当地声速，此时为音速流动阶段，泄漏口处压强p2=pi·vcr，泄漏口处的质量流率为：

(3)

式中：Qm为泄漏口质量流率；A为泄漏口面积；p2为泄漏口压强；pi为容器内压强；vcr为临界压强比；v2为泄漏口处比体积。

容器内压强与环境压强之比小于或等于临界压强比时，泄漏口气流为亚音速流动，泄漏口压强减小到环境压强，即p2=pamb，此后保持不变，此时，泄漏口处的质量流率为：

(4)

式中：pamb为环境压强；Q′m为亚音速阶段泄漏口质量流率；vi为容器内比体积。

1.2 模型验证

为了验证AN-EOS高压气体泄漏模型的准确性，将该模型计算结果与真实气体泄漏数据比对比分析，设定与文献[30]相同的实验条件：泄漏气体内部压强为34 MPa，气体温度300 K，泄漏孔面积为3.17×10-5m2，储气容器容积为2.73×10-2m3，环境压强为标准大气压。将模型参数按照实验条件设置求解，比对模拟结果与参考数据的变化情况，如图1所示。AN-EOS模型与真实气体在泄漏口压强和质量流率的变化趋势一致；表1给出了模拟结果与参考数据的误差分析，泄漏口压强最大相对误差为13.0%，最小为6.5%，质量流率最大相对误差为14.7%，最小为10.0%，随着气体泄漏时间增大，整体误差逐渐减小，说明本文提出的AN-EOS模型可用于开展高压氚气泄漏模拟研究。

图1 氚泄漏行为模型验证Fig.1 Tritium leakage behavior model validation

2 基于Fluent的氚泄漏事故数值模拟

本文采用Fluent软件开展氚泄漏数值模拟研究，主要包括几何模型构建、网格划分与边界条件设置、初始化及求解等步骤，针对高压贮氚容器的泄漏过程，采用用户定义函数UDF函数自定义边界条件，实现UDF与Fluent求解器的耦合计算。

2.1 几何模型构建

对真实尺寸长方体房间进行几何建模，如图2所示。房间长宽高分别为5、3和2.5 m，在正面墙体上方设置一个进风口，右侧墙体上方设置一个排风口，进/排风口尺寸均为0.5 m×0.5 m；高压贮氚容器简化为一个立方体，边长为0.3 m，内部压强为34 MPa，放置在场景地面中心位置，泄漏孔简化为一个圆形孔口，直径为10 mm，泄漏孔位于贮氚容器中心，坐标为(1.5，2.5，0.3)，氚气泄漏方向竖直朝向，泄漏过程中孔径不变；在空间设置8个监测点位T1～T8，各监测点坐标如表2所示。

图2 氚泄漏场景模型Fig.2 Tritium leakage scene model

2.2 网格划分与边界条件设置

在网格划分中，对整个物理模型采用六面体非结构化方法进行网格划分，对泄漏口和通风口区域进行加密处理，对壁面添加边界层，流体区域根据模型尺寸自动设置网格大小，网格最小尺寸为4.24×10-11m3，最大尺寸为2.8×10-4m3，网格数量约130万。划分的网格如图3所示，网格最小正交质量为0.78，满足质量要求。采用质量流率变量进行了网格无关性分析，分别计算泄漏后0.05 s和1 s的质量流率，前后误差分别为7.2%和5.3%，满足网格无关性要求。

图3 氚泄漏场景模型网格划分Fig.3 Tritium leakage scene model divided grid

在模型设置中，选择Realizablek-ε两方程湍流模型，模型中湍动能k的Prandtl数σk为1.0，散率ε的Prandtl数σε为1.2，经验常数C2为1.9，选择全浮力模型；在边界条件设置中，将泄漏口设置为质量入口，进风口设置为速度入口，排风口设置为压力出口，通过UDF函数定义入口质量流率、泄漏口压强和温度随时间的变化关系。设置组分气体为氚气，氚气摩尔质量为6.034 kg/kmol，比热容为7243 J/(kg·K)，热传导系数为0.138 1 W/(m·k)，粘性系数为1.259×10-5kg/(m·s)，扩散系数为7.41×10-6m2/s，设置房间壁面材料为碳酸钙，贮氚容器材料为钢铁，长方体空间内初始温度为300 K，温度为标准大气压。

表1 氚泄漏行为模型误差分析Table 1 Tritium leakage behavior model error analysis

表2 监测点坐标Table 2 Coordinates of monitoring points

2.3 氚气泄漏瞬态行为

为描述氚气泄漏瞬态行为，采用AN-EOS模型在Fluent求解器中进行计算，得到氚气质量、速度、质量流率及贮氚容器内压强随时间的变化曲线，如图4所示。初始时刻容器内氚气质量约为1.37 kg，氚气质量随时间非线性变化，当t=2.3 s时，氚气质量为0.137 kg，贮气容器泄漏量已达90%，当t=4.7 s时，氚气停止泄漏，容器内仍残存质量约为0.036 kg的氚气；初始时刻泄漏速度约为859.23 m/s，在音速流动阶段(0～3.6 s)，泄漏速度随时间非线性变化，在亚音速流动阶段(3.6～4.7 s)，则随时间近似线性变化，氚泄漏停止时，泄漏速度为0；质量流率及贮氚容器内压强随时间非线性变化，前2 s内，压强变化曲线斜率较大，压强由34 MPa降低至1 MPa；在2～4.7 s，斜率较小，变化较慢，由1 MPa降低至环境压强。图4表明：高压氚气泄漏主要集中于泄漏初期约2 s时间内，具有泄漏速度快、时间短的特点。

图4 氚气泄漏各物理参量随时间的变化曲线Fig.4 Time varying curves of physical parameters of tritium gas leakage

2.4 氚气浓度扩散规律

为描述通风环境下氚气浓度在场景空间的扩散规律，采用8个监测点分别测量不同时刻不同位置的氚气浓度，如图5所示。

图5 各监测点氚气浓度扩散规律Fig.5 Tritium gas concentration diffusion law at each monitoring point

图5(a)是场景四周1.5 m高处4个监测点T2～T5的氚气浓度随时间的变化规律，发生泄漏时(0～4.7 s)，各监测点位氚气浓度迅速增加；停止泄漏后(t>4.7 s)，各监测点氚浓度先增加后减小，由于T4监测点位于进风口下方，在17.5 s时，浓度值骤然减小，且明显低于其他各监测点。图5(a)是泄漏口正上方2个监测点T1和T8的氚气浓度随时间的变化规律，发生泄漏时，2个监测点的浓度迅速增加，泄漏停止后，两点浓度迅速减少，并最终与其他监测点的浓度达到基本平衡。图5(c)是通风进出口处氚气浓度随时间的变化规律，发生泄漏时，氚气浓度出现波动，停止泄漏后，排风口处T7的氚气浓度逐渐减小，而进风口处T6的氚气浓度则呈断崖式降低。图5的氚气浓度扩散规律说明：1)短暂的氚泄漏停止后，通风可以加剧了场景内氚气的流动，起到较好的除氚效果；2)通风口及其下方的区域氚浓度水平相对较低，易形成场景内的相对安全区域，x=2 m处YZ平面的氚气浓度分布云图如图6所示。

图6 氚气浓度分布情况Fig.6 Distribution of tritium concentration

3 应急通风优化方法

通风已被证实为氚泄漏事故下的有效应急处置措施，特定场景下的最佳通风位置及相应的最少通风时间是工程实践中关心的重要问题。

3.1 最佳通风位置

为得到实际尺寸场景空间的最佳通风位置，区分顶部单进单出、侧壁双进双出、前壁双进双出3种情况，进行了仿真计算，如图7所示，各通风口位置参数如表3所示，分别从流体速度流线、氚浓度分布和氚浓度降至本底水平的时间评估除氚效率，得到最佳通风位置。为保证相同的通风效率，单进单出和多进多出通风口的进/排风口总面积相同，单进单出进/排风口尺寸为500 mm×500 mm，双进双出进/排风口尺寸为354 mm×354 mm，各种情况的排风量均为1 m3/s。

图7 3种不同通风位置Fig.7 Three different ventilation positions

表3 通风口位置坐标Table 3 Coordinates of air inlet and outlet

图8给出了3种不同通风位置的除氚效果对比图。图8(a)为经进风口至排风口的50条速度流线，从分布来看，有限空间流体域内形成了数量、规模不等的涡流，双进双出比单进单出形成的涡流较多，说明内部气体流动更加剧烈，其中侧壁双进双出时涡流数量更多，尺寸更小，气体流动更剧烈，更易有效将氚气排出；图8(b)为氚气泄漏400 s时各剖面的浓度分布，双进双出比单进单出场景下的浓度水平低5个数量级，除氚效率更好，侧壁与前壁2种情况的差别不大；图8(c)为氚气浓度随时间的变化关系，其中氚气浓度进行了对数处理，依据文献[18]选取本底水平为3×105Bq/m3(1.2×10-12kg/m3)，可见顶部单进单出通风除氚时间为860 s，侧壁双进双出为830 s，前壁双进双出为850 s。由此可见，当高压贮氚容器位于地面中间位置，泄漏孔位于容器上表面中心位置，泄漏方向竖直朝上时：1)在通风量不变的情况下，增加通风口数量，可以促进空间内氚气流动，提高除氚效率；2)长方体场景空间中，最佳通风位置为侧壁双进双出，即2个进风口置于前壁，两个排风口置于侧壁。

图8 3种通风位置除氚效果对比图Fig.8 Comparison of tritium removal effects at three ventilation positions

3.2 最少通风时间

依据图7(b)侧壁双进双出的最佳通风位置，考虑实际通风系统最小功率、平均功率和最大功率3种运行工作模式，设置4、6、8 m/s共3种通风速率，计算空间内氚气浓度随时间的变化关系，得到不同速率下氚浓度降至本底所需的最少通风时间。

图9给出了3种不同通风速率下的除氚效果对比图。图9(a)和(b)分别为泄漏后100和400 s时空间内的氚气浓度分布，通风100 s时，8 m/s通风速率下空间内的氚气浓度比6 m/s低3倍，比4 m/s约低1个数量级；通风400 s时，8 m/s通风速率下的氚气浓度比6 m/s约低70倍，比4 m/s低5个数量级；图9(c)为3种通风速率下，氚气浓度降至本底水平时的分布情况，其中8 m/s通风速率下空间内氚浓度降至本底所需的最少通风时间为430 s，6 m/s通风速率下为550 s，4 m/s通风速率下为850 s。由此可见：1)通风速率越大，空间内涡流分布越多，氚浓度下降速度越快，通风速率每增大1倍，氚浓度将呈数量级下降；2)给出了不同通风速率下氚气浓度降至本底水平所需的最少通风时间，为真实场景氚泄漏事故应急通风处置提供了技术支持。

图9 不同通风速率下的除氚效果对比图Fig.9 Comparison of tritium removal effects under different ventilation rates

4 结论

1)本文提出了基于Fluent的氚泄漏事故数值模拟方法，通过UDF函数定义边界条件，实现了UDF与Fluent求解器的耦合计算，给出了氚气瞬态泄漏行为和特定场景空间的扩散规律，高压氚气泄漏过程主要集中在初期的2 s内，通风可以加剧氚气的流动，起到较好的除氚效果，并在通风口及其下方区域形成相对安全区域。

2)本文提出了应急通风优化方法，优化得到长方体空间的最佳通风位置为侧壁双进双出，给出了不同通风速率下氚气浓度降至本底水平所需的最少通风时间，为真实场景氚泄漏事故应急通风处置提供了技术支持。

本文在应急通风优化方法方面的研究主要侧重于数值模拟，相关基础理论的研究是需要进一步关注的重点问题。