An open-source tool for predictive simulation of diffusion flames in analytical chemistry

This work presents an open-source computational model of gas flow and hydrogen combustion in a miniature diffusion flame atomizer.


Introduction
Hydride generation (HG) combined with atomic uorescence (AF) or atomic absorption (AA) spectrometry for determination of hydride forming elements is a viable alternative to the conventional approaches based on liquid phase sampling inductively coupled plasma mass spectrometry. 1 The most oen employed hydride atomizer for AF is a miniature diffusion ame (MDF) 2-7 but it can also be very useful for AA spectrometry. 4,8 An optimization of the atomization is required to reach an excellent performance of the whole procedure of element determination based on HG. Presently, the optimization of hydride atomization must be performed by the laborious trialand-error approach. Obviously, it could be made in a straightforward and elegant manner based on the knowledge about what really happens in hydride atomizers. The theory describing what happens in atomizersthe radical theory of hydride atomization 4,[9][10][11] is based on the formation of free hydrogen atoms (H radicals) at a concentration several orders of magnitude above the equilibrium. Hydrides are then atomized by interaction with H radicals, so their distribution controls what happens in atomizers. A straightforward optimization of hydride atomizers thus requires knowledge of the distribution of H radicals for given atomization parameters. In principle, two-photon absorption laser-induced uorescence (TALIF) can be employed to determine the distribution of H radicals in atomizers; 12 however, a computation quantitative simulation of a MDF could perform this in a revolutionary manner. A predictive computational model allows the researcher to assess the performance of many geometrical congurations and experimental conditions at a much faster rate than TALIF experiments.
The aim of this was to work was to prove that a predictive computational model of the processes in a MDF can facilitate much easier optimization of the atomizer. To make this effort easier for the community, we do not only describe the physics and chemistry model of the atomizer but we also provide executable model les under an open-source academic license, so that they can be further utilized by the scientic community. Since the model is rather general it is also applicable to the optimization of the design and operating parameters of other hydride atomizers 4 employed for AF and AA spectrometers, such as ames in gas shield atomizers, conventional externally heated quartz tube atomizers or multiatomizers.
Furthermore, we assess the sensitivity of the model results to uncertain experimental parameters (laboratory conditions, tube diameter and air contamination). This analysis provides insight into the uncertainty of both the model and the experiment and reveals factors to which the distribution of H radicals is most sensitive.
Finally, this work conrms that the GRI-Mech 3.0 reaction system, 13 which was originally developed for natural gas combustion, can also be used for reliable predictive modeling of the combustion kinetics in hydrogen non-premixed ames. As such, this publication extends the list of studies that have experimentally validated the GRI-Mech 3.0 reaction system for various types of amessee ref. 13  Numerical simulations similar to this work have been previously presented in the literature, 14 where the fundamentals of hydrogen-air coow ames were studied. Compared to ref.
14, this work offers model validation which is immediately relevant for analytical chemistry applications and, above all, it provides the model itself as a scientic instrument. To our knowledge, this is the rst openly available numerical model of diffusion co-ow ames.

Computational method
The computational model self-consistently solves the continuity equations for the mass fractions of all the species included in the model Y i , coupled to the momentum equation for the velocity of the entire mixture u and the energy equation for the temperature of the gas mixture T.
The model was not made from scratch but rather implemented using the open-source computational framework called laminarSMOKE+. 15 Unlike the freely available reactingFoam solver and the commercially available ANSYS Fluent solvers, 16 which have been designed for the more common case of turbulent combustion, laminarSMOKE+ is much more suited for building models of diffusion ames because it uses binary coefficients and mixture averaging for calculating the diffusion coefficients.

Governing equations
The set of governing equations and the approach to their discretization and solution in laminarSMOKE+ is described in great detail in ref. 15. The set of equations and the solution approach are well known in the combustion community but we list them here nevertheless, because they might not be so well known to the analytical chemistry community. We also provide a more detailed commentary on the meaning of the individual equations, which is typically considered obvious in combustion modeling papers. The presented set of equations also differs from that of the original work 15 upon omitting the terms that were not included in the diffusion ame model. Firstly, the model solves the continuity equation for the mass of the gaseous mixture in the form where r is the mixture mass density and v is the mixture's mean velocity. In other words, the model does not solve for the velocities of the individual components and it replaces them with the velocity of the center-of-mass of the mixture. The density r in this equation is then not only the function of temperature T and pressure p, but also of the local gas composition, expressed by mass fractions The momentum equation for the mixture velocity then takes the form: Again, since the equation above is formulated for the entire gas mixture, the stress tensor T depends on the local gas composition and mixture temperature. The vector g is the acceleration vector due to gravity.
The mixture temperature is obtained from the energy equation of the mixture in the form: where C p is the mixture heat capacity at constant pressure, C p,k are individual components' heat capacities at constant pressure, h k is the sum of formation and sensible enthalpies for species k and U k is the formation rate of species k. The vector q is the heat ux expressed through the thermal conductivity l as q ¼ ÀlVT and V k is the diffusion velocity of species k calculated using the species diffusion coefficient D k using the Fick expression as

