Large Eddy Simulation of Gasoline-Air Mixture Explosion in Long Duct with
Branch Structure

Gas explosion is a process involving complex hydrodynamics and chemical reactions. In order to investigate the interaction between the flame behavior and the dynamic overpressure resulting from the explosion of a premixed gasoline-air mixture in a confined space, a large eddy simulation (LES) strategy coupled with sub-grid combustion model has been implemented. The considered confined space consists of a long duct and four branches symmetrically distributed on both sides of the long duct. Comparisons between the simulated and experimental results have been considered with regard to the flame structure, flame speed and overpressure characteristics. It is shown that the explosion process can qualitatively be reproduced by the numerical simulation. Due to the branch structure, vortices are generated near the joint of the branch and long duct. Vortices rotate in opposite directions in the different branches. When the flame propagates into the branch, the flame front is influenced by the flow field structure and becomes more and more distorted. The overpressure displays a similar behavior in the two branches which have a different distance from the ignition point. It is finally shown that the overpressure change law can directly be put in relation with the shape of flame front.


Introduction
The flammable gas explosion is a common accident type in industrial areas and usually brings large losses in economic and human life [1]. In order to reduce losses, researchers have carried out a lot of indepth research on various combustible gases. It is found that the space structure and shape have great influence on the law of explosion especially in the overpressure time history and flame propagation. The differences between the space structure and the shape contain the L/D ratio [2], vent size or location [3][4][5], obstacle shape or number [6][7][8], and turnings or bends [9,10]. Especially for the turnings or bends, researchers have made some findings and discoveries. Investigations show that 90-degree bend in a duct will result in pressure gradient and vortices, thus enhancing flame speeds and overpressure and to some extent the effect is comparable to baffle type obstacles [9,11]. The effect of directions of the turnings on explosion characteristics are also studied through experiment [12]. The results show that the explosion process is largely depend on the number of turnings and less depend on the directions of the turnings. The pipe network or parallel pipes are also play an important role in flame propagation and overpressure dynamics [10,13]. Experimental studies show that the law of peak overpressure and flame speed are affected by the ignition position and the branch length connected to the parallel pipe. The explosion characteristics in parallel pipe with two unequal branches are more complex than those with two equal branches, and a high-pressure area is created in the middle of the transversal branch.
As a prevalent structure, ducts with branches are very common in oil depot and pipeline. However, explosion in ducts with branches is far different from those with turnings or bends because of the multiple path and end of flow [14]. The duct with branches means two or more paths and end of flow which made the explosion flow structure more complex. Li et al. [15] found that the branch structure has the ability to enhance the peak overpressure, overpressure rise rate and flame speed of gasoline-air mixture explosion and it mainly influences downstream flow of the branch structure. Zhang et al. [14] studied explosion in five different branch configurations through experiment. It is found that deflagration index increased as a function of the branch number. However, there is an optimal branch number for flame speed, which makes the flame speed increasing most. Although researchers have made some research on explosion in duct with branches, unfortunately, these studies mainly focus on phenomenal descriptions, including the law of overpressure dynamics, flame structure and flame speed. Although the flame-vortex is the key process in an explosion description [16], the coupling relationship among flow field, flame propagation and overpressure dynamics is less involved.
As the progress of computer hardware performance and the development of numerical algorithms, more and more researchers introduce numerical simulation method into the study of combustible gas explosion and help to investigate the characteristics of flow field during the explosion [17,18].
In order to study the influence of the branch structure on the explosion characteristics, a numerical simulation model is established in this paper and verified by experiments. Through the simulation, the interaction among flow field, flame speed/shape and overpressure dynamics during gasoline-air mixture explosion in a long duct with four symmetrical settled branches is studied.

