Parametrical investigation for the optimization of spherical jet-stirred reactors design using large eddy simulations

Due to the importance of gas-phase chemical reaction kinetics in low-emission combustion, stirred tank reactors have been used for decades as an experimental tool to study high- and low-temperature oxidation. A Jet-Stirred Reactor (JSR) setup is valuable to determine the evolution of species mole fractions. For the accuracy of the experimental results, it is important that a JSR is designed such that the concentration field is as homogeneous as possible in order to avoid disturbance of the chemical kinetics. In this work, numerical simulations were performed to investigate the mixing in a JSR chamber. The turbulent structures inside the JSR and the nozzles are captured using Large Eddy Simulations. We conducted numerically a parametric study to evaluate the effects of thermodynamic conditions and geometrical parameters on the mixing characteristics. More specifically, the diameter of the spherical chamber is modified together with the diameter of the nozzles through which fresh gases are fed. The characterization of the gas flow inside a typical spherical JSR layout and results derived by the normalized standard deviation of a tracer mass fraction show that a reduction of the JSR diameter at high pressures improves the homogeneity. Further, we propose a new optimized configuration consisting of six nozzles pointing to the center of the reactor which provides a more uniform composition compared to the standard JSR design.


Introduction
In the last decades, a high demand for clean combustion and a high energy conversion efficiency emerged to meet the stringent environmental standards on pollutant emissions. To this end, the chemical kinetics of alternative/renewable fuels have been extensively investigated through lab-scaled model experiments, like flames or chemical reactors [1,2]. The development of the kinetic mechanisms is supported by measurements in continuous and batch reactors [3]. Accordingly, validation and developing extensive chemical reaction mechanisms in such systems reduce total costs and efforts in engine manufactures. A Jet-Stirred Reactor (JSR) as an ideal continuous stirred tank reactor is particularly suitable for studies of gas-phase kinetics. To analyze the chemical reactions in detail, JSRs are usually coupled with diagnostic tools. The homogeneity of the spherical JSRs under different experimental conditions is the main requirement to investigate the simple 0-D kinetic modeling of this reactor type. Therefore, enhancing the homogeneity of the composition is achieved by chaotic advection through stirring, provided by turbulent jets inside the spherical JSR. In fact, the mixing rate and, thus, the realized level of homogeneity is the key factor to obtain chemical kinetic mechanisms of an optimized accuracy. Furthermore, simulations like Computational Fluid Dynamics (CFD) coupled with a developed kinetic mechanism are excellent for providing a general model that can be applied in several applications; consequently, this reduces the expense of experiments.
In a JSR experiment, premixed reactants are introduced into the reactor by a number of nozzles and well-mixed products are continuously removed at the outlet. These nozzles can be of different designs, e.g., inclined, crossed, or disk shape. Also, JSRs have so far been operated in various pressure and temperature ranges [4]. Davani and Ronney [5] even proposed a configuration of eight jets, each surrounded by a concentric outlet aiming to improve the mixing and consequently experimental accuracy. Later, the jet arrangement was improved by Zhang et al. [6] by adding four holes on the reactor wall to remove the exhaust gas. They concluded that the design of an inwardly offcenter shearing chamber generates small vortices over the chamber. These vortices promote the reactor performance in terms of mixture uniformity and the distribution of species compared to a traditional JSR. Esmaeelzade et al. [7] studied the mixing and flow field characteristic inside a JSR with four inclined nozzle through highly-resolved CFD analysis. Their results agreed well with the experimental data from the literature [8], especially for lower flow rates. Further, De Oliveria et al. [9] simulated an enhanced JSR configuration by considering a hemispherical premixed chamber and an outlet pipe close to the equatorial plane of the reactor. They reported that the improved geometry results in an increase of turbulence, mixing efficiency, and homogeneity of the temperature profile. All these discussed studies demonstrate the potential to enhance the mixing quality dependent on the geometrical configuration of the reactor.
Typically, the experimental setup consists of an electrically heated oven and a pressure-resistant jacket which completely enclose the JSR. Due to these enclosures only a limited optical accessibility to the flow field is given and detailed experimental information concerning the complex interaction of the jets is not available.
MacMullin and Weber [10] first proposed to characterize the macro-mixing inside real chemical reactors compared to ideal ones based on their Residence Time Distributions (RTD). The RTD is determined in general by injecting a tracer into the reactor and then monitoring the purged tracer concentration in the effluent stream as a function of time [3]. Wheres the RTD is an integral quantity, a detailed picture of the condition in the reactor is still lacking. Thus, the need to gain deeper insight in the mixing in a JSR motivates CFD simulations. A more detailed understanding of the flow field characteristic would enable to facilitate further improvement of the JSR design [11,12] which would finally lead to the more accurate characterization of chemical mechanisms.
In the context of our study, we analyze the flow generated inside the JSR by means of unsteady three-dimensional LES. Based on these data we attempt to identify flow patterns such as dead zones and short circulations which are responsible for the non-homogeneity inside JSRs. The results of our simulations are compared to the RTD data of the previous experiments by Ayass et al. [8] and CFD analysis by the authors of the present paper [7]. Our results indicate that some of the commonly used reactors could have significant concentration non-homogeneity and therefore may not be suitable for chemical kinetic studies. Further, we perform a parametric study to optimize the JSR in terms of the diameters of the nozzles and the chamber, and the temperature and pressure. Finally, we suggest a new JSR design including six impinging jets which provides a more efficient mixing compared to the typical JSR configuration of four jets.

