基于VOF模型的偏二甲肼泄漏液池扩展过程的数值模拟

詹 翔，崔村燕，周宵灯，赵继广，辛腾达

(1.航天工程大学研究生管理大队，北京 101416；2.航天工程大学宇航科学与技术系，

北京 101416；3.航天工程大学电子与光学工程系，北京 101416)

偏二甲肼(Unsymmetrical Dimethylhydrazine,UDMH)易燃、易爆、剧毒、腐蚀性强，在航天发射场使用的安全问题历来受到人们的高度关注[1]。肼类燃料对非金属的溶胀作用很强烈，因此法兰连接处的密封、阀门和泵填函部位的密封长时间使用后容易失效[2]。实际上在偏二甲肼储存和加注期间，容器阀门法兰腐蚀、设备失效、焊接缺陷以及人员的误操作等原因[3]都会造成推进剂不同程度的泄漏。偏二甲肼泄漏后在地面上形成液池，遇到火星易着火爆炸，也容易蒸发形成有毒气体接触到人员和设备,可能对其造成损害[4]，严重威胁着周围设备和工作人员的安全。而研究偏二甲肼泄漏液池的扩展规律是进一步研究推进剂蒸气扩散规律、开展推进剂安全评价工作、制定推进剂泄漏事故应急预案等的重要前提。在以往对偏二甲肼泄漏液池扩展过程的研究中[5]，对偏二甲肼泄漏液池做了极简处理，即简单地以圆形液池代替，这样并不能真实地反映偏二甲肼泄漏液池的分布规律。为此，本文采用VOF模型对偏二甲肼泄漏形成液池的扩展过程进行了数值模拟，研究了在受限空间内外泄偏二甲肼液池的扩展规律。

1 数值模拟方法

1. 1 VOF模型

数值模拟软件FLUENT中的多相流模型可以模拟不同流体间的相互流动，其中VOF模型通过自由液面的重构来获得自由液面的位置和形状[6]，是处理液气两相交界、互相推动力瞬变问题的重要方法[7]。VOF模型采用k-ε紊流模型，其基本的控制方程包括连续方程、动量方程、湍流方程和湍动能耗散率方程[8]。利用VOF模型追踪液体的自由表面是数值模拟方法中较为有效的方法，目前已有大量的研究，如陈群等[9]通过对阶梯溢流坝建模，对溢流坝面流场的紊流进行了数值模拟，得到的数值模拟结果与实测值相比分布规律近似，数据相差较小；关大玮等[10]采用FLOW3D软件，将RNGk-ε模型和VOF模型结合，对溢流坝三维水流进行了数值模拟，模拟结果对照缩比物理模型试验结果大体上吻合；史晓蒙等[11]通过建立库区地面输油管道泄漏扩散过程的三维CFD模型，采用FLUENT软件中的VOF模型对油品泄漏扩散过程进行了数值模拟，分析了油品泄漏后的扩散趋势，该研究结果有助于分析泄漏事故后果影响范围以及合理制定防护及救援措施。

1. 2 自由液面的重构

利用数值模拟方法研究偏二甲肼泄漏液池的扩展规律，重点在于处理偏二甲肼和空气的相交界面。FLUENT软件中的控制体积公式需要计算穿过控制体积面的对流和扩散通量，并且与控制体积内部的源项平衡[12]。划分好网格后FLUENT软件即可按给定的初始运动界面开始计算，但是能捕捉到的新运动界面是一个模糊的轮廓。为了提高界面的分辨率和精度，FLUENT软件必须在每个时间步长计算完成后对界面进行重构。对于VOF模型，FLUENT软件中有4种方案计算流体界面的通量：几何重构方案、施主受主方案、欧拉显式方案和隐式插补方案[13]。其中，几何重构方案求解简单，本文主要采用几何重构方案。

在几何重建方案中，若单元被某一项填满，那么FLUEN软件将以标准插值方案计算流体界面的通量，如果单元包含两相之间的交界面，则使用几何重建方案。几何重建方案利用分段线性的方法描绘流体之间的界面，见图1。

图1 真实流体界面与几何重构流体界面Fig.1 Real fluid interface & geometric reconstruction

由图1可以明显地看出，真实流体的界面是圆滑的弧形界面[见图1(a)]；而几何重构方案则将原来的圆滑弧形界面在单元格内由线段代替[见图1(b)]。

假定两种流体之间的界面在每个单元内有线性的斜面，并使用这个线性界面形状来计算穿过单元面之间的对流，在计算流体界面通量时，首先根据体积分数和单元的相关信息来计算出线性界面的位置，再通过界面上的切向速度分布信息来计算每个流体界面的通量。

当实际交界面的尺度相对于网格尺度很小时，几何重构出的线性界面将无法反映真实的流体交界面，见图2。

图2 交界面与网格尺寸的对比图Fig.2 Comparison of interface and grid size