Test Case and Experimental Setup
This study was based on experiment about gasoline-air mixture premixed combustion conducted in confined space which contains a long duct and 4 branches as shown in Fig. 1. Both the long duct and 4 branches had the same rectangular cross-section size of 100 mm × 100 mm, made of three sides of steel plate and one side of plexiglass plate. The thickness of both steel plate and plexiglass plate were 200 mm.
The experimental system consisted of a gasoline-air mixture generation and testing system, a pressure acquisition system, a high-speed camera, an electric ignition system and a computer, as shown in Fig. 2. The gasoline-air mixture generation system consisted of a vacuum pump, an oil bottle, four ball valves and a gas transmission pipe. The working principle of the system was to use the airflow generated by the vacuum pump to drive the acceleration of gasoline evaporation in the oil bottle. Driven by the airflow, the gasoline vapor entered the confined space and mixed with air. The gasoline-air mixture concentration was tested by a GXH-1050 infrared gas analysis apparatus. The pressure acquisition system consisted of two pressure transducers, a dynamic testing analyzer and a dynamic signal testing software. The type of the pressure transducer is piezo resistive, and its measurement range was 0-2 Mpa. The dynamic testing analyzer was supported by the signal testing software, and both of them were designed and produced by DongHua Testing Technology company. The high-speed camera (GC-P100BAC, JVC) was used to capture the detailed flame front at the branches. The frame frequency was set at 500 Hz. The electric ignition system consisted of an ignition plug, a capacitor and a transformer. The ignition plug was set at the center of the right-side blind flange and its ignition energy was set at 6 J.
In order to obtain the most remarkable effect of the explosion overpressure and flame speed, the initial gasoline-air mixture of this paper was 1.7%, which was close to the stoichiometric concentration. Due to the variability of gasoline distillation product, its experimental repeatability is not as good as that of the pure substance. In order to reduce the error and ensure the accuracy of the results, the experiment was carried out for three times.

Large Eddy Simulation (LES) and Combustion Model
This paper used LES model to simulate gasoline-air mixture explosion process in a confined space with 4 branches. The LES model has been used to simulate transient flammable gas explosion in some papers previously [18][19][20]. The idea of LES model is based on the separation between large and small scales through the Favre-filter. LES governing equations is obtained by filtering on the time-dependent Navier-Stokes equations. Then, the equation set resolves the large turbulent structure directly, while small scales would be resolved by subgrid-scale models. In LES model, species transport equation is resolved by using a reaction progress variable c which is also used to model the flame propagation [21]. The equation can be expressed as: In the equation, the over bar (-) denotes a filtered quantity and the tilde ($) denotes a Favre-filtered quantity. ρ is the mass density, u is the velocity, l eff is the effective viscosity and derived from the renormalization group (RNG) theory [22], Sc eff is the effective Schmidt number, S c is the reaction progress source term. For the effective Schmidt number, assumed to equal to the effective Prandtl number [22].
For the effective viscosity l eff , its detailed expressions can be written as: where H(x) is the Heaviside function [23], l is the dynamic viscosity, eddy viscosity , V cv is the volume of a computational cell,s ij is the rate of strain tensor.
For the reaction progress source term S c , its detailed expressions can be written as: where q u is the density of unburnt mixture, U t is the turbulent flame speed.
In this paper, in consideration of the wrinkling and stretching of the flame front, the turbulent flame speed was computed by TFC model or wrinkled and thickened flame front proposed by Zimont et al. [24]. The turbulent flame speed U t can be written as: where A is a model constant and the value is taken to 0.5 [25], u 0 is the root-mean-square velocity written as u 0 ¼ l t s À1 sgs , U l is the laminar flame speed, a is unburnt thermal diffusivity, l t is the turbulence length scale.

