汽油作为一种工业制成品，广泛应用于交通运输和工业清洗行业。但是，作为一种易燃易爆易挥发的危险品，在储存、运输及加注的全过程中，汽油时刻有着安全威胁。为有效应对安全威胁，降低人类的生命和财产损失，非常有必要针对汽油蒸气/空气混合气的爆炸特征开展深入研究。此前，已有大量针对不同可燃气体爆炸特征的研究，发现受限空间结构和形状特征对可燃气爆炸超压和火焰的变化规律有较大影响。主要影响因素包括为狭长受限空间长径比^[1-2]、半受限空间开口尺寸及位置^[3-5]、受限空间障碍物形状及尺寸^[6-8]和所含弯管及分支结构特征^[9-11]等。特别地，弯管及分支结构是储油洞库坑道及输油管路中常见的结构形式，研究表明弯管及分支结构的存在对爆炸超压和火焰传播速度有较强的激励作用，并在一定程度上可与挡板式障碍物类比^[9]。Zhu等^[12]发现含弯管结构受限空间可燃气爆炸规律主要受弯管结构数量影响，而受弯管方向影响较小。Niu等^[13]对甲烷在输气管网中的爆炸特征进行了研究，发现爆炸超压峰值及火焰传播速度受连接平行管网的分支长度影响，两分支长度不同时爆炸更为剧烈复杂。

尽管弯管和分支结构在形状特征上有一定相似，但由于分支结构具有更复杂的通道和多个流动终点^[14]，因此确有必要对含多分支结构受限空间油气爆炸特征进行研究。Zhang等^[14]对5种不同分支结构受限空间爆炸进行研究，发现油气爆燃指数的增加与分支数量呈函数关系。Li等^[15]、杜扬等^[16]通过实验研究，发现分支结构可以增大汽油蒸气/空气混合气爆炸超压峰值、升压速率及火焰传播速度。杜扬等^[17]对含双侧分支受限空间油气爆炸进行研究，发现火焰在含有双侧分支结构的管道内会呈现“浪花状”火焰。目前，人们已经对含分支结构受限空间油气爆炸进行了一定的研究，但研究多集中于实验分析和现象描述，对爆炸流场和火焰的耦合关系仍需深入研究。

本文中，在已有实验研究^[17]的基础上，通过大涡模拟方法对含双侧分支受限空间油气爆炸过程进行模拟，并将火焰形态、传播速度和爆炸动态超压等模拟结果与实验进行验证，分析含双侧分支受限空间油气爆炸流场结特征、火焰传播行为和动态超压特性，指出“浪花状”火焰形成的内在原因，揭示流场、火焰和超压的耦合机理。

1.   数学模型

1.1   大涡模拟控制方程

大涡模拟通过对小尺度湍流脉动进行过滤，求解计算量比直接模拟方法小，但空间分辨率比雷诺平均方法高，适合于捕捉油气爆炸的湍流特征^[18]。

大涡模拟对小尺度湍流脉动进行Favre过滤后，利用N-S方程对大尺度湍流脉动进行直接求解。经Favre过滤后，大涡模拟控制方程的连续性方程、动量方程、能量方程、反应进程变量及其守恒方程为：

式中：ρ为密度，t为时间，u_i、u_j为速度分量，x_i、x_j为方向分量，τ_ij为黏性应力张量分量，p为压力，h_s为显焓，λ为热导率，T为温度，D为扩散系数，${\dot \omega _{\rm{c}}}$为归一化化学反应速率。Y_f为局部未燃燃料质量分数，Y_f⁰为初始未燃混合气中燃料质量分数；c为反应进程变量，c=0代表化学反应未进行，c=1代表化学反应完全；符号～表示Favre过滤量，—表示物理空间过滤量。

为使大涡模拟控制方程封闭，需建立亚网格模型。壁面自适应局部涡黏模型（WALE）能够较好地捕捉流场从层流到湍流的转变，且计算量较小^[19]。因此，采用WALE模型作为计算的亚网格模型。涡黏模型方程为：

式中：δ_ij为单位张量分量；μ_S为亚网格黏性系数，为：

式中：L_S为亚网格尺度混合长度，${S_{ij}}$为旋率张量分量。

