Abstract
A method is presented to significantly improve the convergence behavior of batch nonpremixed counterflow flame simulations with finite-rate chemistry. The method is applicable to simulations with varying pressure or strain rate, as it is, for example, necessary for the creation of flamelet tables or the computation of the extinction point. The improvement is achieved by estimating the solution beforehand. The underlying scaling rules are derived from theory, literature, and empirical observations. The estimate is used as an initialization for the actual solver. This enhancement leads to a significantly improved robustness and acceleration of batch simulations. The extinction point can be simulated without cumbersome code extensions. The method is demonstrated on two test cases and the impact is discussed.
1. Introduction
Laminar steady nonpremixed counterflow flames have been investigated in numerous studies in the past [1–3]. Due to their quasi-one-dimensional nature, they are simulated on a computationally very inexpensive grid. Extensive chemical mechanisms can be applied, tested, and investigated. The results can be validated well against state-of-the-art experiments [1].
From counterflow flame simulations, detailed insight is gained into the structure of nonpremixed flames, including temperature and species profiles [3], heat release [3, 4], and flame radiation [5]. This setup is often used as a test bed for novel or reduced combustion mechanisms. Another purpose is the tabulation of flames in the form of flamelet libraries for use in turbulent combustion simulation [6].
Even though counterflow flame simulations usually run on a very small grid, the computation is time consuming because of the stiff equation system resulting from the detailed chemical mechanism. If the initialization is too far from the steady state solution, the solvers are prone to diverge.
A major aim of counterflow flame simulations is to compute extinction strain rates [7]. In this case, a second problem arises. At the extinction point, the numerical solution of a regular counterflow flame simulation becomes singular [2, 7]. The exact extinction flamelet therefore cannot be computed. To circumvent this problem, two sophisticated methods have been developed and applied on counterflow flames: the arc-length continuation method by Keller [8] and Giovangigli and Smooke [9] and the flame controlling continuation method by Nishioka et al. [10]. The idea behind both methods is to compute the unstable burning branch, that is, to calculate pseudostable flames which show maximum temperatures and strain rates below the extinction values. The extinction point can thus distinctly be identified. To realize this, the computational procedure has to be altered severely. A recent implementation is well documented by Wang et al. [7].
For both the computation of extinction strain rates and the creation of flamelet libraries, multiple simulations at varying strain rates are necessary. The operating pressure is often an investigated parameter as well. This can result in a large matrix of input parameters for which separate simulations are required. While the computational effort might still be acceptable for single simulations, it can be very large for the overall simulation campaign. A reduction of the computational effort is therefore favorable.
This paper develops a method to stabilize and accelerate batch counterflow flame simulations. These goals are achieved by improving the initialization for the iterative steps. Depending on the pressure and the strain rate changes, estimations for the change of the velocity, species, and temperature profiles are derived from the literature, theoretical considerations, and empirical observations.
2. Method
To properly deduce the scaling rules, the setup of the counterflow flame and its definitions are briefly described first. The approach to improve the batch computations is then derived theoretically and finally put into an applicable form.
2.1. Counterflow Flame Setup
A schematic view of the considered counterflow flame configuration is shown in Figure 1. Two axis-symmetric jets of fuel and oxidizer impinge on each other from opposed nozzles, thus creating an axis-symmetric flow field which can be described by the axial coordinate and the radial coordinate . Only the axial coordinate has to be resolved to specify the flame [2].
The well-known governing partial differential equations are described by Kee et al. [2]. The problem can be solved numerically by evaluating the continuity equation, the radial momentum equation, a species conservation equation for each species considered, and the energy equation. Species source terms are calculated from detailed chemistry reaction mechanisms. Pressure, density, and temperature are linked by the ideal gas equation of state. Real gas effects were shown to have little influence on the structure of counterflow flames [7, 14].
The flow field is fully characterized by the axial coordinate , the axial velocity , the radial velocity divided by the radial coordinate , the radial pressure curvature (which is an eigenvalue and constant throughout the flame), the temperature , and the species mass fractions for all considered species.
So-called plug-flow boundary conditions are applied on both the oxidizer and the fuel inlet [2, 3]. The mass flow rate area densities and , the inlet temperatures and , and the inlet compositions and are specified on the fuel and oxidizer boundaries. The distance between fuel and oxidizer inlets is defined as the grid size .
The counterflow flame is defined by setting the boundary conditions, the grid size, the operating pressure, and a chemical reaction mechanism. Counterflow flames are commonly characterized by their strain rate . However, different definitions of the strain rate are used in the literature: can be determined from the mean axial velocity gradient [3], the maximum axial velocity gradient [15], the axial velocity gradient at the stoichiometric point [16], or the axial velocity gradient “before the flame” [3, 4, 6]. The latter corresponds to the boundary conditions in the potential flow setup where fuel and oxidizer are assumed to be infinitely far away. Sometimes, radial velocity gradients are used instead [2]. If oxidizer and fuel do not have equal density, the strain rates on either side of the flame are different [3]. All of these definitions are approximately directly proportional to each other. Nevertheless, attention has to be paid when comparing absolute results. In this paper, the maximum axial velocity gradient is used to characterize the strain rate.
In the plug-flow configuration, the potential flow boundary conditions can be approximated by for the fuel and for the oxidizer sides, respectively (this result can be obtained by reshaping from [17] or from [2]).
With all parameters defined, counterflow flames can be computed with standard software like OPPDIFF [18], Cosilab [19], FlameMaster [20], or Cantera [21]. These programs solve the underlying partial differential equations of the one-dimensional problem. To fully resolve the flame structure, an automatic grid refinement is often necessary. The solution procedure can take several minutes because of the stiff system of equations to be solved. The further away the final solution is from the initial solution, the larger the computation time is from the initial solution. If the initial solution is too far from the actual solution, the solver can easily diverge.
2.2. Approach
As pointed out in the introduction, for the calculation of extinction strain rates or the generation of flamelet tables, a large number of counterflow flames have to be simulated. This is generally done by performing a batch simulation.
The idea of improving the convergence behavior of batch simulations is to estimate the solution beforehand. This guess is used as the initialization for the actual computation. The estimation is based on the assumption that the structure of flames with finite-rate chemistry behaves similarly to the structure of flames with infinitely fast chemistry. In nonpremixed flames, the structure is mostly influenced by diffusion which is equally represented in flames with finite-rate and infinitely fast chemistry.
Counterflow diffusion flames with infinitely fast chemistry are known to have a self-similar structure [17]. According to Law [3], the spatial profiles in nonpremixed counterflow flames scale with the term If the inlet temperatures are considered to be constant, the density is directly proportional to the pressure . Assuming that the thermal conductivity and the heat capacity at constant pressure are independent of pressure and strain rate, the spatial profiles are directly proportional to the term .
This formulation is equivalent to the postulation that the diffusion flame thickness (here, full width at half maximum (FWHM) of the temperature profile [3]) , which corresponds to the thickness of the mixing layer, scales with [16, 17] If the diffusion coefficient is assumed to be inversely proportional to the pressure, the above relation can again be transformed to This relation has been confirmed numerically to hold for a wide range of counterflow flames with detailed chemistry [4, 22].
2.3. Application to Batch Computation
To apply this on batch simulations, an existing solution with the strain rate and the pressure is considered. This flame has the known thickness . The thickness of a new flame to be calculated at and can be estimated according to (4): Since all species profiles and the temperature profile are assumed similar with respect to the flame thickness, they can be scaled accordingly.
However, instead of scaling the flame on a fixed grid, the grid size can be modified: This procedure has the advantage that the temperature and the species mass fractions remain constant with respect to the grid points and do not have to be scaled separately. The magnitudes of both the temperature and the species mass fractions are supposed to change little for a batch step. Therefore, the absolute values of these variables can remain unaltered. Additionally, a potential refinement of the grid is adapted as well automatically.
In contrast to the temperature and species profiles, the magnitudes of the velocity profiles, the pressure curvature, and the boundary conditions have to be updated. Using the definition of the mean strain rate and the fact that all strain rate definitions are approximately directly proportional to each other, the axial velocity profile can be assumed proportional to By substituting this relation into the continuity equation [21], a scaling for can be derived: This leads to If the ideal gas densities are assumed to scale with the pressure, the mass flux densities at the boundaries adapt accordingly: Following (1), the pressure curvature can be estimated:
The equations derived above are summarized in Table 1. In addition to the general case, the scaling factors at constant pressure and are attached. The constant pressure case is important for the computation of isobaric flames, for example, for the computation of the extinction strain rate at a specific pressure or the creation of a flamelet table. The last column is added for the case that different pressures should be calculated. With rising pressure, the extinction strain rate increases. If a batch of counterflow flames up to extinction is to be computed for various pressures, starting each pressure batch at the same initial strain rate would lead to an increasing number of simulations with pressure. To compute approximately the same number of simulations for all pressures, the initial strain rate has to be modified. The relation was found empirically to be a good scaling for the initial strain rates. This is based on the analysis of hydrogen-oxygen and ethane-air flames. The exact quantity probably depends on the fuel and oxidizer compositions, but is recommended as a good starting guess for such a scenario.
3. Results and Discussion
The method described above is implemented on top of the Cantera environment (version 2.2a, revision r2777 [21]). The simulation is set up using the counterflow flame class of the Cython-based interface. The batch procedure is formulated as a Python script. In between the simulation runs, the new initial condition is estimated using the relations derived above.
The two benefits of the proposed method, the speed-up and the robustness, are demonstrated using two test cases described in the following.
3.1. Speed-Up Test Case
To indicate the speed-up achieved by the proposed method, a series of counterflow flames at various pressures and strain rates is simulated. This could be the basis for the creation of a flamelet table. Since the computation of the underlying counterflow flames is the most computationally expensive step in the creation of a flamelet library, such a task is assumed to speed up accordingly. The acceleration is shown by comparing a batch simulation which incorporates the proposed method with a state-of-the-art benchmark simulation. This is demonstrated for a hydrogen-oxygen flame and an ethane-air flame. The overview over the simulation parameters is given in Table 2.
Starting from initial counterflow flames at a pressure of 1 bar, flames up to 100 bar should be calculated. The pressure is increased by approximately 15% in each step. Additionally, at 13 of the previously computed pressure levels, isobaric flames should be calculated at increasing strain rates. The strain rates are raised by 25% in each step until the flame is extinguished.
For each fuel-oxidizer configuration, the test case is run two times. The accelerated run makes use of the procedure described in this paper. In the corresponding benchmark run, the initialization procedure is not applied. The pressure and strain rate are set in this case by changing only the operating pressure and the inlet mass fluxes, respectively. Both cases are run sequentially on a standard desktop PC (housing an Intel i5 CPU at 3.2 GHz) under standard load conditions.
For the hydrogen-oxygen flame, the benchmark computation fails to converge for pressures above 64 bar. At this point, the cumulated simulation time already amounts to 14:12 minutes. The accelerated run runs smoothly for all pressure levels; at 64 bar the cumulated simulation time is only 43 seconds. For the ethane-air flame, both variants are able to compute all pressure levels. The simulation time for the entire loop is accelerated from 2:47:35 hours to 15:15 minutes using the proposed method. The values are compared in Table 3.
The loops through the increasing strain rates profit even more from the preconditioning step. For the hydrogen-oxygen flame, the average time for an isobaric loop decreases from 1:49 minutes to only 5 seconds. For the ethane-air flame, the corresponding times are 3:55:29 hours and 4:08 minutes. At pressures above 10 bar for the hydrogen-oxygen flame and above 30 bar for the ethane-air flame, the solver diverges in the strain rate loop in the benchmark run. This problem does not occur for the accelerated run. This fact also indicates the improvement in robustness achieved by using the derived scaling rules.
Figures 2 and 3 show the axial velocity and temperature profiles over the normalized grid for the accelerated case of the hydrogen-oxygen flame. It can be seen that the flame thickness with respect to the grid size indeed remains constant. The axial velocity profiles also remain similar. This proves that the scaling rules derived from infinitely fast chemistry in fact apply to finite-rate chemistry counterflow flames.
3.2. Simulation of the Extinction Point
A second independent test case is evaluated to demonstrate the robustness increase achieved by this method. This makes it possible to compute the extinction point of counterflow flames very accurately without having to apply a cumbersome method to circumvent the singularity at this point.
For multiple flames at constant pressure levels, the strain rate is increased gently until the flame is extinguished. From the last burning flame, the strain rate is increased in smaller steps to get closer to the extinction point. This iteration is repeated until the temperature change is below a threshold value. The strain rate and the maximum temperature of the last burning solutions are taken as the extinction strain rate and extinction temperature.
This is demonstrated for a hydrogen-oxygen flame. The boundary and operating conditions of the initial simulation are given in Table 2. The output strain rates and temperatures at the extinction point are shown in Figure 4. The data are compared to the data by Wang et al. [7]. The same reaction mechanism and the boundary conditions were used as in the reference paper. While the extinction strain rate matches the reference values almost perfectly, the maximum flame temperature at extinction is slightly overpredicted. This is due to the fact that the reference data were created using a sophisticated two-point flame controlling continuation method. The method presented in this paper approaches the extinction point from higher temperatures.
4. Conclusions
In this paper, a method is described which significantly improves the convergence behavior of batch counterflow flame simulations with varying pressure and strain rate. This is achieved by estimating the initial solutions from previously computed simulations. The scaling parameters for the solution vectors with respect to the pressure and the strain rate are derived from existing relations for flames with infinitely fast chemistry, dimensional analysis, and empirical observations.
The method is demonstrated using the Cantera software package [21] on a nonpremixed hydrogen-oxygen and an ethane-air counterflow flame. By computing batch simulations with and without applying the proposed routine, the calculation speed is observed to increase by a factor of 10–20 for a variation in pressure and a factor of 20–60 for variations in strain rate. Additionally, using the improved procedure, several operating points can be computed that fail to converge in the reference run.
The method is also used to compute the extinction point of a hydrogen-oxygen counterflow flame. The result very accurately matches data from the literature. Using the described simple procedure, the extinction point can be calculated without having to use cumbersome continuation methods.
Nomenclature
| : | Diffusion coefficient [m2/s] |
| : | Temperature [K] |
| : | Radial velocity potential [1/s] |
| : | Mass fraction [—] |
| : | Strain rate [1/s] |
| : | Isobaric heat capacity [J/(kg K)] |
| : | Grid size [m] |
| : | Mass flux density [kg/(s m2)] |
| : | Pressure [Pa] |
| : | Axial velocity [m/s] |
| : | Axial coordinate [m] |
| : | Pressure curvature [N/m4] |
| : | Flame thickness [m] |
| : | Heat conductivity [W/(m K)] |
| : | Density [m3/kg] |
| : | Fuel inlet |
| : | Oxidizer inlet. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The Deutsche Forschungsgemeinschaft (DFG) in the framework of Sonderforschungsbereich Transregio 40 provided financial support for this project. This work was supported by the German Research Foundation (DFG) and Technische Universität München within the funding programme Open Access Publishing.