由图2可见，气泡的大小比单个计算单元格还小，若按照几何重构方案，在一个单元格内用一条线段代替交界面，则无法正常刻画气泡的轮廓。Rabha等[14]认为网格数与气泡直径的比值至少为16时才能较好地刻画出气泡的轮廓。本文为了较好地刻画偏二甲肼泄漏过程中出现的液滴，在局部区域做了网络加密处理。

2 几何模型及边界条件设置

以某偏二甲肼容器参数为依据，对容器实际外形条件适当简化后建立几何模型。某偏二甲肼容器高35 m、半径为1.675 m，假定加泄连接阀处法兰腐蚀发生连续泄漏，并通过容器外壁面上2 m处的连接口向外泄漏。由于泄漏口以上的部分在实际计算中对偏二甲肼的泄漏几乎没有影响,因此在不影响数值模拟结果的前提下对几何空间进行了合理简化，只留下偏二甲肼容器的部分弧形壁面和泄漏口正前端的部分空间，见图3。

图3 总体结构与计算区域Fig.3 Structure and calculation area

根据实际要求，受限空间环境参数需控制在温度为15℃，相对湿度为50%，风速为0.5 m/s。以偏二甲肼容器加泄连接阀处法兰接触面发生泄漏为例，偏二甲肼泄漏量的计算公式为[15]

(1)

式中:Q为偏二甲肼泄漏量(cm3/s)；h为间隙的大小(cm)；Δp为设备内外的压力差(MPa)；η为偏二甲肼的黏度(Pa·s)；r1为法兰内径(cm)；r2为法兰外径(cm)。

偏二甲肼容器内部为正常大气压，求得其泄漏速率为45.8 cm3/s。假定忽略泄漏过程中容器内压力的变化，则泄漏速率保持不变，进而可根据加泄连接口面积大小获得出口速度为0.03 m/s。本次模拟假设偏二甲肼容器壁面为绝热、无滑移，且流体区域不可压缩，呈湍流状态，湍流分布符合伯努利假设；由于空间有限，假设空间内空气分布均匀，不考虑空气的分层和稳定度。

3 偏二甲肼泄漏液池扩展过程的数值模拟与分析

3. 1 偏二甲肼泄漏液池的分布

不同时刻偏二甲肼泄漏液池分布的数值模拟结果见图4。

图4 不同时刻偏二甲肼泄漏液池的分布Fig.4 UDMH liquid pool distribution at different times

由图4可见，当t=0.15 s时，偏二甲肼在泄漏口附近聚集,在重力作用下，一部分偏二甲肼从泄漏液主体分离，单独沿壁面下落;当t=1.50 s时，下落的偏二甲肼泄漏液分裂成面积约为3 cm2的液滴，散布在壁面上；当t=3.00 s时，偏二甲肼泄漏液的一部分已经到达地面；当t=4.05 s时，偏二甲肼在地面上聚集形成面积大小约为224.35 cm2的液池。

3. 2 垂直壁面上偏二甲肼泄漏的流态分析

图5给出了不同时刻垂直壁面上偏二甲肼泄漏液池面积的增长曲线图。

图5 不同时刻垂直壁面上偏二甲肼泄漏液池面积的增长曲线图Fig.5 Growth curve of the UDMH liquid pool area on the vertical wall at different times

由图5可见，在垂直壁面上偏二甲肼泄漏液池的扩展过程大致可分为A、B、C三个过程：A过程中，偏二甲肼泄漏的时间较短，泄漏量较少，泄漏液体积小，液体黏性力和表面张力对液池扩展起到了一定的阻碍作用[16]，此时泄漏液受重力的加速作用还不明显，在泄漏口附近聚集并缓慢下落，液池面积增长速度比较慢；B过程中，一部分液体在重力的作用下获得了相对高的下落速度，在壁面摩擦力的作用下，分散成多个小液滴，小液滴在壁面上横向分散，因此该过程中液池扩展速度迅速增加，且在B过程末期，偏二甲肼泄漏液不再形成新的液滴，壁面上偏二甲肼液池呈现下宽上窄的分布规律，此时液池的面积增长速率也逐渐放缓；偏二甲肼泄漏液到达地面时，开始进入C过程，该过程中由于偏二甲肼液滴到达地面且无新的液滴产生，垂直壁面上偏二甲肼泄漏液池的总面积不断减少，且C过程结束后所有分散的液滴达到地面，当垂直壁面上只剩下垂直向下且细长的偏二甲肼泄漏液池时，偏二甲肼泄漏液池的面积将基本不再变化。

图6给出了偏二甲肼泄漏液池具有A、B、C三个过程典型特征的流态图。

图6 不同时刻垂直壁面上偏二甲肼泄漏的流态图Fig.6 Flow pattern of the UDMH leakage on the vertical wall at different times

当垂直壁面上偏二甲肼泄漏液池的面积不再变化后，由数值模拟结果可得到偏二甲肼泄漏液池的面积约为301 cm2。偏二甲肼蒸发速率的计算公式为[17]：