1.2   亚网格燃烧模型

湍流燃烧过程中，火焰锋面厚度较小，为取得精确的结果，必须采用更精细的网格。但这样会使网格量增大，运算成本陡增。因此，采用Zimont模型^[20]，对火焰锋面进行增厚处理，确保在有限运算量下精确模拟火焰锋面。在Zimont模型中，湍流火焰速度公式为：

式中：A^*为模型常数，计算中取值0.5^[20]；${u'_\Delta }$为亚网格脉动速度；S_L为层流火焰速度，计算中取0.43 m/s^[6]；$\chi$为未燃混合物的摩尔传热系数；$L_\Delta ^{}$为亚网格湍流长度尺度。

2.   数值模拟

2.1   物理模型

以文献[6]中的实验数据为基础开展数值模拟。其实验系统如图1所示，由实验管道、点火系统、压力采集系统、高速摄影系统、油气产生系统、油气浓度测试系统、同步触发系统和数据采集系统等部分组成。图2为实验管道尺寸，数值模拟几何模型与实验管道尺寸一致。主管道及两侧分支管道内部横截面尺寸均为100 mm×100 mm，主管道左端面密闭，中心处为点火点，右端面开口率为50%。实验前，用聚乙烯薄膜密封右端开口，但爆炸破膜压力很小，对爆炸发展演变影响较小，因此在数值计算过程中未考虑薄膜的影响^[8]。为保证计算结果的准确性，数值计算中在管道右端开口处设置了一个700 mm×500 mm×800 mm的外部空间模拟实验外场，以抑制爆炸回流对管道内部的影响，如图2所示。沿爆炸传播方向，在分支前、后各275 mm的主管道处分别设置两个压阻型传感器P1、P2，对应数值计算相同位置两个压力监测点。实验时，管道充装汽油体积分数1.73%的汽油-空气混合气，环境压力约为1.013×10⁵ Pa，环境温度约为300 K。高速摄影仪记录爆炸过程火焰发展过程，帧率为500 Hz。

图  1  实验系统

Figure  1.  Experimental system

下载: 全尺寸图片 幻灯片

图  2  实验管道

Figure  2.  Experimental pipe

下载: 全尺寸图片 幻灯片

2.2   网格划分

大涡模拟计算过程中，对较大尺寸涡团直接求解、对较小尺寸涡团过滤的特性决定了利用大涡模拟进行数值计算没有“网格无关性”，越小的网格尺寸和时间步长会获得越精确的计算结果，这与雷诺平均方法不同。但随着网格尺寸趋近更小，计算量将趋近于直接模拟方法，因此在网格划分上要尽可能平衡计算精度和计算成本^[21]。经过调试验证后，本文中网格划分如图3所示，图3(a)为计算域整体网格划分，图3(b)为分支处局部加密网格。因计算区域规则，内场和外场网格均为六面体结构网格。在内场分支结构区域加密网格，最小为2 mm；外场沿爆炸传播方向网格逐渐增大，最大为10 mm：计算域网格总数约为490万个。

图  3  网格划分

Figure  3.  Computational grids

下载: 全尺寸图片 幻灯片

2.3   边界条件和初始条件

为与实验条件相对应，对含双侧分支结构受限空间充装汽油体积分数1.73%的汽油-空气混合气的爆炸过程进行模拟。主管道和分支管道壁面设置为绝热无滑移壁面且温度为300 K，外场空间壁面设置为压力出口且压力为1.013×10⁵ Pa。环境温度为300 K，流场初始速度和反应进程变量为零。在实验管道左端面中心处设置一个半径5 mm的半球形区域，将其内部反应进程变量设为1，用以模拟实验点火过程。

2.4   计算方法

计算过程中，采用二阶迎风格式对对流项进行离散，采用二阶中心差分格式对扩散项进行离散，采用SIMPLE算法求解压力和速度耦合方程组。时间步长为10⁻⁴ s，每时间步迭代20次，保证能量残差小于10⁻⁶，动量残差小于10⁻⁵，质量残差小于10⁻⁶，反应进程变量残差小于10⁻⁴。所有计算在Think Station P700工作站完成，耗时约100 h。

3.   结果分析

油气爆炸关键特征通常包括火焰传播行为和超压变化规律，为验证数值模拟结果的准确性，选取火焰传播形态、火焰锋面传播速度和超压曲线三方面进行数值模拟，并和文献[17]实验数据对比。

3.1   火焰传播行为

图4为火焰形态的实验和数值模拟结果。在初始阶段（32 ms前），火焰受到管道侧壁的约束沿轴向拉伸，以典型的“指形”火焰传播，且此时火焰锋面未受到分支结构的扰动作用，火焰锋面包络面清晰光滑，为层流状。当火焰传播至分支管道附近时，因壁面约束消失和分支管道尖角的障碍物作用，致使流场波动剧烈，火焰从层流状向湍流状转捩，火焰锋面褶皱弯曲。火焰传播至分支管道（54～64 ms）时，火焰锋面形成了三路传播路径，即火焰继续沿主管道向末端传播，同时在两侧分支管道中对称传播。在分支管道中，火焰锋面因受到右侧壁面的约束，在传播中逐渐向左侧壁面弯曲，并形成了实验中对称的“浪花状”火焰。并且，随着反应的不断进行，狭长火焰锋面两侧大量未燃气被同时引燃，火焰锋面严重变形并形成明显褶皱。

图  4  实验和数值模拟的火焰形态

Figure  4.  Experimental and simulated flame structures

下载: 全尺寸图片 幻灯片

由图4可以看出，火焰在主管道中传播时，数值模拟能够直观准确地反映实验中的火焰形态变化。当火焰传播进入分支管道时，数值模拟能够较准确地反映实验中“浪花状”火焰的轮廓和锋面位置，但在细节上和实验仍存在一定差别。究其原因有两点：一方面，高速摄影仪在拍摄较大尺寸物体时像素数会呈指数降低，实验照片对火焰的细节捕捉不够；另一方面，数值模拟采用了Zimont模型计算火焰传播锋面，而该模型为了更好地捕捉到火焰锋面位置，对火焰锋面进行了增厚处理，导致在60和62 ms时对分支管道中火焰锋面末端弯曲细节表达不精准。鉴于本文的研究目的，是为了分析火焰与流场耦合规律，因此可以忽略。

在火焰形态对比基础上，对火焰位置和传播速度进行了对比，如图5～6所示。由图5可以发现，数值模拟的火焰位置能够较好地反应火焰传播位置随时间的变化，最大误差为14.436%，发生在42 ms处。这是因为，此时火焰锋面传播至分支管道附近，爆炸流场复杂，火焰锋面结构出现严重溃散，扰动随机性较大，导致误差增大。由图6可以发现，数值模拟的火焰速度与实验的变化规律整体一致，在爆炸初期（28 ms前），因参与反应的油气较少，反应速率较小，火焰传播速度缓慢升高；随着化学反应持续进行（实验28～40 ms，数值模拟28～45 ms），火焰锋面向末端传播过程中沿轴向不断拉伸，火焰表面积增大，参与反应的单位时间油气增多，化学反应速率增大，进而火焰传播速度迅速增大；当火焰锋面传播到分支管道附近时，初期分支管道的泄压作用较大，火焰传播速度增长趋势出现停滞，其后随着分支管道内油气被引燃，火焰锋面表面积呈指数增加，火焰传播速度继续增大。数值模拟误差最大为24%，发生在45 ms处。

图  5  实验和数值模拟的火焰位置

Figure  5.  Experimental and simulated flame front locations

下载: 全尺寸图片 幻灯片

图  6  实验和数值模拟的火焰速度

Figure  6.  Experimental and simulated flame speeds

下载: 全尺寸图片 幻灯片

3.2   动态超压变化规律

在压力传感器P1处实验和数值模拟的动态超压曲线，如图7所示。表1为实验和数值模拟的3个典型超压峰值的对比。可见，两者超压曲线的规律基本一致，分别在65和68 ms达到爆炸超压峰值，峰值误差为8.6%。但在数值模拟中，未考虑薄膜破裂时的破膜压力对管道内压力场的影响，因此在爆炸超压上升阶段，实验存在3个峰值，而数值模拟只存在2个峰值。

表  1  实验和数值模拟的典型超压峰值

Table  1.  Experimental and simulated typical overpressure peaks

方法	p1/kPa	ε1/%	p2/kPa	ε2/%	pmax/kPa	εmax/%
实验	12	0	25	0	58	0
数值模拟	—	—	27	8.0	67	15.5

下载: 导出CSV 

| 显示表格

图  7  P1处实验和数值模拟的爆炸超压曲线

Figure  7.  Experimental and simulated overpressure curvesat point P1

下载: 全尺寸图片 幻灯片

在不考虑破膜压力的情况下，爆炸超压升压过程可以划分为4个阶段，以数值模拟得到的曲线为例：第Ⅰ阶段为缓慢升压阶段（0～20 ms），该阶段从模拟点火开始至火焰拉伸呈“指形”终止，在该阶段爆炸集中于点火源附近，参与反应未燃气量少，反应释能小，超压随之上升缓慢；第Ⅱ阶段为超压上升阶段（20～50 ms），随着爆炸进行，参与反应的未燃气逐渐增多，单位时间释能迅速增大，超压开始明显上升；第Ⅲ阶段为超压平衡阶段（50～56 ms），主管道泄流效应和分支管道二次泄压作用弱化了内场超压上升，超压上升速率减慢；第Ⅳ阶段为快速升压阶段（56～65 ms），分支管道增大了燃烧面积，反应速率迅速上升，超压再一次快速上升，直至可燃气消耗殆尽转入下降阶段。

在小尺度狭长空间发生油气泄爆时，通常管道内压力波动以静压为主，动压对超压的影响非常小，因此空间内部各点的超压基本一致。而对于含双侧分支的受限空间，当压力波传播经过分支附近时，内部空间横截面的变化必然影响压力波的传播，因此在分支管道前后测点的超压一定会有所差异。为此，对比了分支管道前后测点P1和P2处的超压曲线，如图8所示。由图8可以看出，P1和P2处的超压曲线演变规律基本一致，仅在局部有数值大小等方面的差异，具体有3点：（1）在缓慢升压阶段和超压上升阶段，P2处超压上升滞后于P1；（2）P1处最大超压峰值大于P2，形成时刻略滞后于P2；（3）P2处的超压振荡比P1剧烈。P1和P2处的超压特性差异说明，分支管道导致主管道横截面突扩，在分支管道附近诱导稀疏波导致压力波衰减，分支管道对主管道起到了泄压作用，并且P2位置更靠近泄爆口，受泄流效应的影响也更大。

图  8  P1和P2处数值模拟的爆炸超压曲线

Figure  8.  Simulated overpressure curves at points P1 and P2

下载: 全尺寸图片 幻灯片

3.3   流场-火焰-超压耦合

由文献[22]，在爆炸发展过程中，流场变化会影响火焰传播行为，火焰传播行为的改变亦会对爆炸流场产生影响，火焰传播和湍流流场间存在一定的正反馈耦合机制。此外，火焰结构在一定程度上表征反应速率的快慢，进而影响爆炸超压的波动。

图9为火焰传播形态和流场速度矢量图，白色箭头代表速度矢量，颜色越白速度数值越大，速度矢量所在平面为管道中轴面（y=0.05 m）。由于重点关注分支管道附近的火焰和流场形态，因此火焰传播至分支管道前（20、32 ms）的图片包含点火端至分支管道；而为了清晰显示火焰在分支管道中传播形成“浪花状”的过程，这个阶段（48、60、64 ms）的图片仅包含分支管道区域。

图  9  火焰传播形态和流场速度矢量

Figure  9.  Flame propagations and velocity vectors of flow field

下载: 全尺寸图片 幻灯片

在20 ms时，火焰以层流状态传播，火焰锋面前端附近速度矢量沿半球面指向未燃区，主管道未燃区气体在火焰膨胀作用下沿水平方向向管道末端流动。当流动经过分支管道时，两侧壁面对未燃气的约束消失，未燃气对称地向两侧分支管道膨胀，在主管道与分支管道交界处建立了层流流场，此时分支管道内部流场速度趋近于零。随着反应进行，在32 ms时火焰以“指形”传播，火焰面积增大，主管道轴向速度场迅速增强，由于受到分支管道尖角的扰动流场在分支管道口附近向湍流转捩，且因分支管道右侧壁面的约束，在主管道与分支管道交界处形成了对称的涡旋结构，上侧旋转方向为逆时针，下侧为顺时针，此外分支管道内部流场速度逐渐增大。当反应进行到48 ms、火焰传播至分支管道口时，火焰锋面在复杂湍流流场作用下显著变形，而早期建立的对称涡旋结构在火焰冲击作用下向分支管道内部发展且旋转方向维持不变，但涡旋结构沿分支管道轴向拉伸成为较大的椭圆形，此时火焰面积较小且受到较强的泄压作用，导致火焰传播速度上升趋势减缓。此后，火焰传播进入分支管道，在漩涡结构作用下火焰锋面沿分支管道右侧壁面拉伸，而后向左侧壁面弯曲形成60 ms时的“浪花状”火焰。与此同时，火焰锋面前侧未燃气不断被引燃，火焰膨胀促进分支管道内部速度矢量指向主管道，分支管道内未燃气在此流场作用下进入主管道，导致主管道与分支管道交界处形成新的燃烧面，如64 ms时。

为了便于分析火焰形态和动态超压的耦合关系，将P1处数值模拟的动态超压曲线和典型时刻火焰形态结构一起呈现，如图10所示。根据前述分析，爆炸超压变化过程可划分为4个升压阶段和1个振荡下降阶段。第Ⅰ阶段，在点火初期，只有点火头周围少量油气参与化学反应，而后火焰以半球形传播，参与反应的可燃气有限，爆炸缓慢发展，超压上升趋势不明显。第Ⅱ阶段，爆炸进入快速升压阶段，初期火焰以“指形”结构传播（32 ms），火焰面积随着锋面拉伸不断增大，参与反应的油气增多，超压逐渐上升，当火焰传播至分支管道附近时，因流场复杂导致的火焰锋面褶皱严重进一步导致反应更加剧烈，超压上升趋势更明显。第Ⅲ阶段，火焰传播进入分支管道形成“伪受限空间”燃烧（54、56 ms），初始火焰传播速度较慢，火焰面积变化不明显，参与反应的油气量没有显著增加，而火焰锋面靠近末端泄爆口，主管道泄压作用增强，导致这个阶段升压速率下降。第Ⅳ阶段，当火焰在分支管道中涡旋结构的作用下充分发展时，火焰锋面不规则溃散严重，同时参与反应的油气量也大幅增多，爆炸剧烈进行，升压速率较前3阶段都大。振荡下降阶段，火焰传播到64 ms时，由于主管道末端的泄压效应，超压迅速下降，而后在压力波来回传播和管道内残余油气持续爆炸的两方面作用下形成超压振荡。

图  10  数值模拟的超压曲线和火焰传播形态

Figure  10.  Simulated overpressure curve and flame propagation structure

下载: 全尺寸图片 幻灯片

4.   结　论

利用大涡模拟WALE模型，模拟了双侧分支结构受限空间油气泄压爆炸的超压特性和火焰行为，与实验数据对比了火焰传播形态、火焰传播速度和动态超压曲线，验证了所选取模型对含多分支结构受限空间泄压爆炸计算的有效性。然后，运用该模型对双侧分支结构受限空间油气泄压爆炸的空间流场结构、火焰行为和超压特性进行了研究。得出以下结论：

（1）火焰传播进入分支管道前，火焰行为会影响分支管道内部流场，在主管道和分支管道交界处会产生方向相反的对称涡旋结构，随着火焰进一步传播，涡旋结构会向分支结构内部发展，建立旋转流场；

（2）当火焰传播进入分支管道后，分支管道内部前期已建立流场决定了火焰的形态，火焰锋面在涡旋结构作用下形成“浪花状”火焰，此后火焰和流场相互影响，流场向湍流转捩，火焰锋面溃散；

（3）爆炸超压升压过程可划分为4个阶段，受到火焰锋面面积、主管道泄流效应和分支管道二次泄压共同作用，爆炸流场、火焰行为和动态超压具有显著耦合性。