Finally, a diffusion equation is solved for each species in the form
where the second term on the le-hand side corresponds to species transport due to advection (bulk ow) and the third term corresponds to diffusive transport. Unlike the abovementioned reactingFoam solver, the diffusion coefficients in the laminarSMOKE solver are calculated using mixture-rules 17,18 from binary diffusion coefficients G jk , which is currently the most accurate approach for diffusiondriven combustion. In the mixture-averaged approach, the binary diffusion coefficients are converted to Fick diffusion coefficients by calculating a weighted average according to the formula where X k is the mole fraction of species k. This way, the Fick diffusion coefficient for each species becomes a function of not only temperature, but also the local gas composition. This includes strong non-linearity into the set of equations, but the main advantage of this approach is the possibility to simulate the movement of the mixture as a single uid with the mean dri velocity v.

Species and the reaction set
The aim of the presented model is to simulate argon-diluted hydrogen ames in air. As mentioned in the Introduction, the reaction set in this work is the subset of the GRI-MECH 3.0 reaction mechanism, 13 from which we selected only the reactions which include the species mentioned above. The sub-set of reactions is provided in the ESI † in the CHEMKIN format and is also available in the open source model repository (see Section 2.4 below).
We did not list the reactions and their rates in this paper. Instead, we included both the reaction mechanisms in a CHEMKIN-compatible format as the ESI † to this publication.

Computational geometry and boundary conditions
The computational geometry is a 2D axially symmetrical one, shown in Fig. 1. It corresponds to the geometry of the employed MDF designsee the Experimental section. The computational geometry is discretized using a structured rectangular mesh.
There is a different boundary condition imposed for each boundary. For the inlet boundary, we prescribe the mean velocity of the ow, calculated from the ow rate and tube diameter. We also impose the Dirichlet boundary condition for the mass fractions of all the speciesin the default case without impurities, all mass fractions are set to zero, except for argon and hydrogen.
For the wall boundary, we prescribe a zero-gradient boundary condition for the species' number densities and a no-slip boundary condition for mixture velocity. By assuming zero-gradient at the walls, we effectively neglect quenching of reactive species by surface reactions. This is a reasonable assumption because the area of the wall which is in contact with the ame is small.
The air intake boundary condition serves as a source of fresh air (oxidizer) in the model. Therefore, we set the mass fractions of all the species except for O 2 , N 2 and H 2 O to zero. Regarding velocity, we impose a constant ow velocity of 0.05 m s À1 on the air intake, to provide a steady supply of air. In the experimental system, there is no forced ow outside the MDF tube, but the air circulates around it due to the buoyancy force. It is, however, not practical to assume a zero-gradient boundary condition for velocity on the intake because the system of equations would not have a unique solution.
Finally, the outow boundary condition assumes a zerogradient for both the species mole fractions and for the gas mixture velocity, acting as a free outow for both these quantities.

Model distribution
As mentioned in the Introduction, the model described in this work is provided to the community in the form of conguration les that can be executed using the laminarSMOKE solver. The model les are provided as the ESI † and also in a public git repository located at: https://bitbucket.org/adamobrusnik/mdfatomizer-laminarsmoke.

Experimental
The model was benchmarked against the MDF atomizer. The simple design of this device reduces the number of experimental uncertainties and imperfections.
Unlike many studies focusing on combustion simulation, which chooses the simulated and experimental ame radius or length and ame velocities as the main observable, 19-21 this work goes into greater detail by using the H radical 2D maps as the observable. The rationale behind this is that the distribution of H radicals has a determining inuence on an atomizer's performance, as detailed in the Introduction.

