2. 清华大学 热科学与动力工程教育部重点实验室, 北京 100084
2. Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Tsinghua University, Beijing 100084, China
作为一种极具潜力的可再生能源,氢气具有清洁、高效、单位质量能量密度高和来源广泛等优点。以氢燃料电池汽车为代表的氢能研究和应用越来越受到关注[1]。然而,在氢能的大规模商业化应用之前,必须全面系统地研究氢安全问题,以推动相关安全标准和技术规范的制定[2]。典型的氢安全事故序列往往包括意外泄漏、扩散、点燃、爆燃或爆炸等过程[3]。其中,氢气泄漏和扩散阶段的特性对于后续事故现象的发展和事故应急处理都具有决定性的作用。氢气泄漏一般通过等效直径非常小的裂口(毫米级),以高速射流的形式发生,因此属于小孔泄漏的范畴,目前已有许多针对高压氢气射流的理论和实验研究[4-6],但是理论模型仅适用于自由空间射流分析,难以用于限制空间或有障碍物存在时的情况;实验研究测试工况有限,很难同时测量射流的流场和浓度场,且氢气射流实验所需的安全保障成本较高。因此,在一定的实验验证的基础上进行数值模拟研究就成为完善氢安全研究的重要途径和必然选择。
高压氢气泄漏射流的数值模拟存在两个主要的难点:一是泄漏点外剧烈变化气流参数使得数值计算难以稳定和收敛;二是泄漏点的尺度(毫米级)与整个计算域的尺度(米级)之间的差距较大,且泄漏点附近需要特别精细的网格以捕捉气流参数剧烈变化的激波区域,从而导致整个计算域网格单元数目巨大。虽然目前已有一些完整的高压氢气泄漏模拟研究[7-9],但是由于上述困难的存在,每个算例都需要耗费大量的计算资源,制约了数值模拟方法在高压氢气泄漏研究中的应用。为解决上述问题,文献中存在两种方案:一是Xu等[9]提出的“两步法”,即先计算激波区,再将计算结果作为入口条件计算下游的扩散区,这种方法对计算效率的提升不明显;二是Ruggles等[10]、Xiao等[11]利用虚喷管模型来避免对激波区的计算,从而显著提高了计算效率,但此方法的模型假设不合理,计算结果也不够准确。
为解决上述问题,本文在实验测量高压氢气射流激波结构的基础上,结合守恒方程,提出了分层流动模型以简化高压氢气射流数值模拟,并将简化模拟的计算结果与实验测量值和完整模拟的计算结果进行了对比。
1 实验研究本文的实验研究包括两部分:一是利用纹影法对高压氢气射流激波结构进行测量,以提供模型中所需的激波特征尺寸;二是利用平面激光Rayleigh散射(PLRS)技术对高压氢气射流浓度场进行可视化测量,以检验模型的计算结果。
1.1 激波结构测量实验使用特制的LED灯作为光源,可以产生波长为520 nm的绿色单色光,从而消除光源色差对拍摄结果的影响。在LED灯前方放置凸透镜和平凸透镜,以将光源产生的光线聚焦。在平凸透镜的焦点处安装光阑,从而形成实验所需的点光源。通过光阑的光线随后透过准直镜形成平行光,并通过射流气体,纹影镜将平行入射的光线聚焦到水平放置的刀口上,通过调整刀口的高度可以改变最终进入相机镜头的光强度。最后,利用CCD相机拍摄纹影图像并输入计算机。采用出口直径de为1 mm的喷嘴,分别进行了储存压力p0为1~5 MPa的氢气射流实验。
实验测得的储存压力为5 MPa的氢气射流激波结构如图 1所示。可以清晰地看到呈圆桶状的射流核心区、包围核心区的边界层区和明亮的Mach盘,以及Mach盘下游的反射波。其他储存压力的氢气射流激波结构与之类似。
通过对纹影法实验测量结果的分析,得到Mach盘的位置zm、直径dm和边界层厚度BS的经验公式如下:
$ {z_{\rm{m}}}/{d_{\rm{e}}} = 0.67\sqrt {{p_0}/{p_\infty }} , $ | (1) |
$ {d_{\rm{m}}}/{d_{\rm{e}}} = 0.35\sqrt {{p_0}/{p_\infty }} , $ | (2) |
$ {B_{\rm{S}}}/{d_{\rm{e}}} = 0.30\sqrt {{p_0}/{p_\infty }} . $ | (3) |
其中:p为压力,下标e表示喷嘴出口,下标∞表示周围环境。图像数据处理方法可参考文[4, 12]。
1.2 射流浓度场测量Rayleigh散射是指当光或其他电磁波通过分子直径比其波长小得多(小于入射光波长的1/10)的透明介质时发生的弹性散射。当入射光的强度一定时,散射光强度正比于介质的分子数密度,而在气体流动中,介质的分子数密度就反映了流体的密度分布情况。当气体为混合气体时,混合气体的总体密度取决于各气体组分的体积分数,因此可以利用不同的Rayleigh散射光强度来计算混合气体的密度,进而计算出各组分气体的浓度。
射流浓度场利用PLRS技术进行非接触式测量,从而避免测量设备对气流的干扰。实验采用波长为532 nm的Nd:YAG激光,激光束由激光器生成并经过反射镜调整方向,随后相继通过一系列透镜系统形成平面激光束。激光束在通过实验段后进入特制的收集器,以消除反射对测量结果的影响。采用CCD相机来收集Rayleigh散射的光信号,同时由计算机来记录摄像机拍摄的图像。随后,采用自行编写的图像处理程序对获取的光强度图像进行处理,以获得射流的气体的摩尔分数,具体的图像处理方法可参考文[13]。实验所测得的气流浓度值将用于对后文数值模拟计算结果的检验。
2 数值计算模型 2.1 分层流动模型高压氢气射流激波结构如图 2所示。储罐内(位置0)的气体经过喷嘴出口(位置1)泄漏到环境空气中,形成了核心区和边界层区两个流区。在核心区内的气流经过加速膨胀后达到最大流速(位置2a),然后经过Mach盘后减速为亚声速流(位置2b),而边界层区内的氢气与空气混合气流在经过Mach盘后仍保持很高的流速(位置3)。模型假设气体由滞止状态到喷嘴出口之间为绝热膨胀过程,器壁与气流之间的热量传递与气流的焓相比非常小,可以忽略不计。同理,气流由喷嘴出口绝热膨胀到Mach盘处,气流经过Mach盘后压力和温度恢复到环境压力和温度[14-15]。边界层区内的氢气和空气混合气流在到达位置3处压力等于环境压力。此外,为了简化计算,假设位置2和3处的气流速度均匀。此外,从氢安全研究的角度考虑,假定泄漏为稳态泄漏,即考虑泄漏浓度扩散范围最大、最不利于安全的情形。
为了考虑高压气体物性对理想气体的偏离,模型中使用如下Noble-Abel实际气体状态方程[16]:
$ p = \frac{{\rho {R_{\rm{g}}}T}}{{1 - b\rho }}. $ | (4) |
其中:p、ρ和T分别为气体的压力、密度和温度;Rg为气体常数;b为与气体种类有关的常数,对于氢气b=7.691×10-3 m3/kg。绝热膨胀关系式为
$ {p_0}{\left( {\frac{1}{{{\rho _0}}} - b} \right)^\gamma } = {p_1}{\left( {\frac{1}{{{\rho _1}}} - b} \right)^\gamma }. $ | (5) |
其中:γ为绝热指数,数字下标表示位置。结合式(4)和(5)即可得到喷嘴出口处的气流参数。
在喷嘴出口处的气流速度u1即为当地声速,可以根据声速公式和气流温度T1计算得到。射流核心区内的气流在紧靠Mach盘的上游处达到最大速度,此处的Ma可以通过求解如下方程得出
$ \frac{{{p_0}}}{{{p_\infty }}} = \frac{{{{\left( {1 + \frac{{\gamma - 1}}{2}Ma_{{\rm{2a}}}^2} \right)}^{\frac{\gamma }{{\gamma - 1}}}}}}{{\frac{{2\gamma }}{{\gamma + 1}}Ma_{{\rm{2a}}}^2 - \frac{{\gamma - 1}}{{\gamma + 1}}}}. $ | (6) |
利用正激波关系式可得紧靠Mach盘下游处的Ma
$ M{a_{{\rm{2b}}}} = {\left( {\frac{{2 + \left( {\gamma - 1} \right)Ma_{{\rm{2a}}}^2}}{{2\gamma Ma_{{\rm{2a}}}^2 - \left( {\gamma - 1} \right)}}} \right)^{1/2}}. $ | (7) |
进而可以根据当地声速计算公式得到u2b。Mach盘处边界层流区的截面积为
$ {A_3} = \pi {\left( {\frac{{{d_{\rm{m}}}}}{2} + {B_{\rm{S}}}} \right)^2} - {A_2}. $ | (8) |
其中A2为Mach盘的面积。
从喷嘴出口到Mach盘下游(位置2b和3)的质量、动量和能量守恒方程为:
$ {\dot m_1} = {\dot m_2} + (1 - {\omega _{{\rm{air}}}}){\rho _3}{A_3}{u_3}, $ | (9) |
$ {p_1}{A_1} + {\dot m_1}{u_1} = {p_{{\rm{2b}}}}{A_1} + {\dot m_2}{u_{{\rm{2b}}}} + {\rho _3}{A_3}u_3^2, $ | (10) |
$ \begin{array}{l} {{\dot m}_1}\left( {{c_{{\rm{p, gas}}}}{T_1} + \frac{{u_1^2}}{2}} \right) = {{\dot m}_2}\left( {{c_{{\rm{p, gas}}}}{T_{{\rm{2b}}}} + \frac{{u_{{\rm{2b}}}^2}}{2}} \right) + \\ \;\;\;\;\;\;{\rho _3}{A_3}{u_3}\left( {{c_{{\rm{p, 3}}}}{T_3} + \frac{{u_3^2}}{2} - {\omega _{{\rm{air}}}}{c_{{\rm{p, air}}}}{T_\infty }} \right). \end{array} $ | (11) |
其中:
利用分层流动模型可以计算得到Mach盘下游(位置2b)和Mach盘处边界层截面(位置3)等处的气流条件,并将这些条件结合相应的几何尺寸作为数值模拟的入口,可以避免计算喷嘴出口到Mach盘之间的复杂激波区,从而可以大大提高模拟效率。
由于混合气流和氢气流入口压力都与环境压力接近,整个计算区域内的压力变化也不大,所以模拟中采用基于压力的求解器和隐式算法。压力-速度耦合使用SIMPLE算法,动量和能量方程中的对流项使用二阶迎风差分格式,混合气体物性根据理想气体混合定律计算。湍流模型使用标准k-ε模型。Mach盘和边界层入口的边界条件均采用质量流量入口。计算采用二维轴对称几何模型,在整个计算域内使用结构化网格(四边形),网格总单元数为8万。此外,还对每个算例建立了总单元数为6万和10万的网格进行了网格独立性验证。由于计算使用基于压力的求解器,且网格总单元数较少,因此可以在普通的个人计算机上进行并行计算,每个算例仅需要花费几个计算机时即可完成。
2.2 完整模拟完整模拟的计算域包括图 2所示的所有部分,即入口条件为高压储氢罐内的压力入口。模拟采用基于密度的求解器和隐式算法,动量和能量方程中的对流项使用二阶迎风差分格式,混合气体物性根据气体混合定律计算。湍流模型使用标准k-ε模型。计算采用二维轴对称几何模型。由于喷嘴出口附近存在复杂的激波结构,因此这部分计算域需要非常精密的网格,而在距离出口较远区域使用比较粗糙的网格,以减少总的网格数量,提高计算效率。在计算开始时,采用非常小的Courant数(0.01~0.1)和松弛因子(0.1),以防止计算发散;当计算稳定后,特别是激波区已经形成并不再发生明显变化时,逐渐增大Courant数到2,增大松弛因子到0.8,以加快收敛速度。即使采用高性能计算机进行多核并行计算,完成每个算例也要花费几百个计算机时。采用总单元数为14万、17万和19万的3套网格进行了网格独立性检验。
2.3 虚喷管模型虚喷管模型的核心思想是假设一个流动参数均匀的等效出口,在这个等效出口处的射流质量流量与真实的喷嘴出口处的流量相等,根据守恒方程和模型假设计算等效出口的直径和流动参数[17-19]。目前常用的虚喷管模型往往假设在等效直径处气流速度为当地声速。然而,在实际的射流过程中并不存在这样一个位置,因此虚喷管模型的基本思路仅是通过假设流量相等来近似地评估射流远场的流场和浓度场,对于实际喷嘴附近存在障碍物时的情况,则难以应用。
作为对比,本文采用虚喷管模型中比较典型和常用的Ewan和Moodie模型[19]进行模拟计算。该模型假设等效出口处流速为当地声速,等效出口处气流温度与实际出口处相同。模拟设置与分层流动模型基本相同,采用总单元数为5万、7万和10万的网格进行了网格独立性验证。
3 计算结果 3.1 射流速度场完整模拟计算得到的以Ma云图表示的激波结构与纹影图像的对比如图 3所示。由于喷嘴的半径已知,通过调节图像的尺寸使模拟图像和纹影图像的像素/空间比一致,从而比较两者的激波结构尺寸。从图 3可见,完整数值模拟准确地预测了射流的激波结构,并准确地计算了Mach盘的位置和直径。由此可见,完整数值模拟的计算结果是较为可信的,实验中无法测量射流的速度时,可以使用完整模拟的计算结果来验证两种简化模型的计算结果。此外,由图 3中的Ma云图可见,Mach盘后射流核心区和边界层区存在明显的速度分层现象,即射流核心区内的流动在经过Mach盘后减速为亚声速,而此时边界层区内仍保持超声速流动。与轴向流动相比,两个区域间的速度掺混可以忽略,因此也说明了分层流动模型假设的合理性。对于射流速度分层情况的详细描述,可参考文[14]。
采用完整数值模拟和两种简化模拟计算得到的射流中心线的气流速度ucl如图 4所示。因为完整模拟中包含了喷嘴出口附近激波区,所以在距离出口较近的范围内计算得到的气流速度很高,因此在图 4中被纵轴的上限所截断。采用分层流动模型简化模拟计算得到的射流中心线速度先升高后降低,这是由于流速较低的核心区气流与流速较高的边界层区气流先相互混合使得流速增加,随后混合气流逐渐在空气中减速所造成的。采用虚喷管模型的计算结果明显偏离了完整模拟的计算结果,而采用分层流动模型的计算结果则与完整模拟基本一致。因此,在Mach盘之后,即可采用分层流动模型所计算的气流速度。
3.2 射流浓度场
采用完整数值模拟和两种简化模拟计算得到的沿射流中心线的氢气摩尔分数Xcl与实验测量值的对比如图 5所示。可以看到,完整模拟的计算结果略低于实验测量值,但整体上比较准确;而简化模拟方法倾向于高估沿射流中心线的氢气浓度。这可能是由于在简化模型中不存在实际的喷嘴,故在模型中没有考虑管道阻力造成的流量损失,因此模型计算的氢气流量略高于实际泄漏的氢气流量。
完整数值模拟和两种简化模拟计算得到的氢气摩尔分数为4%的最低可燃轮廓线与实验测量值的对比如图 6所示。由于射流激波区内涡结构不断生成和脱离,使得完整模拟的最低可燃轮廓线有较大的波动;而在分层流动模型中,由于模型所计算的出口平面上流速仅有轴向速度,忽略了射流在径向的掺混和流动,因此没有得到类似的涡结构。但是图 6仍能够说明完整模拟和分层流动模型简化模拟的计算结果都与实验测量值一致,特别是分层流动模型简化模拟的计算结果与测量值非常接近;而虚喷管模型简化模拟的计算结果则明显偏离了实验测量值和完整模拟的计算值。由于本文所采用的浓度测量数据都在距离喷嘴出口大于100de处,远大于Mach盘所在的位置,在此范围内采用分层流动模型所计算的浓度场是比较准确的。因此,综合考虑射流的速度场和浓度场,在距离喷嘴出口大于100de时,应用分层流动模型所计算的数据是较为合理的。
4 结论
本文分别利用纹影法和激光Rayleigh散射技术测量了高压氢气射流的激波结构和浓度场。在激波结构测量的基础上,结合实际气体状态方程和守恒方程建立了分层流动模型,以简化高压氢气射流的数值模拟研究。为了验证模型的有效性,将分层流动模型的计算结果与完整模拟、虚喷管模型和实验测量结果进行了比较。结果表明:完整模拟可以较为准确地计算射流的激波结构和浓度场;采用分层流动模型计算得到的射流中心线速度与完整模拟一致,浓度场计算结果与实验测量值也基本一致;而采用虚喷管模型计算得到的速度和浓度值都较为明显地偏离了实验测量值和完整模拟计算值。因此,综合考虑模型计算结果的准确性和计算效率,采用分层流动模型的简化模拟方法可以在保证结果准确性的基础上显著提高计算效率。
[1] |
MAZLOOMI K, GOMES C. Hydrogen as an energy carrier:Prospects and challenges[J]. Renewable and Sustainable Energy Reviews, 2012, 16(5): 3024-3033. DOI:10.1016/j.rser.2012.02.028 |
[2] |
NAJJAR Y S H. Hydrogen safety:The road toward green technology[J]. International Journal of Hydrogen Energy, 2013, 38(25): 10716-10728. DOI:10.1016/j.ijhydene.2013.05.126 |
[3] |
SAFFERS J B, MOLKOV V V. Hydrogen safety engineering framework and elementary design safety tools[J]. International Journal of Hydrogen Energy, 2014, 39(11): 6268-6285. DOI:10.1016/j.ijhydene.2013.06.060 |
[4] |
RUGGLES A J, EKOTO I W. Experimental investigation of nozzle aspect ratio effects on underexpanded hydrogen jet release characteristics[J]. International Journal of Hydrogen Energy, 2014, 39(35): 20331-20338. DOI:10.1016/j.ijhydene.2014.04.143 |
[5] |
李雪芳, 毕景良, 柯道友. 高压氢气储存系统泄漏的热力学模型[J]. 清华大学学报(自然科学版), 2013, 53(4): 503-508. LI X F, BI J L, CHRISTOPHER D M. Thermodynamics models of leaks from high pressure hydrogen storage systems[J]. Journal of Tsinghua University (Science and Technology), 2013, 53(4): 503-508. (in Chinese) |
[6] |
BONELLI F, VIGGIANO A, MAGI V. A numerical analysis of hydrogen underexpanded jets under real gas assumption[J]. Journal of Fluids Engineering, 2013, 135(12): 1894-1901. |
[7] |
HAN S H, CHANG D, KIM J S. Release characteristics of highly pressurized hydrogen through a small hole[J]. International Journal of Hydrogen Energy, 2013, 38(8): 3503-3512. DOI:10.1016/j.ijhydene.2012.11.071 |
[8] |
VUORINEN V, YU J, TIRUNAGARI S, et al. Large-eddy simulation of highly underexpanded transient gas jets[J]. Physics of Fluids, 2013, 25(1): 281-319. |
[9] |
XU B P, EL HIMA L, WEN J X, et al. Numerical study on the spontaneous ignition of pressurized hydrogen release through a tube into air[J]. Journal of Loss Prevention in the Process Industries, 2008, 21(2): 205-213. DOI:10.1016/j.jlp.2007.06.015 |
[10] |
RUGGLES A J, EKOTO I W. Ignitability and mixing of underexpanded hydrogen jets[J]. International Journal of Hydrogen Energy, 2012, 37(22): 17549-17560. DOI:10.1016/j.ijhydene.2012.03.063 |
[11] |
XIAO J, TRAVIS J R, BREITUNG W. Hydrogen release from a high pressure gaseous hydrogen reservoir in case of a small leak[J]. International Journal of Hydrogen Energy, 2011, 36(3): 2545-2554. DOI:10.1016/j.ijhydene.2010.05.069 |
[12] |
LI X, CHRISTOPHER D M, HECHT E S, et al. Comparison of two-layer model for hydrogen and helium jets with notional nozzle model predictions and experimental data for pressures up to 35 MPa[J]. International Journal of Hydrogen Energy, 2017, 42(11): 7457-7466. DOI:10.1016/j.ijhydene.2016.05.214 |
[13] |
LI X, HECHT E S, CHRISTOPHER D M. Validation of a reduced-order jet model for subsonic and underexpanded hydrogen jets[J]. International Journal of Hydrogen Energy, 2016, 41(2): 1348-1358. DOI:10.1016/j.ijhydene.2015.10.071 |
[14] |
李雪芳, 王俞杰, 罗峰, 等. 欠膨胀氢气射流激波结构数值模拟研究[J]. 工程热物理学报, 2018, 39(4): 880-886. LI X F, WANG Y J, LUO F, et al. Numerical simulation of shock structures of under-expanded hydrogen jets[J]. Journal of Engineering Thermophysics, 2018, 39(4): 880-886. (in Chinese) |
[15] |
WOODMANSEE M A, IYER V. Nonintrusive pressure and temperature measurements in an underexpanded sonic jet flowfield[J]. AIAA Journal, 2004, 42(6): 1170-1180. DOI:10.2514/1.10418 |
[16] |
SCHEFER R W, HOUF W G, WILLIAMS T C, et al. Characterization of high-pressure, underexpanded hydrogen-jet flames[J]. International Journal of Hydrogen Energy, 2007, 32(12): 2081-2093. DOI:10.1016/j.ijhydene.2006.08.037 |
[17] |
YUCEIL K B, OTUGEN M V. Scaling parameters for underexpanded supersonic jets[J]. Physics of Fluids, 2002, 14(12): 4206-4215. DOI:10.1063/1.1513796 |
[18] |
BIRCH A D, HUGHES D J, SWAFFIELD F. Velocity decay of high pressure jets[J]. Combustion Science and Technology, 1987, 52(1-3): 161-171. DOI:10.1080/00102208708952575 |
[19] |
EWAN B, MOODIE K. Structure and velocity measurements in underexpanded jets[J]. Combustion Science and Technology, 1986, 45(5-6): 275-288. DOI:10.1080/00102208608923857 |