Numerical Details
The structure and dimension of the computational domain were the same as the confined space described in Fig. 1. It was a 4000 mm-long duct contained four 500 mm-long branches. The whole computational domain consisted of about 1,685,000 hexahedral cells and was characterized with the size of 2 mm to 10 mm through the BiGeometric meshing law. The grid size is mainly 10 mm in the domain, meanwhile the size is set to 2 mm at the decussations. Both ratio 1 and ratio 2 from the elements of 2 mm on the curve are set to 1.05. Due to the flame spreads faster, the heat transfer between the flame and the wall is less. Then all faces of the confined space were set to walls and adiabatic and no slip boundary conditions. The physical properties of the combustion medium were set to be the same with the gasoline-air mixture. The laminar flame speed was considered as a constant with pressure and temperature, and equal to 0.43 m/s [26]. The molecular viscosity was calculated through the Sutherland's law. The specific heat was considered approximately following the piecewise fifth-powder polynomial function of temperature [27]. The initial concentration of the fuel is set to 1.7%, which is about the equivalent concentration of the gasoline-air mixture.
The initial temperature was set to 300 K, while the initial pressure and reaction progress were set to 0. A hemisphere with a radius of 5 millimeter was marked at the left end of the duct and the reaction progress was patched as 1 in this hemispherical to simulate the ignition progress [21]. Equations of the model were solved by the ANSYS Fluent software, and the whole computational domain was discretized by the finite volume method. The standard SIMPLE algorithm was chosen as the pressure and velocity coupling method. The second-order upwind scheme was used to calculate the face values of convection terms and the second-order-accurate central-differencing discretization scheme was used for diffusion terms. The time step size was set to 1 × 10 −4 seconds, and each time step contained 20 iterations. During the iteration, the residual for energy equations, momentum equations, mass equations, progress variable equations were less than 1 × 10 −6 , 1 × 10 −5 , 1 × 10 −4 , 1 × 10 −4 . The calculation was computed on a Thinkstation workstation with 12 processors and 16G RAM. The time for a calculation was about 100 hours.

Model Validation
The most important factors to evaluate gas premixed combustion are the flame structure, flame speed and overpressure time histories. In this work, these three parameters were also used to assess the accuracy of the simulation model through comparing the results between the experiment and simulation. Because of the main goal is to analyze the influence of branch structures on flame propagation, so the comparisons of the flame structure and flame speed were emphasized at the area near branches.
Comparison between experimental and simulated flame structures at four different instants are shown in Fig. 3. The parameter used to simulated flame structures in simulated results is progress variable. Although some of the details of flame structure in the experiment were not fully captured due to the limitation of the camera, we can also find that the simulation reproduced the flame structure effectively. Fig. 4 shows the comparisons of flame speed versus time in the pipe and branch structure between experiment and simulation. The comparison shows that the numerical simulation can reproduce the trend of flame propagation speed in the experiment, but the value obtained by numerical simulation is slightly larger than the experiment. This is because the wall boundary condition in the numerical simulation is set as the adiabatic wall, but in the actual experiment, there is heat dissipation.
The flame speed S f is calculated by the equation: where x nþ1 À x n means the distance between two flame fronts with time intervals of Dt n , x n was defined as the distance from the leading point of the flame front to the ignition end at the distant t n .  Fig. 5 show that the numerical simulation reproduced the overpressure change law especially the gentle increase in the initial stage of flame propagation and the oscillation as the flame propagates through branch structures. The difference is that the overpressure obtained by numerical simulation is larger than the experimental value at the end of the explosion. This is due to the larger flame speed and combustion rate obtained by numerical simulation, which leads to greater overpressure. Fig. 6 shows the 2D slices of flame front and flow field variation. The slice lies in the middle of geometry model. The white lines represent the streamlines in the confined space and the arrows indicate the direction of flow field. In order to clearly show the flame structure and the flow field at the branch, the flow field structure at the branches at each moment is separately amplified. Since the geometry model is axially symmetrical, only the branches on one side of the geometry are studied.  Interaction between the flame front and flow field structure around the branch can be divided into two stages. The first is the time before flame propagates to the first branch. At this stage, the flame propagation helps to establish the flow field in the branches. The second stage is the flame propagates into the branches. At this stage, the flame front is influenced by the flow field structure.