Miniature diffusion ame atomizer
The MDF atomizer described e.g. in ref. 22 is a simple device consisting of a quartz glass tube (quartz, inner diameter 6 mm, outer diameter 8.5 mm), through which a mixture of hydrogen (99.95%, SIAD Ltd., Czech Republic) and argon (99.996%, SIAD Ltd., Czech Republic) was fed to ambient air. The ame was ignited at the end of the tube, where the hydrogen/argon mixture mixes with ambient air. The gas ow rates were controlled by using mass ow controllers (FMA 2400 Series, Omega Engineering, Inc., USA). The gases were pre-mixed before entering the quartz tube.

Laser-induced uorescence measurement
In order to measure the distribution and concentration of hydrogen radicals, 23,24 they were excited by 8 ns long laser pulses with a wavelength of 205 nm and a spectral width of 0.06 cm À1 . The laser pulses were produced with a 30 Hz repetition rate by using a laser system consisting of a Q-switched pumping laser (Quanta-Ray PRO-270-30), a dye laser (Sirah, PRSC-D-24-EG) and a frequency tripling unit. Laser radiation was focused on the center of the studied atomizer by using a spherical silica lens (focal length approx. 50 cm). The energy of each pulse (typically hundreds of mJ) was detected by using a pyroelectric energy sensor (Ophir PE9). The hydrogen uorescence (H a , 656 nm) was recorded by using an intensied CCD camera (PI-MAX 1024RB-25-FG43) with an interference lter. The integration time of the ICCD was 10 s in all the experiments.
In order to obtain the absolute concentration of atomic hydrogen, the sensitivity of the detection was calibrated by measurement of Kr TALIF. 25,26 Kr atoms with a known Fig. 1 Computational geometry of the numerical model with dimensions and annotated boundaries (a) and the mesh used for discretization of the partial differential equations (b). concentration (1.5-8%) in Ar gas were excited by using a laser with a wavelength of 204 nm and the generated Kr uorescence was detected. Then, the atomic hydrogen concentration was calculated by means of the formula N H ¼ N Kr (S H /S Kr )(E Kr /E H ) 2 (n H /n Kr ) 2 (s (2) Kr /s (2) H )(q Kr /q H )(K Kr /K H ), where N denotes concentration, S the measured TALIF signal that was generated by laser pulses with energy E and that was Fig. 2 . Steady-state spatial distribution of species' number densities (plotted in log 10 scale), temperature and reaction heat source term in the miniature diffusion flame at Q Ar ¼ 700 sccm and Q H 2 ¼ 300 sccm. The quantity gas_Q is the heat released by the chemical reactions per unit time.
spectrally integrated over the whole excitation line prole, s (2) two-photon excitation cross section, q quantum efficiency of the uorescence and K the sensitivity of the detection system for the uorescence wavelength. The indexes H and Kr differentiate between the measurements of atomic hydrogen and krypton. The details of the method are described in ref. 27 and 28. 4 Results and discussion

Species' spatial distribution
The spatial distribution of the different species in the MDF atomizer is shown in Fig. 2. As expected, the combustion of hydrogen and oxygen occurs in the interfacial region, where the argon/hydrogen mixture mixes with ambient air. This is illustrated especially by the subplots gas_Q, showing the heat released by combustion and T being the gas temperature. As expected, the most reactive H and OH species are also present only in the combustion region.

Experimental validation
We benchmarked the model against the spatially resolved density of H radicals measured by TALIF. It is a good measure of the ame's properties because H radicals will only be present in the active combustion zone. Most importantly, according to the radical theory of hydride atomization the distribution of H radicals controls what happens in atomizers (see Introduction). Two types of data sets have been obtained experimentallya 1D radially resolved measurement of the H radical number density which has been absolutely calibrated and 2D spatially resolved maps of H radicals. With the 2D data, the quantum yield of the H uorescence could not have been measured with sufficient accuracy, so we only plot the data in arbitrary units. The 1D radially resolved data allow us to quantitatively compare the model predictions with the experiment. The data are plotted in Fig. 3 for 2 mm above the top of the atomizer tube. It is apparent that, in terms of absolute values, the number density of H is within the same range. The experimental data are, however, more diffuse compared to the simulation data. This is most likely caused by the jittering of the ame, which effectively changes the position over the 10 s integration time of the ICCD camera.
The more diffuse nature of the experimentally measured data is also apparent in the 2D spatially resolved data, shown in Fig. 4. However, the model tends to capture very accurately the  and from the model (right half) obtained for various argon/hydrogen mixtures. The experimental data lack absolute calibration and the scale has been modified to match the simulation. The maximum value of the H number density was obtained from the simulation. Note that we take advantage of the axial symmetry of the flame, so r < 0 shows the experimental data while r > 0 shows the simulation data.
length of the ame and the overall shape of the H-rich regions in several ratios of argon/hydrogen supplied to the atomizer. Apparently, H radicals are present especially in the interfacial mixing region, where the argon/hydrogen mixture supplied to the atomizer mixes with ambient air.
The model and experiment consistently show that the length of the ame can be controlled by the argon ow ratecompare Fig. 4(b) and (c). This is because increasing the argon ow rate effectively increases the ow velocity and the combustion of hydrogen occurs over a greater distance. At the same time, however, increasing the argon ow rate decreases the absolute density of H radicals. On the other hand, increasing the fraction of H 2 (under the same ow rate of the Ar/H 2 mixture) does increase the ame length and H radical densitycompare Fig. 4(a) and (c).

