基于虚拟点的反应流Navier-Stokes特征边界条件及其应用.pdf

基于虚拟点的反应流 Navier-Stokes特征边界条件及其应用赵俊琪1，钱琛庚2，王成1，孙远翔1（1.北京理工大学 爆炸科学与技术国家重点实验室，北京100081；2.中国航空工业空气动力研究院高速高雷诺数气动力航空科技重点实验室，辽宁，沈阳110034）摘 要：针对通风条件下大尺度加氢站事故模拟，从一维局部无黏特征分析出发，发展了使用虚拟点的三维反应流 Navier-Stokes 无反射特征边界条件，通过引入横向项和化学反应源项，有效消除了火焰面、亚声速流场与计算域边界相互作用时产生的非物理反射，实现了计算域边界处通风条件的无反射施加，提高了开放空间数值模拟的计算效率.通过开展 51 m51 m10 m 加氢站中长管拖车、加注机意外氢气泄漏扩散的数值模拟，给出了多种风况下可燃气云扩散结果，探究了通风条件及复杂环境对可燃气云发展规律的影响，对潜在危险区域进行了定量分析，并选取最危险氢气泄漏扩散结果开展高度非均匀气云爆炸的数值模拟，对不同设备处接收的超压和冲量进行了定量转化，完成了加氢站典型事故的风险评估.关键词：无反射边界条件；高压氢气泄漏；非均匀气云爆炸；数值模拟中图分类号：O362 文献标志码：A 文章编号：1001-0645(2024)04-0348-11DOI：10.15918/j.tbit1001-0645.2023.109Navier-Stokes Characteristic Boundary Conditions Using Ghost Cellsfor Reacting Flows and ApplicationsZHAO Junqi1，QIAN Chengeng2，WANG Cheng1，SUN Yuanxiang1（1.State Key Laboratory of Explosion Science and Technology,Beijing Institute of Technology,Beijing 100081,China;2.Aviation Key Laboratory of Science and Technology on Aerodynamics of High Speed and High Reynolds Number,AVIC Aerodynamics Research Institute,Shenyang,Liaoning 110034,China）Abstract：The large-scale hydrogen station accident was simulated under ventilation conditions.Based on one-dimensional local inviscid characteristic analysis,a three-dimensional Navier-Stokes non-reflection characterist-ic boundary condition for reacting flow using ghost cells was developed.By introducing transverse term andchemical reaction source term,non-physical reflection generated by flame and subsonic flow interacting withboundary was effectively eliminated,which realized the non-reflection application of ventilation conditions at thecomputational boundary and improved the computation efficiency of open space numerical simulation.Numeric-al simulation of accidental hydrogen leakage and diffusion of long tube trailers and filling machines of 51 m51 m10 m hydrogen station was carried out,and the distribution result of flammable gas cloud under variousventilation conditions was reported.The impact of complex environments on the development of flammable gasclouds was quantitatively analyzed.The most dangerous hydrogen leakage and diffusion result was selected tocarry out numerical simulation of highly inhomogeneous gas cloud explosions,the overpressure and impulse re-ceived at different places of the equipment was analyzed,and the risk assessment of typical accidents in hydro-收稿日期：2023 05 12基金项目：爆炸科学与技术国家重点实验室探索性资助项目（YBKT22-05）作者简介：赵俊琪（1998），女，硕士，E-mail：junqizhao_.通信作者：王成（1972），男，博士，教授，E-mail：.第 44 卷第 4 期北 京 理 工 大 学 学 报Vol.44No.42024 年 4 月Transactions of Beijing Institute of TechnologyApr.2024gen refueling stations was completed.Key words：non-reflecting boundary conditions；high-pressure hydrogen release；inhomogeneous gas cloud ex-plosion；numerical simulation 氢能是一种清洁、制取来源丰富、燃烧热值高的二次能源，但是氢气的物理化学性质决定了氢气在制取、运输、储存、加注和使用等过程均存在泄漏和爆炸的风险，使得公众对于氢能安全问题充满疑虑，研究加氢站内的氢能安全事故，包括高压氢气泄漏扩散、气云爆炸等，对于制定加氢站相关安全规范和标准、进行风险评估及事故调查具有重要意义1.高压氢气的危险性较高，实验开展较为困难、可重复性差，已有的研究成果大多局限于较低的压力或者在小尺度空间内2 4.而理论模型，如用于预测气体扩散的 Gauss 烟羽模型、Sutton 模型等及用于评估爆炸破坏作用的 TNT 当量模型、TNO 多能模型等，由于假设较多，仅能对简单的事故过程进行再现，无法给出详细的流场信息，也不能直观地展现动态变化.因此对于加氢站氢能安全事故的大尺度空间、时间研究，基于计算流体力学（computational fluid dynam-ics,CFD）的数值模拟成为一种主流趋势.目前常用的 CFD 商业软件包括 Fluent、FLACS、CFX 等，国内在工程尺度上有关分析事故发展规律、预测事故结果的研究几乎基于以上软件展开.但是 CFD 软件也存在缺点，国外商业软件开发较早，数值方法精度最高为二阶，依赖国外商业软件难以实现大规模问题的高精度模拟仿真，因此自主研发复杂环境下可燃气体泄漏扩散爆炸大规模高精度仿真软件具有重要的工程应用价值.对于自然通风情况下，加氢站高压氢气泄漏扩散及气云爆炸的数值模拟，需要在计算域边界施加风，即亚声速入流边界条件，同时计算域边界可能会接收到火焰面或多组分气流，即亚声速出流边界条件.边界条件通常被视为一个非关键的部分而被简化处理，例如采用零阶、一阶外推方法.由于忽略了亚声速流动的迎风性问题，这种数值边界条件通常不能等效于真实的物理边界条件，导致高精度格式在大尺度、长时间模拟中产生的数值波与边界相处作用而产生非物理振荡5，边界处的非物理反射还将向上游传播及其他的边界产生强数值耦合.最初POINSOT 等5提 出 了 N-S 特 征 边 界 条 件（Navier-Stokes characteristic boundary conditions,NSCBC），通过对一维特征波进行分解，对进入、离开的特征波进行控制，可以显著改善亚声速入流、出流在边界处的反射问题.后续，BAUM 等6通过将 POINSOT 等5开发的一维特征波动关系的边界公式扩展到用真实热力学模型及详细化学反应模型描述的真实气体混合物，给出了一种基于湍流反应流的准确、稳定的NSCBC.目前广泛应用的 NSCBC 应用研究多数考虑变量在边界上的时间推进，以计算边界处的通量.例如，KIM 等7提出了在广义坐标系下高阶格式的NSCBC 扩展保守形式，SUTHERLAND 等8主要研究了理想多组分气体流动的边界条件，DAVILLER 等9讨论了入射特征波波幅的指定方法.在本文计算背景中，数值解设置在网格单元中心，计算域边界位于两个节点之间，如图 1 所示，若要计算流过边界的数值通量，需要通过计算域内节点值以及域外虚拟点值进行推导表示，无法直接使用传统的 NSCBC.目前关于使用虚拟点的 NSCBC 研究较少，在研究中也未全面考虑附加项.例如，GROSS 等10针对单组份气体推导了一种使用虚拟点的 NSCBC，通过描述虚拟点内的流动变量，使其产生的通量梯度满足非反射条件，PIROZZOLI 等11建立了可压缩流动的 NSCBC，通过引入计算域外第一个虚拟点的流动信息参考值，使用一阶偏心差分近似估计入射波特征波波幅，MOTHEAU 等12针对无化学反应流动的亚声速流动问题提出了一种基于 NSCBC 给虚拟点赋值的方法.边界虚拟点I,JI1,JI+1,JI+2,JFI+1/2,J图 1 通量离散方法Fig.1 Discretization procedure of the flux 因此，本文针对多组分反应流，从特征变量的一维局部无黏方程出发，推导多组分反应流 Navier-Stokes 方程的特征波幅表达式，引入基于特征空间的反应源项、横向项，分析边界处特征波的运动方向，使用计算域内信息计算离开计算域边界的特征波幅，第 4 期赵俊琪等：基于虚拟点的反应流 Navier-Stokes 特征边界条件及其应用349使用目标值松弛法计算进入计算域边界的特征波幅.通过边界处特征波幅计算原始变量梯度，结合高精度插值计算虚拟点处的守恒变量，发展基于虚拟点的反应流 N-S 特征无反射边界条件方法.通过该程序可以更高效、精确、稳定地实现通风条件下大尺度加氢站氢气泄漏扩散爆炸全过程的数值模拟，以给出环境风对可燃气云形成的影响规律以及可能出现的气云爆炸事故后果.1 基于虚拟点的 Navier-Stokes 特征无反射边界条件方法N-S 特征边界条件的核心思想是，在边界处不进行强赋值以求解守恒变量，而是根据特征波的方向控制波的传播，以更缓和的方式求解.在特征分析的过程中，为可以同时应用于非反应流和反应流，引入了基于特征空间的化学反应源项，为解决高维问题的模拟，引入了基于特征空间的横向项.本文程序基于单元中心型均匀网格的有限差分方法展开，使用的插值模板为 5 点模板，为求解计算域边界通量梯度，需要在计算域外引入 3 个虚拟点，而虚拟点的赋值方法与真实物理域中的流场信息有关，对于亚声速入流、出流边界应用 N-S 特征边界条件方法给虚拟点赋值.1.1 控制方程NsNrNe=5+Ns1对于包含个组分，可能会产生个化学反应的气体混合物，在笛卡尔坐标系下，其非定常 3D 可压缩 NavierStokes 控制方程组含有个方程，表述如下Ut+Fv(U)=Fv(U)+S（1）U=uEYiT（2）F(U)=u(uu+p)(E+p)uYiuT（3）Fv(U)=0(uq)Vi,YiT（4）S=000 iT（5）U式中：为守恒向量，与其对应的原始变量向量为Q=uvwpYiT;（6）F(U)Fv(U)S为对流通量项；为黏性通量项；为化学反应源项.其中化学反应速率的计算考虑一步化学vxyzNs1u=(u,v,w)uuEYiipqVi,i ii反应，采用 Arrhenius 公式得到；下标 表示黏性/扩散项；下标、为、3 个方向的空间索引；下标 i为组分索引，此处对应 1 到个组分；为克罗内 克 符 号.在 控 制 方 程 中，为 流 体 质 量 密 度；为速度向量；、为速度分量；为比总能量；为组分 的质量分数；为压力，为黏性应力张量；为热通量分量；为组分 的质量扩散速度分量；为组分 的化学反应生成速率.本文基于多组分反应流 Navier-Stokes 控制方程组，采用大涡模拟，使用 5 阶 WENO 格式离散对流通量项、守恒型 4 阶中心差分格式离散扩散通量项，使用 3 阶 RK TVD 进行时间推进，给出高精度有限差分程序用于模拟.1.2 数值方法Tx对于 N-S 控制方程组，仅考虑一维问题并忽略扩散通量项，并将另外 2 个方向上的空间偏导项合并简写为横向项.以下为 方向上的推导，另外 2个方向的结果完全对称.Ut+F(U)x=T+S（7）xJB上式为基于守恒变量的一维控制方程组，给出方向上守恒变量和原始变量之间的关系矩阵以及对流通量项和原始变量之间的关系矩阵，它们之间关系为F(U)x=FUx=FJ1Qx=BxQx（8）因此可改写出基于原始变量的一维控制方程为J1Qt+BxQx=T+S（9）Qt+FQx=TQ+SQ（10）FF=VxxV1xxVxV1xLx将非守恒雅可比矩阵进行特征分解，得到特征值的对角矩阵以及对应的左特征矩阵、右特征矩阵.其中，特征值矩阵中的元素特征值在物理上表示特征波的传播速度.定义特征波波幅为Lx=xV1xQx（11）若忽略黏性项和源项，并以特征形式来重新描述控制方程，附加横向项和反应源项为V1xQt=LxTC,x+SC,x（12）其中，基于特征变量的横向项表示为350北 京 理 工 大 学 学 报第 44 卷TC,x=V1xTQ=12c2v(pycuy)+pvy+12c2w(pzcuz)+pwz1c2v(c2ypy)+1c2w(c2zpz)(vvy+py)+(wvz)(vwy)+(wwz+pz)12c2v(py+cuy)+pvy+12c2w(pz+cuz)+pwzvc2Yi(c2ypy)+pYiy+wc2Yi(c2zpz)+pYiz（13）spsYi由于忽略体积力，只有对应压力和组分的源项不为 0，并分别表示为、.基于特征变量的反应源项表示为SC,x=V1xSQ=12sp1c2sp0 012spsYiT（14）LL的简单物理含义可以描述为，除去基于特征变量的横向项、反应源项，为特征变量随时间变化的相反数.为无化学反应的 N-S 方程组所计算出的特征波波幅，由 POINSOT 等5给出 Lx=V1xQx=12(pxcux)2c2(c2xpx)3vx4wx52(px+cux)5+kYixT（15）LL将横向项、反应源项扩展至特征波幅，称为，因L此可表示出边界处原始变量的空间梯度与的关系为Qx=F1xdx=Vx1Lx=L22+1c2(L11+L55)1c(L55L11)L33L44L11+L55L5+i5+iT（16）xx1/2n=1 sgn=1xxnx+1/2n=nxsgn=1方向左侧边界有最左侧节点，方向右侧边界有最右侧节点，.通过梯度插值可得到虚拟点的值为Qnsgn=Qn+sgnsgn2xQx,Qn2sgn=2Qn+sgn3Qn+6Qnsgn+sgn6xQx,Qn3sgn=3Qn+sgn+10Qn18Qnsgn+6Qn2sgnsgn12xQx（17）x1xnx计算左、右边界外虚拟点的值时，分别对应使用计算域内、处的原始变量梯度.1.3 特征波幅 1.3.1 亚音速无反射出流边界条件pppxLsgn=1sgn=14+Ns1Lp p指定无穷远处的静压，当计算域流出边界的局部压力与不同时，仅部分无反射并将出口平均压力松弛到.考虑一个三维笛卡尔坐标系下 方向边界的流出，对于该边界，(当 n=1 时=5、，当 n=nx时=1、)进入计算域，而其他条特征波动离开计算域.对于，施加了一个物理边界条件，无穷远处的静压，将波动振幅与压差关联，并扩展包括源项、横向项为L=12sp+K(p p)(1)T（18）K式中：为压力松弛系数；为横向项贡献系数.根据 RUDY 和 STRIKWERDA13的建议，压力松弛系数可以表示为K=c(1Ma2)/2l（19）lc1 式中：Ma 为边界处流动的平均马赫数；为计算域的特征尺寸，在本文中取为网格尺寸；为平均当地声速；为人工系数，与 的取值一样高度依赖于流动特性.当较小时，平均压力可能出现漂移，而越大时，边界上的反射强度会更大.针对不同的计算情况，应该取不同的值，主要依靠经验选取.SELLE 等14提出，为避免平均压力漂移或较大反射，系数范围可以取为.1.3.2 亚音速无反射入流边界条件4+Ns1Lut,vt,wtTtYti对于亚音速无反射入流，有条特征波进入计算域，仅有一条特征波离开计算域.直接施加目标入流速度可能会导致声波反射并与入流边界相互作用，为了在边界上维持想要施加的目标速度（）、温度（）和质量分数（），可以在入流边界处通过对目标值进行松弛，来控制特征波幅并反映部分非反射的特性,第 4 期赵俊琪等：基于虚拟点的反应流 Navier-Stokes 特征边界条件及其应用351L=12sp+(uut)(1)T（20）L2=1c2sp+2(T Tt)(12)T2（21）L3=3(vvt)(13)T3（22）L4=4(wwt)(14)T4（23）L5+i=sYi+5+i(YiYti)(15+i)T5+i（24）m(m=1,2,5+Ns1)m式中：松弛系数用于控制施加到计算域的波动振幅阻尼；用于控制横向项的贡献量.1.4 典型算例验证 1.4.1 非反应流xL=13x=y=0.1Tref=300pref=101.325该案例为在 方向均匀流场上叠加的单个涡旋，为评估边界条件方法性能的常用案例，此外在非反应流的案例中不引入基于特征空间的化学反应源项.正方形计算域长度为 mm，网格尺寸为 mm，初始时为常温 K、常压 kPa 下的氮气.初始流函数的定义式为=exp(xx0)2+(yy0)22R2v)（25）(x0,y0)=0.11m2/sRv=0.1L=1.3式中：为涡旋中心；和 mm 分别为涡旋强度和半径.因此速度场定义为u=u0+yv=x（26）u0=100 m/s式中：为初始流速.初始压力场表示为p=prefexp2(cRv)2exp(xx0)2+(yy0)2R2v)（27）=0.25=0.34 个计算域边界均分别使用一阶外推方法、亚声速无反射出流边界条件方法，无反射边界条件中参数，横向项参数为.图 2 明显表现出来 2 种方法的差异性，当涡旋开始与边界相互作用时，一阶外推方法使涡旋在计算域边界变形，通过边界的过程中涡旋变形扭曲更加严重，开始发散，并在边界处产生强烈的压力反射返回计算域，涡旋在离开计算域时非物理反射依然存在于计算域中.而使用虚拟点的 NSCBC 使涡旋在离开计算域边界过程中几乎没有变形，也达到了几何没有反射波回到计算域的效果.可见亚声速无反射边界条件可以良好地抑制高维效应带来的压力、声波反射.1210864y/mmx/mm201210864y/mm201210864y/mm201210864y/mm201210864y/mm201210864y/mm201210864y/mm201210864y/mm2002468 10 12x/mm02468 10 12x/mm02468 10 12x/mm02468 10 12x/mm02468 10 12x/mm02468 10 12x/mm02468 10 12x/mm02468 10 12(a)t=4 s,一阶外推方法(b)t=26 s,一阶外推方法(c)t=40 s,一阶外推方法(d)t=54 s,一阶外推方法(h)t=54 s,使用虚拟点的NSCBC(g)t=40 s,使用虚拟点的NSCBC(f)t=26 s,使用虚拟点的NSCBC(e)t=4 s,使用虚拟点的NSCBCmin(p*)=0.00max(p*)=0.98min(p*)=0.56max(p*)=1.05min(p*)=0.88max(p*)=1.11min(p*)=0.65max(p*)=0.70min(p*)=0.03max(p*)=0.03min(p*)=0.06max(p*)=0.63min(p*)=0.08max(p*)=0.82min(p*)=0.00max(p*)=0.98图 2 t=4,26,40,54 s 速度等值线和归一化压力场对比结果Fig.2 Velocity contours and normalized pressure field at t=4,26,40,54 s 1.4.2 反应流L=3.0 x=y=12.5m该案例计算域长度为 mm，网格尺寸为.初始时计算域有球心（1.7 mm,1.7 mm），半径 1 mm 的点火源，为当量比为 1 的氢气空气混合物的完全反应产物，混合物温度为 2 250 K，计算=0.25域其他部分均为均匀预混的当量比为 1 的氢气空气混合物，混合物温度为 300 K，全场压力为 101.325 kPa.3 个方向上每个边界均分别使用一阶外推方法、亚声速无反射出流边界条件方法，在亚声速无反射出流边界条件中参数.352北 京 理 工 大 学 学 报第 44 卷图 3 显示了二维预混火焰流出计算域不同时刻温度等值面对比结果.非反射的边界条件方法保持稳定，而一阶外推给出了不符合真实火焰流动的结果.在 2 s 时，火焰还未到达计算域边界，两种方法给出的结果一致；在 7 s 时，火焰面接近计算域边界，一阶外推对火焰具有拉伸效果，而无反射条件对于火焰传播无影响；在 15 s 时，一阶外推的温度场受到边界的反射而产生反转；在 28 s 时，无反射条件的火焰面已经接触到右侧和上侧的计算域边界，温度场结构保持稳定不变，而一阶外推导致温度场被破坏，产生较大的逆温度梯度，气流在边界的流动受到阻滞，产生了非物理的结果.332y/mmx/mm2113x/mm213x/mm213x/mm213x/mm213x/mm213x/mm213x/mm21032y/mm1032y/mm1032y/mm1032y/mm1032y/mm1032y/mm1032y/mm10(h)t=28 s,使用虚拟点的NSCBC2 400T/K2 0001 6001 2008004002 400T/K2 0001 6001 200800400(g)t=15 s,使用虚拟点的NSCBC(f)t=7 s,使用虚拟点的NSCBC(e)t=2 s,使用虚拟点的NSCBC(a)t=2 s,一阶外推方法(b)t=7 s,一阶外推方法(c)t=15 s,一阶外推方法(d)t=28 s,一阶外推方法图 3 t=2,7,15,28 s 温度等值线对比结果Fig.3 Temperature contours at t=2,7,15,28 s 2 加氢站氢气泄漏扩散爆炸全链条数值模拟 2.1 计算模型对某真实加氢站进行合理简化，建立几何模型.以加氢站的西南角为原点，建立三维直角坐标系，向东为 X 轴正方向，向北为 Y 轴正方向，垂直向上为Z 轴正方向，将计算域取为 51 m51 m10 m，几何模型如图 4 所示.选择发生频率更高的小规模连续泄漏事件作为数值模拟分析的对象，设定可能发生高压氢气泄漏的地方包括长管拖车 A 末端面板的泄气软管和加注机 B 的加注软管.其中，均考虑配件或管道的 100%损坏破裂.表 1 为数值模拟考虑的泄漏场景15 16.对于高压氢气的泄漏扩散，除了内在因素如泄漏孔径、泄漏压力及泄漏高度等，外部因素也对氢气可燃气云的形成有较大影响.其中，风向会影响泄漏扩散的主要方向，风速会影响扩散速度，因此泄漏气体的扩散过程与环境风况密切相关.参考当地的真实气候条件17，考虑采用平均气温、气压，将模拟环境参数列于表 2.AB图 4 加氢站几何模型Fig.4 Simplified geometric model of hydrogen station 表 1 泄漏场景参数Tab.1 Parameters of release scenario临氢设备编号泄漏原因泄漏压力/MPa泄漏孔径/mm最大泄漏量/kg事故频率/aA卸气软管破裂162090105B加注软管破裂43.820904106第 4 期赵俊琪等：基于虚拟点的反应流 Navier-Stokes 特征边界条件及其应用353 表 2 环境参数Tab.2 Parameters of environment conditions环境编号风向风速/（ms1）气温/K气压/kPa1无风0.0293.15100.02+Y5.0293.15100.03+Y12.0293.15100.04X5.0293.15100.05X12.0293.15100.0 假设所有设备以及地面在事故过程中不发生破损，采用绝热无滑移固壁边界条件.通风边界采用NSCBC 亚声速入流边界条件，其他方向的计算域边界均采用 NSCBC 亚声速出流边界条件.将长管拖车的泄漏方向考虑为水平向西朝向附近的防火墙，泄漏位置为(7.4,25.5,1.6)m，将加注机的泄漏方向考虑为垂直向上朝向罩棚，泄漏位置为(29.1,32.5,2.5)m.由于不同临氢设备 A、B 的泄漏口大小不同，分别建立不同网格，采用 Molkov 虚喷管模型对高压氢气射流进行简化模拟18，在泄漏口处给出等效的入流条件，因此计算域网格尺寸放大至分别为 90、140 mm，计算域网格总量分别为 3 568 万、947.8 万，采用 256核并行计算.2.2 高压氢气泄漏扩散高压氢气泄漏扩散的气流运动十分复杂，在泄漏口附近，气流主要受动量主导，沿泄漏口轴线方向氢气浓度高速上升，而远场区域，周围障碍物形状以及通风条件等多种环境因素造成的湍流对氢气扩散影响更大19.环境风对可燃气云形成可以产生两种截然不同的影响，既可以削弱自然扩散、抑制氢气与空气的混合并缩短扩散距离，也可以推动可燃气云流动、增强湍流、稀释氢气并增大扩散距离.因此，泄漏形成的高速气流在障碍物和环境风共同作用下形成许多小尺度涡旋，不仅具有横向脉动，还存在与气流主方向相反的运动方向，使气流动量、质量传递速率急速上升.给出 10 种工况下可燃气云危险距离、可燃气云体积的变化趋势如图 5 所示，其中考虑人体高度，危险距离为 1.5 m 高度水平方向估计.可看出水平方向危险距离总体满足随风速增大而增大的趋势，但 B2在低风速时危险距离反而比无风时更低，这是环境风对自然扩散的抑制造成的.对于泄漏口 A/B，不同风向时可燃气云体积的变化趋势相反，可见风向对于环境风稀释可燃气云的效果起到极大影响.302520151050510风速/(ms1)危险距离/m气云体积/m3150510风速/(ms1)15A1A2/A3A4/A5B1B2/B3B4/B5A1A2/A3A4/A5B1B2/B3B4/B51 0008006004002000(a)水平方向危险距离(b)可燃气云体积图 5 t=8 s 时所有工况可燃气云 z=1.5 m 的水平方向危险距离及气云体积对比Fig.5 Horizontal dangerous distance at z=1.5 m and gas cloud volume of all simulation conditions at t=8 s 选取危险距离及可燃气云体积综合最大的工况B5进行分析，图 6 为不同时刻的可燃气云浓度云图.泄漏初期，射流由动量控制，可燃气云分布不受环境风影响，竖直向上的射流与罩棚顶部相互作用，射流速度方向改变.泄漏后期，可燃气云在罩棚下大量聚集，随着风速增加，可燃气云扩散路径被改变，Y 方向上的扩散距离小幅度减少，X 方向上的扩散距离大幅度增加.东风使可燃气云扩散到罩棚外，倾向从压缩机和长管拖车之间的通道运动.YH2z=1.5y=32.5图 7 为不同时刻氢气体积分数的水平切片（m）及垂直切片（m）.由于风速高，且下风向存在复杂障碍物，从流线可看出各个障碍物354北 京 理 工 大 学 学 报第 44 卷之间存在大气扰动并产生涡旋.风场使得气流基本朝向X 方向流动，受到障碍物的影响，主要分布于障碍物之间的通道中，造成了压缩机附近大量氢气聚集，因此潜在危险更大.后续因为泄漏氢气总量不断增加，罩棚下的可燃气云浓度升高，可燃气云在空中被延长，由于湍流扩散能力强，浓度较低呈现出不连续的碎片状，最高到达 7.18 m.总的来说，环境风速对泄漏扩散的影响较为复杂.更高的风速将抑制可燃气云在某些方向的自然扩散，如垂直风向的水平扩散以及竖直向上的运动等.但是风速越高，大气湍流越强，空气的稀释能力增强，环境风对可燃气云的推动作用也越强，可燃气云加速扩散，在下风向同一位置处的氢气浓度越低.由于这两种作用相互矛盾，导致风速对可燃气云扩散分布的影响规律呈现非线性，在较高风速下第二种作用会占主导因素，而其余情况下则需要结合风向、障碍物分布综合考虑.2.3 非均匀气云爆炸相对于常见的可燃气云爆炸研究，本文由真实场景下泄漏扩散模拟得到的可燃气云由于复杂障碍物分布和环境风影响，爆炸极限范围内的气云在空间上分布高度不连续、不均匀，在泄漏口附近空气夹带较少，不容易燃烧，在远离泄漏口的地方存在复杂湍流流动以及不均匀的初始速度分布，因此其爆炸结果无法用常见的经验规律来预测.使用上述数值方法开展加氢站内大规模、非规则、非均匀气云爆炸模拟.以 B5的可燃气云分布为初始条件，设置弱点火源位置为(22.0,31.0,1.5)m，并设置一系列压力监(a)t=2 s(b)t=4 s(c)t=6 s(b)t=8 s图 6 不同时刻 B5可燃气云轮廓Fig.6 Flammable gas cloud of B5 at different time 0.1 0.13H2体积分数0.265 0.4 0.535 0.67 0.8050.1 0.13H2体积分数0.265 0.4 0.535 0.67 0.8050.1 0.13H2体积分数0.265 0.4 0.535 0.67 0.8050.1 0.13H2体积分数0.265 0.4 0.535 0.67 0.80555010101520 x/mz/m5010z/m2530355101520 x/m2530355101520 x/my/m25303550454035302520(a)t=4 s水平方向切片(c)t=4 s垂直方向切片(d)t=8 s垂直方向切片5101520 x/m253035y/m50454035302520(b)t=8 s水平方向切片图 7 不同时刻可燃气云 z=1.5 m 水平方向切片、x=29.1 m 垂直方向切片的体积分数分布对比Fig.7 Volume fractions contour of horizontal slice at z=1.5 m,vertical slice at x=29.1 m at different time第 4 期赵俊琪等：基于虚拟点的反应流 Navier-Stokes 特征边界条件及其应用355测点(21.2,26.4,1.5)、(39.6,41,1.5)、(25.2,52.8,1.5)、(29.1,32.6,5.5)m 用于分析各设备所接收到的超压.t=40 ms 时压力波已传播至计算域边界，图 8 为该时刻 B5工况下气云爆炸模拟结果.工况 B5在风向上存在较多障碍物，可燃气云在环境风与复杂障碍物群的耦合作用下产生了剧烈的湍流，且障碍物附近存在大气涡旋，使得火焰面形状变得不规则，从而增大火焰面积，导致火焰传播速度急剧上升并进一步促进湍流的发展，同时，当量浓度区域由于化学反应更加剧烈，可在切片中看到相对高温的“带状”区域20 21.在 40 ms 时，火焰的高温辐射范围不再扩展，但火焰面褶皱较多.气云的尺寸对爆炸强度有显著影响，气云初始尺寸越大，火焰加速距离越长，对于前方未燃气体的压缩程度增加，最终正向反馈激励了超压上升速率，在爆心处可产生约 600 kPa 的超压.40 ms 之后，压力波传播至控制室，防火墙的存在显著减弱了压力波的危害，有防火墙遮蔽的控制室前比无防火墙时的最大超压低约 70 kPa.考虑使人重伤的温度为 450 K，则火焰产生的危险范围约为 29 m16 m.考虑使人轻伤、重伤、死亡的超压为 20、50、400 kPa，则压力波产生的轻伤、重伤和死亡危险范围约为 51 m51 m、48 m46 m、14 m8 m.1200温度/K500800 1 100 1 400 1 700 2 000 2 3001200超压/kPa500800 1 100 1 400 1 700 2 000 2 300温度/K100.0 15.0 37.590.0142.5 195.0 247.5 300.0超压/kPa100.0 15.0 37.590.0142.5 195.0 247.5 300.00102030 x/my/m405050403020100102030 x/m4050y/m5040302010(a)温度云图(b)z=1.5 m水平方向温度分布(c)超压云图(d)z=1.5 m水平方向超压分布图 8 t=40 ms 时气云爆炸火焰面、压力波阵面Fig.8 Flame and pressure wave front of gas cloud explosion at t=40 ms 为精确分析压力波的危害，根据超压冲量破坏准则对各设施进行危害评估.图 9 为工况 B5测点P1P4的超压时程曲线.考虑到 B5总体破坏范围更大，选取了距离爆心更远的监测点，而这组监测点满足随爆心越远爆炸超压越低、超压变化速率更慢的趋势.同一高度处，P1最先接收到超压，由于处于障碍物之间，超压峰值远高于其他监测点，达到 620 kPa，第一次正压作用时长为 17.9 ms，冲量是第二次正压356北 京 理 工 大 学 学 报第 44 卷较弱，主要破坏为一次正压带来.由于 P2、P3和爆距相差不大，而又距离爆心较远，气云初始不规则分布带来的差别开始减弱，超压曲线相似，超压峰值约为60 kPa，正压作用时间分别为 27.1 ms 和 15 ms.P4位于泄漏口正上方，罩棚的约束导致了此处可燃气云大量聚集，火焰传播到该处后超压约为 200 kPa，正压持续了约 22 ms，将在更短时间内造成更大范围内的严重破坏.根据超压冲量准则22，监测点 P1P4对应的设备破坏等级均为 6，对于加氢站有着毁灭性的全场破坏.800600400200002040时间/ms超压/kPa60P1P2P3P4图 9 工况 B5监测点 P1P4超压时程曲线Fig.9 Evolution of overpressure at P1P4 of B5 3 结论本文基于高精度有限差分程序，发展了基于虚拟点的无反射特征边界条件方法，实现了大尺度加氢站氢气泄漏点火爆炸完整过程的数值模拟，并分析了可燃气云形成以及气云爆炸的影响因素，完成了加氢站氢气爆炸事故的定量危害评估.具体结论如下：传统的 Navier-Stokes 特征边界条件方法考虑边界上的时间梯度以求解边界处通量，不适用于使用虚拟点的中心型网格，针对此问题，将一维无黏特征方程在边界处离散，获得并使用边界处的空间梯度对虚拟点进行赋值，构造了使用虚拟点的 N-S特征边界条件方法.通