Mathematical model and numerical setup
In this section, the details of the underlying mathematical model, numerical methods, the studied JSR geometry, and the employed grid resolution are described.

Applied mathematical model
The numerical simulations were performed using the open source software package OpenFOAM [13]. To solve numerically the transient, turbulent flow of a compressible non-reacting mixture, we have used the reactingFOAM solver. Therein, the Favre-filtered form of mass, momentum, energy conservation for Newtonian fluids, species transport equation for each involved species together with the equation of state for ideal gases are solved in Eulerian framework given by [14] where the symbol ̄ is used for Favre-filtered variables in LES; the ∼ is applied to the filtered density-weighted variables. In the above equations, velocity u and mass fraction Y k are functions of both space and time, t. In addition, density, effective kinematic viscosity, pressure, effective mass diffusivity of the species k and identity matrix are denoted by , eff , p, D k and I respectively. On the right hand side of Eq. (3) the source term ̇k represents the source terms of species recombination during chemical reactions. Since we consider a non-reacting flow, in the present study these terms are neglected. Furthermore, the pressure-velocity coupling is achieved using the PIMPLE algorithm which is a hybrid between PISO [15] and SIM-PLE [16]. We chose PIMPLE loop in order to increase the Courant number limitation and to ensure the simulation stability.
In the present study, the solver uses a thermophysical model based on the compressibility of the fluid for a system with a fixed composition. Moreover, the specific heat capacity of pure gas is obtained by the JANAF thermochemical tables [17]. Therefore, the thermophysical properties of the mixture are calculated by the mixing law [18].
Moreover, the governing equations are discretized using the FVM method. Since the generation of a fine grid compensates the necessity of the higher order discretization, a first order approximation is applied for time advancement and second order for all spatial derivatives. The linear equation system for the pressure is solved using the PCG (Preconditioned Conjugate Gradient) method and the preconditioned DIC (Diagonal-based Incomplete Cholesky) in which the reciprocal of the preconditioned diagonal and the upper coefficients divided by the diagonal are calculated and stored [19]. The absolute tolerance for pressure was set to 10 −6 and the relative tolerance to 0. Furthermore, the PBiCG (Preconditioned Bi-conjugate Gradient) solver together with the pre-conditioner DILU (Diagonal-based Incomplete LU) are chosen for the velocity and turbulent kinetic energy. The absolute tolerances of the mentioned variables are equal to 10 −6 and the relative tolerances to 0. Turbulent flows are characterized by a wide range of length and time scales. The kinetic energy in the flow transfers from large to the small scales, eventually being dissipated at the smallest scale. The large scales are most energetic and transport of the momentum and other conserved fluid properties are mostly carried out by these scales. Because of that, in the LES philosophy large scale motions are solved directly on the grid and small scales are modeled. In the filtering procedure of the LES, each variable is decomposed into a filtered component which will be resolved by Eqs. 1, 2, 3 and a sub-grid component. Therefore, the effective kinematic viscosity is the sum of the molecular viscosity and the residual or Sub-Grid Scale (SGS) eddy viscosity sgs . In this study, to reduce the computational costs and achieve a reasonable accuracy, we have chosen the standard Smagorinsky model. Thus, residual (sub-filter) scale stresses are accounted for via the Smagorinsky model [20], where the SGS viscosity is defined as where |S| , Δ , and C s represent the filtered magnitude of the rate of strain tensor, the characteristic size of the filter, and the static Smagorinsky model coefficient, respectively. The classical Smagorinsky constant is calculated as C s = C k √ C k ∕C where the default constants of the Smagorinsky model implementation in OpenFOAM C k and C , are 0.094 and 1.048 respectively [21]. Having the dynamic viscosity of the pure gases available based on Sutherland's law, the mixture viscosity of CO 2 -N 2 composition is calculated according to the mixing laws [18].
Likewise, the effective mass diffusivity consists of the resolved mass diffusivity ∕Sc and the sub-grid scale turbulent mass diffusivity t ∕Sc t where the turbulent Schmidt number Sc t is assumed to be equal to unity.
The fluid elements entering the JSR continuously withdraw from the reactor and the time which various elements stay inside can be different. Some molecules may remain in the reactor indefinite while others almost immediately leave the reactor. The average time which the elements have spent within the boundaries of the reactor is defined as mean residence time, = V ∕Q , where V is the reactor volume, and Q the volumetric flow rate [1]. Numerical calculation of residence time is based on the tracer concentration function C(t) as follows [3,22] Equation 5 relies on the exponential decay of the tracer concentration leaving the outlet C out at the time, t, which is related to the initial steady state concentration of the tracer C 0 . It is calculated based on the tracer concentration at t divided by the area under the concentration curve.