Sensitivity to experimental uncertainties
To model an actual experimental setup, one must take into account all the manufacturing and experimental uncertainties that could inuence the simulation results. In this section, we quantify how these uncertainties inuence the simulation result. We have chosen Q Ar ¼ 700 sccm and Q H 2 ¼ 300 sccm as the default working point with a laboratory air humidity of 40% at 300 K and 0 ppm air admixture in the atomizer. In accord with the dimensions of the MDF, the inner diameter of the quartz tube was set to 6 mm. However, in various reproductions of the experiment, these values might not always be maintained. Therefore, we have carried out a series of simulations to assess the sensitivity of the model to the key parameters. These simulations are summarized in Table 1. The variations to the model parameters correspond to actual instrument or manufacturing tolerances.
Two observables were chosen for quantifying the uncertainty. One is the absolute yield of H radicals produced by the atomizer Y H ¼ Ð V n H ðrÞdVi.e. the number density of hydrogen radicals integrated over the simulation domain. The second observable is the length of the ame L f which is dened as the position on the z-axis, where the mass fraction of H 2 drops below 0.1% of the maximum value. In both cases, we plot the change in these observables with respect to the default simulationthe value of DL f ¼ +10% corresponds to a 10% increase of the ame length w.r.t the default simulation.
Firstly, in Fig. 5, we plot the sensitivity of the total hydrogen yield. Apparently, the main sources of uncertainty in the presented setup are the tube diameter. Having a tube with a diameter increased by 0.4 mm reduces the total H yield by approximately 10% while decreasing the tube diameter increases the H yield by approximately 6%. The humidity of ambient air plays only a small role in the total H yield, with the H yield decreasing with increasing humidity. This result is reasonable, considering that H 2 O molecules will deplete some of the H radicals, ultimately converting them to H 2 O 2 . Finally, the 0.1% air impurity in the gas inlet has only minor impact on the total H yield.
In Fig. 6, we plot the sensitivity of the simulated ame length to experimental uncertainties. In this case, only the tube diameter has non-negligible inuence within the range of experimental uncertainties specied in Table 1. The result is to be expected as a larger tube diameter decreases the initial velocity of the gas ow while a smaller diameter increases it, thereby stretching the ame longer.

Conclusions
In this publication, we have introduced an open-source computational model of non-premixed hydrogen combustion in atomizers. The model is based on an open-source (for academic use) solver laminarSMOKE and the reaction mechanism used is a sub-set of the GRI-Mech 3.0 reaction mechanism.
The model developed has been used for simulating H radical production in a simple MDF atomizer. By comparing the simulated H density prole with TALIF measurements of the   same quantity, we conrmed that the predictions of the model are correct, both in terms of the shape of the ame and in terms of the absolute values of H radical concentrations. The main difference between the simulation and the experimental data lies in the fact that the experiment shows H density proles which are more diffuse. However, this observation may be attributed both to the imperfections of the model (underestimated diffusion) and experimental imperfections (moving or oscillatory ame position).
We have also used the model to assess the sensitivity of the MDF atomizer to uncertain laboratory conditions and geometrical parameters. The ability to quantify how laboratory conditions and other experimental parameters inuence the performance of an atomizer is another example of how the computational model can be utilized.
Since the H radical distribution controls what happens in atomizers, it has become obvious that the presented model is a powerful predictive tool that can be utilized for the design and optimization of not only MDF atomizers but other atomizer congurations as well. The model supports both 2D axially symmetrical setups, e.g. ames in gas shield atomizers, and setups which require full 3D geometry, e.g. conventional externally heated quartz tube atomizers or multiatomizers. 3 By providing the model simulation les under an opensource license, we allow and invite other members of the analytical chemistry community to modify the geometry or operating conditions and use the tool for their research.

Conflicts of interest
There are no conicts to declare.