Interaction between Flame Front and Flow Field Structure
At 20 ms, the flame front propagates in sphere, streamlines are concentrated in the middle of the long duct and are horizontal at the region far from the flame front. Near the two branches, due to the sudden expansion of the cross section, streamline bends to the side of the branch structure respectively. However, as can be seen from the figure, the distribution of the streamline in the two branches is different. The reason for this difference is that the velocity field is mainly concentrated in the middle of the duct before the first branch. After the first branch, the flow velocity in the duct is reduced but the distribution of streamlines is more uniform, and the streamline is more uniform when propagating into the second branch. As the flame propagating, the flame front is finger-shaped and the streamlines distributed horizontally throughout the cross section. Small vortices begin to appear at the branch corners near the ignition point side, as shown in Fig. 6 at 52 ms. When flame front gets flat at about 92 ms, at the four branch inlets, two pairs of corresponding vortices are generated. The upper vortex rotates counterclockwise and the lower vortex clockwise. The vortex velocity at the first branch is greater than the vortex velocity at the second branch. When the flame propagates in tulip shape, the flame propagation speed is small. Since the flow field in the duct is weak at this time, the distribution density of the streamlines in the branches is relatively weak and the vortex in the second branch moves toward the center of the duct under inertia. When the flame propagating into the branch as shown at 214 ms and 238 ms, the flame is stretched along the branch wall by the influence of the flow field and the flame front becomes more distorted. The distorted flame fronts further lead to complex flow fields.

Interaction between Flame Front and Overpressure Dynamics
Comparison of overpressure between branch 1 and branch 2 is shown in Fig. 7. Through the comparison, it can be seen that in the development of explosion, the changes of the overpressure in the two branches are similar. The only difference is the overpressure oscillations at the end of the rapid rise of explosion overpressure. The overpressure oscillations in the first branch are weaker than those in the second branch. According to the analysis of the flow field, the reason is that the first branch enhances the turbulence of the flow field near the second branch. This enhancement makes the flame propagates into the second branch more distorted than the first branch and then influences the overpressure oscillation.
Taking the overpressure time histories in the second branch as an example, Fig. 8 shows the corresponding relationship between the flame front and the overpressure curve. Before 16 ms, the flame propagates in hemisphere and the overpressure rises steadily. Then, due to the restricted role of the wall, the flame stretches along the axial direction and develops into the finger shape. The increase of flame surface area leads to a high combustion rate and the overpressure starts to rise. As the flame skirt touches the wall and the flame front gets flat at 92 ms, the overpressure reaches the first peak of 0.0218 MPa and then starts to reduce due to the decrease of flame surface area. After that, in the process of evolving into tulip flame, the overpressure continues to rise until 164 ms when the overpressure reaches the second peak of 0.0509 MPa. As the flame front arrives at the first branch at 212 ms, the flame becomes more and more distorted and enhancing the combustion rate and heat release rate, and therefore the overpressure rises rapidly. While the flame arrives into the second branch, the combustion in the pseudoconfined space constructed by flame front and the wall of the second branch resulted into the overpressure oscillation like 280 ms and 282 ms [28].

Conclusion
In this paper, an experimental study and LES numerical simulation are performed to investigate characteristics of the gasoline-air mixture explosion in a long channel containing four branch structures. Comparisons between experimental and simulated results on flame structure, flame speed and overpressure time history are made in this paper. Through these comparisons, it is found that LES simulation has the ability to reproduce the process of the explosion evolution. Besides, the interaction between the flame front and flow field structure and the overpressure dynamics are studied through the simulated results.
In the early stage of explosion, the flame front affects the flow field distribution around the branches. According to the whole process of the explosion, the flow fields are different in branches that have different distances to the ignition point. Due to the branch structure, vortices are generated near the joint of the branch and long duct. The vortices in the opposite two branches rotate in opposite directions. Of the two branches at different distances from the ignition point, the vortices in the branches closer to the ignition point are stronger, but the vortices in the two branches rotate in the same direction. When the flame propagates into the branch, the flame front will be influenced by the flow field structure and become more and more distorted. Overpressure dynamics are similar in the two branches which have different distance to the ignition point, while the only difference is the strength of overpressure oscillation when the flame propagating into both of the two branches at one side. This oscillation is caused by the pseudo-confined space combustion phenomenon. The overpressure change law is positively correlated with the shape of flame front. Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.