Geometry and conditions of the parametric study
The base configuration of the JSR in this study follows the design of a spherical reactor with four inclined nozzles proposed by Dagaut et al. [11]. The exemplary computational domain illustrated in Fig. 1 is composed of a spherical chamber of a diameter of 40mm with four nozzles each of a diameter of 1mm in the equatorial plane. The upward and downward pointing nozzles together with the nozzles pointing horizontal generate two circulating streams, parallel to the x-z and x-y planes, respectively (see the arrows at the nozzles' outlet in Fig. 1). The outlet of a diameter of 10mm is placed in the surface parallel to the x-y plane. The internal field of this geometry was meshed. The generated (4) sgs = (C s Δ) 2 |S| , Vol:.(1234567890) Research Article SN Applied Sciences (2021) 3:766 | https://doi.org/10.1007/s42452-021-04743-w mesh was unstructured and refined close to the walls, the inlets, and the outlet where higher flow gradients were expected.
The simulated conditions were oriented on the experimental setup of Ayass et al. [8]. A constant temperature of 298K is applied for the inlet boundary and the initial condition, while a zero-gradient condition for the temperature and the species mass fraction is considered at the walls. For the pressure, an atmospheric condition was considered at the outlet and zero-gradient at the walls. Regarding the velocity field, we applied zero initial conditions together with no-slip boundary conditions at the walls, a constant value with a top-hat profile at the entrance of nozzles' orifices and zero-gradient at the outlet. The Reynolds number (for CO 2 having kinematic viscosity of 10.3 E − 6m 2 ∕s ) based on the averaged inlet velocity of 26.9m/s and the orifice diameter of 1mm is calculated as 2620 and the reactor diameter of 40mm as 10400. Furthermore, the time-averaged Mach number reaches 0.18 at the nozzles' orifice where the gas velocity is the highest.
In accordance with the experiment, a constant flow of CO 2 is first injected inside the spherical chamber by aid of four inclined and crossed nozzles situated in the equatorial plane. Then, a solenoid valve sharply cuts off the CO 2 stream and injects suddenly a continuous flow of N 2 inside the chamber with different constant volumetric flows ranging from 18 cm 3 /s to 180 cm 3 /s. In accordance with the experiment by Ayass et al. [8], a constant flow of CO2 is first injected inside the spherical chamber by aid of four inclined and crossed nozzles located in the equatorial plane. The outlet concentration of the tracer is then calculated as a function of time. In our previous research [7], the results of of the numerical calculations of the RDT were validated versus experimental data from literature [8] for different flow rates and ambient conditions. In this numerical study, we have kept the value of mean residence time constant ( = 0.37s ) in order to understand the effective parameters which improve the mixing quality inside the JSR. That means the volumetric flow rate and the reactor volume were both considered constant in terms of varied thermodynamic conditions (i.e. Sect. 3.2). However, these variables were changed in terms of varied geometrical parameters (i.e. Sect. 3.2).
Aiming to find an optimized JSR design, we performed a parametric study considering geometrical and thermophysical parameters. More specifically, we investigated the effects of varying the reactor diameter and the nozzle diameter together with an increased operation pressure and various temperatures on the mixing. Table 1 presents a list of conditions considered in the parametric study. Figure 2 shows an overview of the unstructured grid over the 3D computational domain. A combination of hexahedral, polyhedral, and prism-shaped cells were generated non-uniformly by snappyHexMesh [23]. Since the impinging jets provide significant flow gradients the grid is refined in the vicinity of nozzles orifices and the chamber walls and also in the region of the outlet. In total, the mesh contains approximately one million cells each of a volume between 3.8 × 10 −15 m 3 and 3.5 × 10 −9 m 3 , see Fig. 2c.  Table 1 Operating conditions in terms of the diameter of the spherical chamber (D), the diameter of the nozzles (d), temperature (T), pressure (P), and volumetric flow rate (Q) In order to determine whether the wall-adjacent cells are close enough to the reactors solid boundaries we refer to the non-dimensional quantity y + = yu * ∕ . In this equation, y and are the distance of the center of the closest mesh cell from the wall and the fluid kinematic viscosity, receptively. Furthermore, u * represents the friction velocity at the wall which is defined as a function of wall shear stress, p , and density as u * = ( p ∕ ) 1∕2 . In this study, the closest grid point to the nozzles' walls are located at y + = 0.005 − 0.032 and to the wall of the spherical part of the chamber at y + = 0.002 . Due to the grid refinement close to the walls generated by "SnappyHexMesh" tools, the distance y + is smaller than 1. Therefore, the grids and filter are enough fine to resolve the large-scale motions in the near-wall regions [24].

Grid resolution
Additionally, the characteristic size of the Smagorinsky model filter which is defined as a cube root of the minimum cell volume, 0.015mm , shows that the viscous sublayers of the wall-bounded flows were well resolved on the numerical grid [25].
Further, a commonly used quantity in the characterization of turbulence are the Taylor microscale eddies which are intermediate in size between the smallest eddies of Kolmogorov scale and largest eddies of the integral scale [25]. The Taylor microscale, l , is computed as where the turbulence Reynolds number, Re T , is calculated as u 0 l 0 −1 . Here, l 0 and correspond to the reactor diameter (because of the generation of the largest eddies) and kinematic gas viscosity, respectively. The turbulence intensity u 0 is approximated by 5% of the mean velocity at the nozzle's orifice [26]. Therefore, the length scale of the Taylor eddies can be estimated to 2.1mm . This value is 3.36 times larger than the average cell size used in the computations which indicates an appropriate grid resolution.
In order to control the numerical uncertainty and to estimate the accuracy of the discretization, a grid-refinement study on progressively reduced mesh size is performed. According to Roache [27], we first select three significantly different sets of grids and then run simulations to determine the values of key variables important to the objective of the simulation study. It is desired that the grid resolution increment leads to the convergence of the solution.
As reported in Table 2, the solution variable g i which represents the magnitude of the average velocity and mean residence time has been calculated using LES for three grids consisting of 1M, 490k, and 240k cells, respectively. The mean residence time is obtained for the simulated reactor defined through Eq. 6. The apparent order PO can be directly extracted from the three solutions g 1 , g 2 , and g 3 on the three meshes of an average cell size of h 1 , h 2 , and h 3 with a constant grid refinement ratio of r 32 = h 3 ∕h 2 and r 21 = h 2 ∕h 1 using Eqs. (6,7,8).   (6) is solved using fixed-point iteration, with the initial guess equal to the first term [28]. Further, an estimation of the error e a between the fine grid solution and the unknown exact solution can be approximated by From this, the numerical uncertainty on the finest grid can be evaluated by the Grid Convergence Method (GCI). The fine-grid convergence index, GCI fine [28] is estimated as follows, According to Table 2, the apparent order of discretization for the average velocity is 1.07 and for the mean residence time it is 1.71 which are in the expected range [28]. Since the grid acts as an implicit filter in this study, the GCI represents both the numerical and modeling errors. Even though we cannot separate the contribution of the errors, their sum has a minor effect on our key results [29].
Another method for discretization error estimation is to illustrate the numerical error e (subtracting the solution on the fine and the course mesh) regarding calculation of averaged CO 2 mass fraction versus total number of cells N (240k, 490.4k, and 1.2M) which is shown in Fig. 3. Furthermore, the slope of the line depicted in this figure shows the order of discretization accuracy which is earlier specified as second order for the spatial derivatives.

Numerical results
In this section, we first characterize the unsteady flow field and the species distribution inside the spherical chamber of the JSR. Since we have assumed a constant temperature on the boundaries, a significant temperature inhomogeneity is not expected. Then we present the results of the effects of the reactor geometry together with the elevated pressure and temperature on mixing enhancement.

Flow characteristic inside the JSR
In our previous LES [7], a top-hat profile , i.e., a constant uniform velocity of 26m/s for the incoming gas, was assumed at the outlets of the jets. In the present simulations we also model the in-nozzle flow in order to capture the flow field in the reactor more realistically. Thus, a tophat profile is assumed at the beginning of each nozzle. Then, in downstream direction turbulent velocity profile is developed inside the nozzle. The top-hat and the developed velocity profiles are illustrated in Fig. 4a and b, respectively. The time-averaged streamlines extracted from the instantaneous velocity field are visualized in Fig. 5 indicating the direction of the flow. The averages in this study are derived from instantaneous data that correspond to time instances where the flow is already statistically stationary but the injection of N 2 has not started yet. In this case (which corresponds to the operating condition given in the first row of Table 1), the averaging time scale of 30ms is much larger than the time scale based on physical boundaries of the flow of 0.8ms (e.g. diameter of the reactor divided by the highest velocity inside the chamber).
According to Fig. 5a and b, the flow impinges to the in-front wall and bends in the counter clockwise direction due to the location of the entry of the nozzles pointing upwards or downwards. Similarly, a circumferential stream is generated by the two other nozzles in x-z plane as it can be inferred from Fig. 5c. Therefore, two circulating flow pattern formed by the orientation of the nozzles aims to efficiently support the mixing within the reactor. Further, the distribution of streamlines using the surface-LIC (Line Integral Convolution) technique is illustrated in appendix A [30,31].
In order to identify the vortical structures of the flow field, the Q-criterion is used which is derived based on the second invariant of velocity gradient tensor (Q = |Ω 2 |-|S| 2 ). Here, S and Ω refer to strain rate tensor and   [32]. The positive value of Q represents the existence of vortex inside the chamber since the rotation rate is greater than the strain rate. The distributions of instantaneous Q-criterion is presented in the cross section of the x-z plane in Fig. 6. As expected, a high vorticity is generated at the walls of the nozzles due to the high velocity gradient caused by the assumption of no-slip boundary conditions. Subsequently, the inertial flow at the exit of the nozzles create a vortical structure. As soon as the flow stream starts to bend, it grows in a roller structure and then breaks down into small-scale motions along the wall until it almost decays. In fact, the change of the direction results in the formation of a large circulation bubble downstream of the location of the impingement of the jet to the inner surface of the vessel as discussed previously by Hodzic et al. [33]. It is to be noted that the  LES model helps to predict the eddy structures with different length scales which is not captured by previous RANS simulations [34]. According to Fig. 6, a lower vorticity region is formed in the center of the chamber because of the interaction of the two circumferential flows growing towards the spherical walls and also being away from the walls. As a consequence, a dead zone is created in this location which inhibits a perfectly mixed composition.
The total specific turbulent kinetic energy k tot consists of the specific turbulent kinetic energy in both the resolved scales k R , and the sub-grid scales k sgs , which is given by [35] The specific turbulent kinetic energy of SGS scales can be calculated from the contribution of the SGS turbulence model (Smagorinsky) as follows, where C k is the static model coefficient. It is suggested that the LES using eddy viscosity SGS models is a good measure for estimating the grid resolution [36]. Further, the resolved specific turbulent kinetic energy can be evaluated as (11) k tot = k R + k sgs .
In the above equation u ′ j denotes the Root Mean Square (RMS) value of the filtered velocity components in the x, y, and z directions [25].
According to Pope [25], the resolved turbulent kinetic energy contains 80% of the total turbulent kinetic energy in a well-resolved LES simulation [37]. Therefore, we investigate the distribution of the resolved specific turbulent kinetic energy in Fig. 7a. As expected the highest turbulence fluctuations which enhance the mixing inside the JSR are present where the jets impinge on the inner surface of the vessel as discussed previously by Hodzic et al. [33]. This behavior indicates the necessity of impinging jets for the mixing inside a JSR. We specified the value of resolved specific turbulent kinetic energy over three points in Fig. 7a for further comparison of cases mentioned in Table 1.
In order to determine how the specific turbulent kinetic energy is distributed over eddies of many sizes (wave numbers) or over many frequencies eddies, we plot the energy spectrum for one point located at the center of the chamber with the coordinate of [x:0.0054mm, y:0mm, z:0mm ] in Fig. 7b. According to Sisemore et al. [38], the dissipated energy spectra can be derived as a function of a smoothed Fourier amplitude spectrum. Therefore, in order to plot the energy spectrum, the Fast Fourier Transform (FFT) of the velocity field is calculated over a time domain [0.01s-0.06s] which is larger than the integral time scale [39]. As shown in Fig.7b, the specific kinetic energy density shows a maximum value for the frequency numbers f corresponding to the large energy-containing eddies. Further, the spectra of the velocity follows the −5∕3 law for a wide range of frequencies which indicates a state of fully-developed turbulence. It should be noted that the energy spectrum is shown for frequencies smaller than 10 5 s −1 . That means for those frequencies which are properly resolved by the numerical time step. The Strouhal number ( St = f ⋅ D∕U 0 ) varies in the range of [4.9-790] where D is the reactor diameter and U 0 is the averaged value of the velocity at the nozzles.

Parametric investigation of the thermodynamic conditions
In order to study the influence of thermodynamic conditions on mixing inside the JSR, we elevated the initial and boundary conditions of the temperature to 400K and 700K at atmospheric pressure, and also the pressure to 10bar at atmospheric temperature, respectively. Since we kept the theoretical mean residence time constant, the inflow velocity is also kept constant at the entrance of nozzles' orifices. To quantify the mixture homogeneity in the JSR, where Y i is the CO 2 mass fraction in cell i and Y mean is its mean value. The total number of cells is N and V i is the volume of each cell. The normalized standard deviation represents the deviation of the tracer mass fraction from its mean value divided by the mean value in the complete reactor chamber Fig. 8 shows the normalized standard deviation of CO 2 mass fraction for the cases of different temperatures. In general, the normalized standard deviation of all cases rises in the beginning until reaching a maximum value and then decreases thereafter. This trends relates to the experimental procedure where N 2 is suddenly injected in a JSR chamber being filled with CO 2 . Also, in the beginning ( t < 0.05 s) a high deviation between the CO 2 concentration of the different cases is observed. This deviation becomes less ( t > 0.05 s) as fluid leaves the outlet. According to Fig. 8, higher temperature effectuates a higher normalized standard deviation compared to the lower temperature cases. This trend can be explained by the calculated Reynolds number. In this study, molecular viscosity is dependent on temperature through Sutherland's law. Therefore, increasing the temperature causes an increase of the dynamic gas viscosity [40]; because an increase in temperature elevates rapid motion of gas molecules and the collision between them and consequently enhances the gas viscosity. Additionally, increasing temperature leads to decrease density of the gas mixture (CO 2 -N 2 ). As a result, the decreased Reynolds number of the gas mixture inside the chamber tends to damp turbulence and consequently decrease the efficient mixing and the homogeneity of the composition inside the JSR. This fact can be better understood by comparing the distribution of the resolved turbulent kinetic energy in Fig. 7a and 9 (see the values of the specific turbulent kinetic energy of resolved motions over three identical points for both cases). The total and local specific turbulent kinetic energy at a temperature of 298K have greater values compared to the 700K case. Research Article SN Applied Sciences (2021) 3:766 | https://doi.org/10.1007/s42452-021-04743-w Figure 10 shows that a pressure of 10bar results in a lower normalized standard deviation of the CO 2 mass fraction compared to the case considering ambient pressure. This means an efficient mixing and a better homogeneity are achieved at higher pressures inside the reactor. According to the ideal gas law, higher pressure increases the gas density inside the chamber. This prevents the short circulation of the fluid elements. That means the fluid elements are convected on a short path from the inlets to the outlet without mixing with the rest of the gas inside the JSR's chamber. In fact, higher pressure creates more vortexes over the reactor which supports mixing. The time-averaged velocity contours captured with the LIC method compares the flow field of different pressures in appendix A.

Parametric investigation of the geometrical parameters
Further, we studied the effects of geometrical parameters on the mixing inside the JSR considering a constant residence time. We changed the diameter of the nozzles to 1.5mm and also the diameter of the chamber from 35mm to 60mm at ambient thermodynamic conditions. Because of the practical limitations, a JSR with a nozzle and chamber diameter of less than 1mm and 35mm respectively, cannot be built up. In Fig. 11 the normalized standard deviation of the CO 2 mass fraction for different nozzle and reactor diameters is plotted versus time. As can be seen in Fig. 11a, a more homogeneous composition after t = 0.15s is achieved by decreasing the chamber diameter. However, the reactor having D= 60mm reaches the maximum of standard deviation later (at t = 0.15s ) compared to the diameter of D = 40mm (at t = 0.05s ). Because N 2 is injected to a JSR  which is already filled with a homogeneous composition of CO 2 , the fluid needs to travel a longer distance and time to the outlet in the case of larger chamber with the same outlet's diameter of 10mm and therefore the maximum value is achieved later. The smaller chamber in the case of D = 35mm with the same outlet's diameter of 10mm leaves the flow faster (the time lag between the inlets and outlet is shorter than the case of D = 40mm ) and does not let to generate an intensive dead zones inside the JSR, therefore it provides a well-stirred mixture. The normalized CO 2 standard deviation of the case of a nozzle diameter of d = 1.5mm is compared to d = 1mm for a constant chamber diameter of D = 40mm in Fig. 11b. The result shows that the larger diameter of the nozzle cannot improve the mixing inside the JSR. The reason can be understood as follows. In order to keep the theoretical mean residence time constant, the volumetric inflow rate is also kept constant. Therefore, the initial velocity modeled at the entry of the nozzles in the case of d = 1.5mm is less than the initial velocity in the case of d = 1mm . Smaller velocity for the case of larger nozzle diameter of d = 1.5mm generates a more intensive dead zone at the center of the JSR as illustrated in Fig. 12. In the dead regions a very small exchange of material happens inside the reactor with insufficient molecular mixing [41].

Proposed improved geometry
In order to improve the mixing quality, we added two injectors of a diameter of d = 1mm to the standard four nozzle design introduced by Dagaut et al. [11]. An overview of this suggested setup having six nozzles is given in Fig 13. One of the advantages of this geometry is that the two nozzles pointing to the center breakdown the dead-zones at the center of the JSR which, based on the discussion above, are mainly responsible for the nonhomogeneity. Figures 14a and b compare two cross sections of the time-averaged tracer mass fraction for cases of four and six nozzles in the x-z and x-y planes, respectively.
Comparing the values of the averaged tracer mass fraction over three points shows that the amount of non-homogeneity in the case of six nozzles is decreased compared to the case of four nozzles. In fact, the dead-zone is removed from the reactor by applying two more jets pointing to the center.
Comparing the resulting curves for the normalized standard deviation of the CO 2 mass fraction in Fig. 15 indicates that a more uniform composition is achieved for the new design of the JSR compared to the typical JSR. The lower normalized standard deviation in the case of six nozzles might be explained by the vortexes at the center of the JSR induced by the high inertia of the jets.  We are aware that the JSR design of Dagaut [11] requires a premixed chamber of the gases before they enter the reactor spherical chamber; however, adding two nozzles will improve the mixing inside the reactor. The time-averaged streamlines using LIC technique and resolved kinetic energy are further discussed in the appendix A.

Conclusions
The mixing process and spatial homogeneity inside a spherical Jet-Stirred Reactor (JSR) have been studied via well-resolved large eddy simulations. First, the grid resolution study is investigated and then the characteristic of flow field created by four impinging jets have been captured through a numerical method. Analyzing stream-lines and velocity distributions with LIC method introduced the details of the flow field pattern which could not be obtained by limited optical accessibility of experiments. Further, studying the time-averaged distribution of the resolved specific turbulent kinetic energy and vorticity shows that the highest turbulence fluctuations and, therefore, best mixing occur where the jets impinge on the inner surface of the vessel. Results of the energy spectrum following −5∕3 law for a wide range of frequencies represents a state of fully developed turbulence. Second, a parametric study in terms of geometrical parameters and thermodynamic conditions has been performed in order to improve mixing inside the JSR. According to the results of the normalized standard deviation of the tracer mass fraction and the resolved specific turbulent kinetic energy, higher temperatures might not follow the ideal reactor. However, increasing the pressure above atmospheric results in a lower value of the normalized standard deviation of CO 2 which indicates an improvement to the mixing. Therefore, the four-jets JSR with nozzle diameters of 1mm and chamber diameter 40mm at atmospheric pressure up to 10bar create a nearly homogeneous composition. Finally, we proposed a new JSR design by adding two nozzles pointing to the center of the reactor. Our simulations predicted that this design would breakdown the dead-zones at the center of the JSR and significantly improve the mixing inside the reactor.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Appendix A: LIC method
Illustration of streamline using LIC Figure 16 displays the time-averaged velocity contour of the rotating confined flow over the time of [ 0.1s−0.4s] after the injection of CO 2 . The streamlines shown in Fig. 16 are extracted using the surface-LIC (Line Integral Convolution) technique which was introduced by Cabral et al. [30] and was subsequently developed by Loring et al. [31]. This technique uses an image vector field algorithm which is gradually extended in several software packages, e.g. paraview [31], to visualize the flow traces efficiently. The obvious lack of vortices in the center of reactor motivates the development of an optimal configuration of multiple pairs of impinging jets.

Comparing the LIC contours in different pressures
The time-averaged velocity contours captured with the LIC method compares the flow structures close to the outlet for a pressure of 1bar in Fig. 17a and pressure of 10bar in Fig. 17b. According to this figure, three vortexes appeared over the reactor in the case of high pressure which improve the mixing. The LIC visualization illustrated in Fig. 18a shows that these vortexes intensify the turbulent structures over the reactor which leads to an efficient mixing. This suspicion is further supported by the resolved specific turbulent kinetic energy field in Fig. 18b that exhibits a high resolved kinetic energy in the center of the reactor.