vu=2.603 6Aγ0.553×(0.001 935T2+0.018 988T+0.511 81)×(-2H2+4.46H+0.026 5)

(2)

式中：vu为偏二甲肼的蒸发速率(mg/min)；A为偏二甲肼泄漏液池的面积(m2)；γ为风速(m/s)；T为温度(℃)；H为相对湿度。

由以上环境条件，可计算出在垂直壁面上偏二甲肼的蒸发速率约为0.119 mg/min。

3.3 偏二甲肼泄漏在地面上形成的液池的扩展过程分析

偏二甲肼泄漏到达地面后紧靠容器壁面形成半圆形的液池，图7给出了不同时刻地面上偏二甲肼泄漏液池面积的增长曲线。

图7 不同时刻地面上偏二甲肼泄漏液池面积的增长曲线图Fig.7 Growth curve of the liquid pool area of UDMH leakage on the ground

由图7可见，最初少量接触到地面的偏二甲肼泄漏液汇聚在一起，由于液体自身黏性力和地面摩擦力的阻碍作用,此时偏二甲肼泄漏液池面积的增长速度较慢；随着偏二甲肼在地面上聚集量的增加，偏二甲肼泄漏液池在自身重力和下落滴液冲击力的作用下克服液体自身黏性力和地面摩擦力，形成面积更小的圆形液池向四周扩散，此时偏二甲肼泄漏液池面积的增长速度也有所增加。图8给出了不同时刻地面上偏二甲肼泄漏液池的扩展过程。

图8 不同时刻地面上偏二甲肼泄漏液池的扩展过程Fig.8 Expansion process of the liquid pool of UDMH leakage on the ground at different times

3. 4 泄漏速率对偏二甲肼漏泄漏液池扩展的影响

为了研究不同泄漏速率对偏二甲肼泄漏液池扩展规律的影响，将泄漏速度加倍即出口速度改为0.06 m/s，其他条件不变进行了数值模拟计算，得到同时刻(3.00 s)不同泄漏速度条件下垂直壁面上偏二甲肼泄漏液池的扩展情况对比，见图9。

图9 同时刻不同泄漏速度下垂直壁面上偏二甲肼 泄漏液池的扩散情况对比Fig.9 Comparison of the expansion of UDMH liquid pool on container wall at the same time with different leakage velocity

由图9可见，不同泄漏速度对垂直壁面上偏二甲肼泄漏液池的形状影响较小，在泄漏速度较大的情况下，垂直壁面上分散的小液滴数量明显减少。由此可见，不同的泄漏速度条件下，垂直壁面上仍然是形成垂直向下且细长的偏二甲肼泄漏液池。

图10为不同泄漏速度条件下地面上偏二甲肼泄漏液池的扩展面积对比。

由图10可见，泄漏速度对地面上偏二甲肼泄漏液池扩展面积的影响十分明显： 当泄漏速度为原来的两倍时，偏二甲肼泄漏液池的扩展面积是原速度情况下的4倍左右；当泄漏量较大时，偏二甲肼泄漏液池更加接近半圆形，并且液池向四周扩散时分裂形成小面积液池的过程进行得更快。这主要是因为泄漏量的增加，导致泄漏液能获得更大的动能，以克服地面摩擦力和液体自身黏性力向四周扩散。

图10 不同泄漏速度下地面上偏二甲肼泄漏液池扩展面积的对比Fig.10 Comparison of the expansion area of UDMH liquid pool on the ground with different leakage velocity

4 结 论

本文通过对偏二甲肼泄漏液池的扩展规律进行数值模拟与分析，得到如下结论：

(1) 偏二甲肼泄漏时，在垂直壁面上主要形成垂直向下且细长的液池；液池前端部分在壁面摩擦力的作用下分散成小的液滴，由于液滴向左右扩散，垂直壁面上偏二甲肼泄漏液池呈现下宽上窄的分布规律；当垂直壁面上所有液滴均到达地面后，壁面上的液池垂直向下流动且宽度基本一致，其面积约为301 cm2。根据环境条件，可求出该液池的偏二甲肼蒸发速率为0.119 mg/min。由此可知，在对偏二甲肼泄漏进行危害性分析时，不能忽视垂直壁面上液池的蒸发作用。

(2) 偏二甲肼泄漏时，地面上的泄漏液在偏二甲肼容器壁面周围聚集，形成近似半圆形的液池；随着泄漏液的不断积累，半圆形液池向周围扩展的过程中会分裂出面积较小的圆形液池。

(3) 偏二甲肼泄漏时，当泄漏速度加大后，垂直壁面上零散的液滴明显减少，地面上形成的液池更加接近规则的半圆形;当泄漏速度为原来的两倍时，地面上偏二甲肼泄漏液池的扩展面积是原速度情况下的4倍左右。