A fast, low-memory, and stable algorithm for implementing multicomponent transport in direct numerical simulations
Aaron J. Fillo, Jason Schlup, Guillaume Beardsell, Guillaume Blanquart, Kyle E. Niemeyer
AA fast, low-memory, and stable algorithm forimplementing multicomponent transport in directnumerical simulations
Aaron J. Fillo a , Jason Schlup b , Guillaume Beardsell c , Guillaume Blanquart c ,Kyle E. Niemeyer a, ∗ a School of Mechanical, Industrial, and Manufacturing Engineering, Oregon StateUniversity, Corvallis, OR 97331, USA b Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA USA c Department of Mechanical Engineering, California Institute of Technology, Pasadena, CA,USA
Abstract
Implementing multicomponent diffusion models in reacting-flow simulations iscomputationally expensive due to the challenges involved in calculating diffu-sion coefficients. Instead, mixture-averaged diffusion treatments are typicallyused to avoid these costs. However, to our knowledge, the accuracy and ap-propriateness of the mixture-averaged diffusion models has not been verified forthree-dimensional turbulent premixed flames. In this study we propose a fast,efficient, low-memory algorithm and use that to evaluate the role of multicom-ponent mass diffusion in reacting-flow simulations. Direct numerical simulationof these flames is performed by implementing the Stefan–Maxwell equations inNGA. A semi-implicit algorithm decreases the computational expense of invert-ing the full multicomponent ordinary diffusion array while maintaining accu-racy and fidelity. We first verify the method by performing one-dimensionalsimulations of premixed hydrogen flames and compare with matching cases inCantera. We demonstrate the algorithm to be stable, and its performance scalesapproximately with the number of species squared. Then, as an initial studyof multicomponent diffusion, we simulate premixed, three-dimensional turbu-lent hydrogen flames, neglecting secondary Soret and Dufour effects. Simula-tion conditions are carefully selected to match previously published results andensure valid comparison. Our results show that using the mixture-averaged dif-fusion assumption leads to a 15 % under-prediction of the normalized turbulentflame speed for a premixed hydrogen-air flame. This difference in the turbulentflame speed motivates further study into using the mixture-averaged diffusionassumption for DNS of moderate-to-high Karlovitz number flames.
Keywords:
Turbulent flames, Direct numerical simulation, Multicomponent ∗ Corresponding author
Email address: [email protected] (Kyle E. Niemeyer)
Preprint submitted to Journal of Computational Physics November 7, 2019 a r X i v : . [ phy s i c s . f l u - dyn ] N ov iffusion, Mixture-averaged diffusion
1. Introduction
Implementing full multicomponent mass diffusion transport in direct numer-ical simulation (DNS) can be memory intensive and computationally expensive.This is because calculating diffusion fluxes requires point-wise knowledge of themulticomponent diffusion coefficient matrix, which scales with the number ofchemical species squared [1]. The unity Lewis number, non-unity Lewis num-ber, and mixture-averaged diffusion assumptions have been used to reduce thecosts associated with mass diffusion by approximating the full diffusion coeffi-cient matrix as a constant scalar value, a constant vector, and a matrix diag-onal, respectively. In addition, several approaches further reduce the system’scomplexity by approximating multicomponent diffusion processes in terms ofequivalent Fickian processes, such as those used by Warnatz [2] and Coltrin etal. [3]. However, to our knowledge, the accuracy and appropriateness of theseassumptions have not been evaluated in turbulent reacting flows against mul-ticomponent diffusion transport due to its high computational expense and adearth of affordable computing tools.As further motivation for this study, Lapointe and Blanquart [4] recentlyinvestigated the impact of differential diffusion on simulations using unity andnonunity Lewis number approximations. They reported that methane, n -heptane,iso-octane, and toluene flames have similar normalized turbulent flame speedsand fuel burning rates when neglecting differential diffusion, but flames usingthe nonunity Lewis number approximation underpredict the normalized flamespeed when including differential diffusion due to reduced burning rates [4].Building on these results, Burali et al. [5] evaluated the relative accuracy of thenonunity Lewis number assumption relative to mixture-averaged diffusion forlean, unstable hydrogen/air flames; lean, turbulent n-heptane/air flames; andethylene/air coflow diffusion flames. They demonstrated that the relative errorassociated with the nonunity Lewis number assumption could be minimized withcareful selection of the Lewis number vector for a wide range of flames [5]. Simi-larly, Schlup and Blanquart [6] examined the impact of multicomponent thermaldiffusion on DNS of turbulent, premixed, high-Karlovitz hydrogen/air flames.They showed that simulations using the mixture-averaged thermal diffusion as-sumption underpredict the normalized flame speeds compared with results fromsimulations using full multicomponent thermal diffusion. In addition, includingmulticomponent thermal diffusion results in increased production of chemicalsource terms in regions of high positive curvature [6]. These observed discrep-ancies in similar flame simulations with different diffusion models warrant adetailed investigation of the fundamental transport phenomena involved.While data are sparse from three-dimensional reacting-flow simulations withmulticomponent transport, several groups have investigated the effects of mul-ticomponent transport in simpler configurations. These studies include one-dimensional [7–12] and two-dimensional flames [13–15] at various unburnt con-2itions. These works compared the multicomponent model with various lev-els of diffusion and transport property models, from constant Lewis numberto mixture-averaged properties. In general, prior studies found some errorsbetween multicomponent and mixture-averaged formulations for simplified hy-drogen/air and methane/air flame configurations, such as unstretched laminarflames. However, these studies did not assess flames where diffusion effectsmay be more important, such as two- and three-dimensional, unsteady laminarand turbulent flames. Moreover, advancing clean and efficient combustion tech-nology requires incorporating realistic fuel chemistry in large-scale turbulentsimulations relevant to practical applications. Thus, there is a clear need fora computationally efficient algorithm capable of modeling full multicomponentdiffusion transport [16].The studies by Lapointe and Blanquart [4], Burali et al. [5], and Schlup andBlanquart [6] each took care to isolate the diffusion assumptions in question byneglecting higher-order terms that may affect diffusion transport. For example,with the exception of Schlup and Blanquart [6], these studies neglected Soretand Dufour diffusion, as it would be difficult to determine the direct cause of anobserved effect when including both molecular and thermal diffusion. However,despite this methodical approach, the results of these studies were presentedwith reference to mixture-averaged diffusion, rather than full multicomponentdiffusion. This further highlights the need for a computationally efficient methodfor implementing full multicomponent transport, and a subsequent examinationof the differences between its “true” results and those resulting from the approx-imations conventionally used.In this direction, several studies have examined the impact of full multicom-ponent transport on simplified three-dimensional flame configurations. Gio-vangigli [14] demonstrated that multicomponent Soret effects significantly im-pact a wide range of laminar hydrogen/air flames. Specifically, they notedthat multicomponent Soret effects influence laminar flame speeds and extinc-tion stretch rates for flat and strained premixed flames, respectively. For high-pressure systems, Borchesi and Bellan [17] developed and analyzed a multi-species turbulent mixing model for large-eddy simulations. They focused onturbulent crossflow mixing of a five-species combustion-relevant mixture of n -heptane, O , CO , N , and H O. The multi-species transport model signifi-cantly improves the accuracy and fidelity of the solution throughout the mixinglayer; however, this study only considered non-reacting flows and, as a result, didnot assess the impact of multicomponent transport on the chemistry inherentin turbulent combustion. In addition, these simulations implement a simplifieddiffusion model to approximate multicomponent diffusion but do not directlysolve the diffusion terms present in the generalized conservation equations forspecies and energy [18].Motivated by the dearth of affordable three-dimensional multicomponenttransport models, Ern and Giovangigli [9, 19, 20] developed the computation-ally efficient Fortran library EGLIB for accurately determining transport coef-ficients in gas mixtures. More recently, Ambikasaran and Narayanaswamy [21]proposed an efficient algorithm to compute multicomponent diffusion velocities,3hich scales linearly with the number of species. This significantly reducescomputational cost compared with previous methods that directly invert theStephan–Maxwell equations and thus scale with the number of species cubed.Although both libraries reduce the computational cost of determining the mul-ticomponent diffusion coefficients, they do not provide a method for reducingthe resulting large memory requirements for multidimensional simulations.Overall, these prior studies provide compelling evidence that multicompo-nent transport is important and can affect the accuracy of combustion mod-els. However, none assessed how multicomponent transport impacts three-dimensional turbulent systems with detailed chemistry. In this article, wedemonstrate and analyze an efficient, dynamic algorithm that reduces the com-putational expense of calculating the multicomponent diffusion fluxes. Wedemonstrate the approach is accurate and stable for a wide range of time-stepsizes. In addition, we present a comprehensive assessment of the numericalcosts associated with this method. To verify the proposed algorithm we presentone-dimensional freely propagating, laminar hydrogen/air flames and comparewith results from Cantera. Finally, we simulate three-dimensional, turbulent,premixed, hydrogen/air flames. As a preliminary comparison of the mixture-averaged and multicomponent diffusion models, we perform an a posteriori as-sessment of how the choice of diffusion model impacts the turbulent statisticsof the three-dimensional hydrogen simulation.
2. Governing equations
This section presents the low-Mach number reacting Navier–Stokes equa-tions used in this study. In addition, this section outlines the method used todetermine the mass diffusion fluxes for both the mixture-averaged and multi-component approaches, abbreviated here as MA and MC, respectively.
In this work we solve the variable-density, low-Mach number, reacting-flowequations [22, 23]. The conservation equations are ∂ρ∂t + ∇ · ( ρ u ) = 0 , (1) ∂ρ u ∂t + ∇ · ( ρ u ⊗ u ) = −∇ p + ∇ · τ , (2) ∂ρT∂t + ∇ · ( ρ u T ) = ∇ · ( ρα ∇ T ) + ρ ˙ ω T − c p (cid:88) i c p,i j i · ∇ T + ραc p ∇ c p · ∇ T , (3) ∂ρY i ∂t + ∇ · ( ρ u Y i ) = −∇ · j i + ˙ ω i , (4)where ρ is the mixture density, u is the velocity vector, p is the hydrodynamicpressure, τ is the viscous stress tensor, T is the temperature, α is the mixturethermal diffusivity, c p,i is the constant-pressure specific heat of species i , c p is the constant-pressure specific heat of the mixture, j i is the diffusion flux of4pecies i , Y i is the mass fraction of species i , and ˙ ω i is the production rate ofspecies i . In Equation (3), the temperature source term ˙ ω T is given by ˙ ω T = − c − p (cid:88) i h i ( T ) ˙ ω i , (5)where h i ( T ) is the specific enthalpy of species i as a function of temperature.The density is determined from the ideal gas equation of state ρ = P o WRT , (6)where P o is the thermodynamic pressure, R is the universal gas constant, and W is the mixture molecular weight determined via W = (cid:16)(cid:80) Ni Y i /W i (cid:17) − , where W i is the molar mass of the i th species and N is the number of species.The diffusion fluxes are calculated with either the mixture-averaged [1] ormulticomponent [24] models, which are both based on Boltzmann’s equationfor the kinetic theory of gases [24, 25]. The baro-diffusion term is commonlyneglected in reacting-flow simulations under the low Mach-number approxima-tion [26]. We have also neglected thermal diffusion because our objective inthis work is to investigate the impact of mass diffusion models; Schlup andBlanquart [6, 27] previously explored the effects of thermal diffusion modeling. The i th species diffusion flux for the mixture-averaged diffusion model isrelated to the species gradients by a Fickian formulation and is expressed as j i = − ρD i,m Y i X i ∇ X i + ρY i u (cid:48) c , (7)where X i is the i th species mole fraction, D i,m is the i th species mixture-averaged diffusion coefficient as expressed by Bird et al. [1]: D i,m = 1 − Y i (cid:80) j (cid:54) = i X j / D ji , (8)where D ji is the binary diffusion coefficient between the i th and j th species.Finally, u (cid:48) c is the correction velocity used to ensure mass continuity: u (cid:48) c = N (cid:88) i =1 D i,m Y i X i ∇ X i . (9)The expression for species diffusion flux can be re-stated in terms of mass frac-tion Y i as j i = − ρD i,m (cid:32) ∇ Y i − Y i N (cid:88) k =1 ∇ Y k WW k (cid:33) + ρY i u (cid:48) c , (10)where D i,m corresponds to the i th element of the diagonal mixture-averageddiffusion coefficient matrix, defined herein as D MA .5 .3. Multicomponent (MC) species diffusion flux The multicomponent diffusion model for the i th species diffusion flux is j i = ρY i X i W N (cid:88) k =1 W k D ik ∇ X k , (11)where D ik is the ordinary multicomponent diffusion coefficient (computed usingthe MCMDIF subroutine of CHEMKIN II [28] with the method outlined by Dixon–Lewis [29]). Equation (11) can be restated in terms of mass fraction as j i = ρ (cid:88) k − D MC ik ∇ Y k , (12)where D MC ik = − W i W D ik − WW k (cid:88) j D ij Y j . (13)The diagonal of the ordinary multicomponent diffusion matrix, D ii , is zero. Aswill be shown later, the D MC matrix is singular with a kernel of dimension one.Interestingly, the vector of species mass fractions is in the kernel: N (cid:88) k =1 D MC ik Y k = − W i W (cid:32)(cid:88) k D ik Y k (cid:33) − (cid:32) W (cid:88) k Y k W k (cid:33) (cid:88) j D ij Y j = 0 . (14)This property will be important later (in Section 3.4) for the stability analysis.The multicomponent diffusion coefficients, thermal conductivities, and ther-mal diffusion coefficients are computed by solving a system of equations definedby the L matrix, composed of nine sub-matrices: L , L , L , L , L , L , L , a a a = XX , (15)where the right-hand side is composed of the one-dimensional mole fractionarrays X . Based on this system of equations, the inverse of the L , blockprovides the multicomponent diffusion coefficients: D ij = X i T P WW j ( q ij − q ii ) , (16)where q = (cid:0) L , (cid:1) − . (17)The L , sub-matrix block is given by L , ij = 16 T P N (cid:88) k =1 X k W i D ik { W j X j (1 − δ i,k ) − W i X i ( δ i,j − δ j,k ) } , (18)where δ i,j is the reduced dipole moment corresponding to the i th component ofthe vector of dipole moments. 6 . Methods As discussed previously, multicomponent mass diffusion has not yet beenincorporated into three-dimensional turbulent flame simulations due to its highcomputational expense. This section presents the discretized equations, nu-merical algorithm, and preconditioner proposed. The method is based on thesemi-implicit time-marching scheme for species mass-fraction fields proposed bySavard et al. [23].
This work was completed using the structured, multi-physics, and multi-scalefinite-difference code NGA [22, 23]. NGA can solve a wide range of problems,including laminar and turbulent flows [30–32], constant- and variable-densityflows [22, 33, 34], large-eddy simulation [31, 35], and DNS [33, 34, 36]. NGAdiscretely conserves mass, momentum, and kinetic energy with an arbitrarilyhigh-order spatial accuracy [22].NGA’s variable-density flow solver uses both spatially and temporally stag-gered variables, storing all scalar quantities ( ρ , P , T , Y i ) at the volume centersand velocity components at their respective volume faces [22, 37]. The con-vective term in the species transport equation is discretized using the bounded,quadratic, upwind-biased, interpolative convective scheme (BQUICK) [38]. Thediffusion source term is discretized using a second-order centered scheme andthe variables are advanced in time using a second-order semi-implicit Crank–Nicolson scheme [39].An iterative procedure is applied to fully cover the nonlinearities in theNavier–Stokes equations and the species diffusion terms. Prior studies demon-strated this iterative process to be critically important for stability and accu-racy [22, 23, 39, 40]. Savard et al. [23] fully detailed the numerical algorithmsequence; we summarize this method here. This summary is independent of thepreconditioning strategy employed in NGA, to which propose modifications inSection 3.2.The algorithm for advancing one time step follows, using a uniform time-step size ∆ t . The density, pressure, and scalar fields are advanced from timelevel t n +1 / to t n +3 / , and the velocity fields are advanced from time t n to t n +1 , where t n is the current time. Each iteration (i.e., time step) consists of Q sub-iterations and follows this procedure:0. Upon convergence of the previous time step, the algorithm stores the den-sity ( ρ n +1 / ), pressure ( P n +1 / ), velocity fields ( u n ), and scalar fields( Y n +1 / ), where Y represents the vector of species mass fractions ( Y , . . . , Y N ) .The solutions for pressure, species mass fraction, and momentum from theprevious time step are used as an initial guess for the iterative procedure: P n +3 / = P n +1 / , Y n +3 / = Y n +1 / , and ( ρ u ) n +10 = ( ρ u ) n . (19)7n Adams–Bashforth prediction evaluates the initial density: ρ n +3 / = 2 ρ n +1 / − ρ n − / , (20)which ensures that the continuity equation is discretely satisfied at thebeginning of the iterative procedure.1. For the sub-iterations k = 1 , . . . , Q , the semi-implicit Crank–Nicolsonmethod advances the scalar fields in time [39, 41]: ρ n +3 / k Y n +3 / k +1 = ρ n +1 / Y n +1 / + ∆ t ( C ∗ k + Diff ∗ k + Ω ∗ k )+ ∆ t (cid:18) ∂ C ∂ Y + ∂ Diff ∂ Y + ∂ Ω ∂ Y (cid:19) n +1 k · (cid:16) Y n +3 / k +1 − Y n +3 / k (cid:17) , (21)where Diff = −∇ · j i and Y ∗ k , C ∗ k , Diff ∗ k , and Ω ∗ k are the mass fraction,convection, diffusion, and chemical terms evaluated on the mid-point (orhalf time-step) scalar field Y ∗ k : Y ∗ k = Y n +1 / + Y n +3 / k . (22)To simplify the discrete notations for spatial differentiation, the operatorscorresponding to the convective and diffusive terms in Equation (4) arewritten as C and Diff , respectively [23]. ∂ C ∂ Y and ∂ Diff ∂ Y are the Jacobianmatrices corresponding to the convective and diffusive terms with respectto the species mass fractions, respectively. C and ∂ C ∂ Y are functions of thedensity and velocity, while Diff and ∂ Diff ∂ Y are functions of the density,diffusivity, and molar weight. They are consistently updated at each sub-iteration [23].2. The density field, ρ n +3 / k +1 , is evaluated from the new scalar fields usingEquation (6). We do not rescale the scalar fields as proposed by Shunn etal. [40]. However, upon convergence of the sub-iterations, this method isequivalent to the density treatment they proposed [23].3. The momentum equation is advanced in time using a similar semi-implicitCrank–Nicolson method for the scalar fields as described by Savard etal. [23].4. A Poisson equation is then solved for the fluctuating hydrodynamic pres-sure using a combination of HYPRE [22, 42], BICGSTAB [43], and/orFFTW [44]. The predicted velocity field is then updated.5. Upon convergence of the sub-iterations, the solutions are updated.The procedure summarized above becomes equivalent to the fully implicit Crank–Nicolson time-integration scheme upon convergence of the sub-iterations [39]. We expand the above numerical procedure to incorporate multicomponentdiffusion by modifying the time-marching step for species mass fraction fields.Specifically, this method modifies the treatment of the mass-diffusion sourceterm in the species mass fraction fields. All other intermediate steps are un-changed. 8 .2.1. Preconditioning iterative method
For simpler implementation, Equation (21) is solved in its residual form: (cid:20) ρ n +3 / k I − ∆ t (cid:18) ∂ C ∂ Y + ∂ Diff ∂ Y + ∂ Ω ∂ Y (cid:19) n +1 k (cid:21) · (cid:16) Y n +3 / k +1 − Y n +3 / k (cid:17) = ρ n +1 / Y n +1 / − ρ n +3 / k Y n +3 / k + ∆ t (cid:0) C n +1 k + Diff n +1 k + Ω ∗ k (cid:1) . (23)This equation can be restated as Y n +3 / k +1 = Y n +3 / k − ∆ t J − · Θ k , (24)where the matrix J is J = ρ n +3 / k I − ∆ t (cid:18) ∂ C ∂ Y + ∂ Diff ∂ Y + ∂ Ω ∂ Y (cid:19) n +1 k (25)and the vector Θ k = ρ n +3 / k Y n +3 / k − ρ n +1 / Y n +1 / ∆ t − (cid:2) C n +1 k + Diff n +1 k + Ω ∗ k (cid:3) (26)is the residual of the species transport equation at the previous sub-iteration,which asymptotes to zero as the sub-iterations fully converge.Written in its residual form, the time advancement of the species transportequations described here resembles the standard preconditioned Richardson-type iterative method [23, 45], where the matrix J acts as a preconditioner.The choice of J as a preconditioner is arbitrary and only affects the convergencecharacteristics of the iterative method [23]. For example, J = ρ n +3 / k I (27)is equivalent to the fully explicit integration of the convective, diffusive, andchemical source terms in the species transport equations. Alternatively, J = ρ n +3 / k I − ∆ t (cid:18) ∂ C ∂ Y + ∂ Diff ∂ Y + ∂ Ω ∂ Y (cid:19) n +1 k (28)is equivalent to fully implicit integration of the convective, diffusive, and chem-ical source terms [23].There is a clear tradeoff in selecting the preconditioner. Since precondi-tioning is applied to each step of the iterative methods, the form of matrix J should be optimized for low computational and inversion cost while maintain-ing strong convergence. The fully explicit preconditioner provides the cheapestoption but in our experience results in poor convergence performance, requir-ing extremely small time steps. Alternatively, the fully implicit preconditionerwould provide excellent convergence criteria and unconditional stability; how-ever, the Jacobian matrices for the chemical and diffusion source terms are9ypically dense [1, 28, 46]. Thus, constructing a fully implicit preconditioner isprohibitively expensive for large kinetic models.To achieve strong convergence while maintaining a low-cost form for thepreconditioner, we propose an approximation of the diffusion Jacobian thatlies between the fully implicit and fully explicit extremes: a semi-implicit pre-conditioner. Savard et al. [23] previously implemented a similar approach forpreconditioning the chemical Jacobian. In Equation (28), the Jacobian of the diffusion source term depends on themulticomponent diffusion flux, which is proportional to the multicomponent dif-fusion coefficient matrix, D MC . However, D MC is a dense matrix and wouldbe a computationally expensive approximation for the Jacobian. Alternatively,the mixture-averaged diffusion coefficient matrix, D MA , is a simplified approx-imation of D MC and thus may provide a reasonable, low-cost approximation ofthe fully implicit Jacobian.The mixture-averaged diffusion coefficient matrix, D MA , and the multicom-ponent diffusion coefficient matrix, D MC , are of a similar order and dependon the underlying species diffusivities. In addition, since D MA is computedfrom the local species and temperature values rather than global changes, itis inexpensive to compute. Finally, since D MA is strictly diagonal and thusinexpensive to invert, it provides a low-cost approximation to the diffusion Ja-cobian. In practice the approximate diffusion Jacobian is a tri-diagonal blockmatrix, where each block is the diagonal D MA matrix. In other words, for eachspecies the part of the Jacobian corresponding to that species is tri-diagonaland described by D MA . As mentioned previously, high-fidelity simulations with full multicomponentmass diffusion will have a high computational expense. Thus, to facilitate a cost-effective implementation of full multicomponent diffusion we propose a simpledynamic memory algorithm that significantly reduces the computational re-sources needed for such simulations.The cost of simulating full multicomponent diffusion comes from evaluatingthe D MC matrix. Thus, we can reduce computational cost significantly by lim-iting the evaluation of D MC to strictly once per grid-point. (In contrast, a naiveimplementation would involve repeated and redundant evaluations when calcu-lating the species diffusion flux vector and its gradient.) This is possible becausethe central-difference scheme used is linear and thus additive and commutativeby nature. In other words, the terms in the discretized equation are simplyadded together, and thus are strictly independent of each other and require noinformation from the surrounding grid points.Recognizing this, it follows that the order of addition does not matter solong as all of the appropriate terms are included in the discretization. Thus, wecan calculate the D MC matrix once per grid point, and calculate and store for10ach species the discrete terms of the discretized scalar field corresponding onlyto the information available at that grid point. The process then repeats at thenext grid point and fills in the remaining information. This approach is simplya memory-efficient rearrangement of the floating-point operations and does notalter the final result. Moreover, this dynamic memory scheme avoids the needto calculate local gradients at each grid point.In practice, we calculate and store the portions of the enthalpy and species-diffusion source terms (in Equations (3) and (4), respectively) that can be com-puted from the information available at the i th grid-point for the ( i − / and ( i + 1 / flux vectors. For example, the discretized form of the diffusion sourceterm is Diff i = −∇ · j i = − j i +1 / + j i − / ∆ x = (cid:20) ( ρ i D i + ρ i +1 D i +1 ) Y i +1 − Y i ∆ x − ( ρ i − D i − + ρ i D i ) Y i − Y i − ∆ x (cid:21) x , (29)where the diffusion source term contributions from the i − , i , and i + 1 gridpoints are Source i − = ρ i − D i − x ( Y i − Y − i ) , (30)Source i = ρ i D i x ( Y i +1 + Y i − − Y i ) , and (31)Source i +1 = ρ i +1 D i +1 x ( Y i +1 − Y i ) , (32)respectively.At the i th grid point, information on the diffusion coefficients at the i − and i + 1 grid points is not available; thus, only the diffusion coefficients for the i thgrid point can be stored. However, by recognizing that D i at the i th grid pointis equal to D i +1 and D i − at the i − and i + 1 grid points, respectively, it ispossible to solve Equation (30), Equation (31), and Equation (32) for the i + 1 , i , and i − grid points, and store them in their respective memory locations. Atthe next grid point ( i + 1 ) the process repeats and the remaining informationfor the i th grid point is calculated and added to the previously stored partialsolution, thus completing the information needed at the i th grid point. Figure 1summarizes this process; fluxes are located at cell faces while source terms areat cell centers.This approach reduces the number of D MC evaluations from once per speciesper grid point to strictly once per grid point. Finally, it reduces temporarymemory requirements from an array sized n x × n y × n z × N to a × arraycorresponding to only the information needed at the current grid point ( i, j, k )and its six surrounding points, where n x , n y , and n z are the numbers of gridpoints in the x , y , and z directions. This optimizes performance by reducingcache calls for both the species mass fractions and species diffusion coefficients.11 i − i + 1 Diff in Diff out for i=1:X do Calculate diffusion coefficient matrix; for isc=1:N do Flux ( i − /
2) +=
Diff in ,isc ;Flux ( i + 1 /
2) +=
Diff out ,isc ; end Source ( i ) += influence from Diff ( i − / and Diff ( i + 1 / ;Source ( i −
1) += influence from Diff ( i − / ;Source ( i + 1) += influence from Diff ( i + 1 / ; end Figure 1: Dynamic algorithm for calculating multicomponent enthalpy and species diffusionsource terms. Fluxes are located at cell faces while source terms are at cell centers. N is thenumber of species. The algorithm is most efficient for a structured grid, but the proposedmethod is easily extendable to finite-volume discretizations on unstructuredmeshes with scalars located at the cell centers. In such schemes, the diffusionterm is written as the sum of fluxes on each cell surface. In turn, these fluxesare written as differences of cell-averaged scalar values. The regrouping of thecontributions of the diffusion term to each cell in Equations (30)–(32) wouldfollow a similar approach.
To evaluate the theoretical stability of the proposed treatment of the dif-fusion source terms, we will perform a one-dimensional von Neumann stabilityanalysis. First, we decompose the vector of species mass fractions into the exactsteady-state solution ( Y ◦ ) and a small perturbation vector. Then, we expandthis perturbation in a Fourier series by assuming a solution of the form Y ( x, t ) = Y ◦ ( x ) + f ( t ) e iκx , (33)where κ is the wavenumber and f ( t ) is the time-varying amplitude of the per-turbation. Under small deviations from a steady-state solution, we can makethe simplifying assumption that ρ n +3 / k ≈ ρ n +1 / = ρ ◦ . (34)Similarly, all diffusion coefficients are evaluated from the steady-state solution.12rom here, we rewrite Equation (21) in a point-wise form neglecting boththe chemical source term—demonstrated to be stable by Savard et al. [23]—andthe convective transport term, which is integrated explicitly in this stabilityanalysis (i.e., not modified by sub-iterations). This transforms the set of N partial differential equations into a set of N ordinary differential equations,where N is the number of species. Equation (23) reduces to the form (cid:18) I + ∆ t D MA κ (cid:48) (cid:19) (cid:16) f n +3 / k +1 − f n +3 / k (cid:17) = f n +1 / − f n +3 / k − ∆ t D MC κ (cid:48) (cid:16) f n +3 / k + f n +1 / (cid:17) , (35)where κ (cid:48) is the modified wavenumber, and D MA and D MC are the mixture-averaged and multicomponent diffusion coefficient matrices calculated from Equa-tions (8) and (13), respectively. For the second-order central differencing schemeused, κ (cid:48) takes the form κ (cid:48) = 2∆ x [1 − cos( κ ∆ x )] . (36)While here we apply this to a second-order central difference scheme, the stabil-ity analysis holds for any spatial discretization of the diffusion terms in Equa-tion (23). In the present case, the most unstable mode manifests as cell-to-celloscillations corresponding to κ = π/ ∆ x and κ (cid:48) = 4 / ∆ x .Recall that f n +1 / is the value at the previous time step as defined in step 0of Section 3.1 and f n +3 / ≡ f n +1 / . (37)Dropping the superscripts for clarity, we can reduce Equation (35) to f k +1 = Af + Bf k , (38)where A = (cid:18) I + ∆ t κ (cid:48) D MA (cid:19) − (cid:18) I − ∆ t κ (cid:48) D MC (cid:19) (39)and B = (cid:18) I + ∆ t κ (cid:48) D MA (cid:19) − (cid:18) ∆ t κ (cid:48) (cid:0) D MA − D MC (cid:1)(cid:19) . (40)Inspecting Equation (38), matrix A is multiplied by the constant value of theprevious time step ( f ) and therefore does not contribute to the stability of thesub-iterations. We focus on the properties of the B matrix, which acts as theamplification/growth factor. Theoretically, the stability of the sub-iterations isensured if the spectral radius of matrix B , defined as the largest absolute valueof the eigenvalues, is less than one: ρ ( B ) ≤ . (41)The matrix B has some interesting properties that deserve further discussion.13 MC† D MA† B imp B exp Table 1: Eigenvalues for the multi-component (left) and mixture-averaged (center-left) diffu-sion matrices ( D MC and D MA ) and absolute values of the eigenvalues for the amplificationmatrix ( B ) for the implicit formulation (center-right) and explicit formulation (right) evalu-ated on the burned side of the lean hydrogen premixed flame (see section 4.1). † units are10 −3 m /s. A time-step size of ∆ t = −5 s was used for B . First, this matrix is proportional to the difference between the two diffusionmatrices D MA and D MC . Recall that the D MA is a purely diagonal matrix.Table 1 compares the eigenvalues of these matrices on the burned side of thelean hydrogen premixed flames (see Section 4.1 for details on the flame). Theburned side of the flame is characterized by the largest diffusion coefficients andis expected to be the most unstable location within a flame as far as diffusion isconcerned. The two sets of eigenvalues are extremely close, which is expected asthe mixture-averaged diffusion model approximates the multi-component diffu-sion model. As a result, the norm of the difference of the two matrices isexpected to be much less than the norm of either matrix. In other words, weanticipate that ρ ( B ) (cid:28) regardless of the time-step size ( ∆ t ) and grid spacing( ∆ x ).One noticeable difference between D MC and D MA is the presence of a nulleigenvalue for multi-component diffusion. As described in Section 2.3, the multi-component matrix is singular, and its kernel is spanned by the species massfraction matrix, here Y ◦ . Consider the special case of f = Y ◦ . Leveraging thefact that f lies in the kernel of D MC , the amplitude at the next sub-iterationwill be f = ( A + B ) f = f . (42)By recursive reasoning, one can show the property holds for all sub-iterations.In other words, this mode is unaffected by the iterative process and remainsthe same between time steps. This time-invariant mode is nothing more thanthe steady-state solution Y ◦ . To avoid “double counting” in Equation (33), theeigenvalue analysis of the matrix B should be performed on the linear space notincluding the vector Y ◦ .Table 1 provides an example of the eigenvalues of the amplification matrix B for a time-step size of ∆ t = 10 − s . Practically, the eigenvector associatedwith the largest eigenvalue of B forms a small angle ( ∼ . deg) with the14pecies mass fraction vector ( Y ◦ ). Hence, most of it is in the kernel of D MC .Following the previous discussion, this eigenvalue (indicated in italics) shouldnot be considered, and the overall stability is controlled by the second-largesteigenvalue, in this case . . As expected, this eigenvalue is much less thanunity, thus proving the stability of the iterative procedure. This should becompared to the stability of the explicit formulation obtained by setting D MA =0 in Eq. (40). Under these conditions, the spectral radius of B becomes ρ ( B ) = 2 ∆ t ∆ x ρ (cid:0) D MC (cid:1) ≈ t ∆ x max (cid:0) D MA (cid:1) , (43)which resembles a Fourier number. As shown by the large eigenvalue of theexplicit-method amplification matrix, B exp , in Table 1, solving the system ofequations would not be stable at ∆ t = 10 − s without the proposed implicitformulation. Section 5.1 presents an in-depth comparison of this theoreticalstability criterion against practical numerical convergence results for a one-dimensional freely propagating flame.
4. Test cases
We will evaluate the performance of the proposed iterative method and therelative cost of the implemented memory algorithm in Section 5. We base ourevaluation on two flow configurations: a one-dimensional, unstretched, laminarflame and a three-dimensional, statistically stationary, turbulent flame; bothare premixed hydrogen/air flames. All simulations used the same nine-specieshydrogen mechanism of Hong et al. [47] with updated rate constants from thesame group [48, 49]. This section describes the configuration and conditions usedfor the one- and three-dimensional simulations used for this study. Appendix Bincludes additional method verification.
To verify the implementation of the multicomponent mass-diffusion modeland evaluate its accuracy, we performed one-dimensional, unstretched (flat),laminar flame simulations and compared these with similar mixture-averagedand multicomponent results computed using Cantera [50]. We selected the one-dimensional flat flame configuration because it restricts all transport to thestreamwise direction. As a result, the spanwise fluxes are zero by definition forthis geometry. This condition may not hold in a multidimensional flow simu-lation where the multicomponent diffusion fluxes may be misaligned with thespecies gradient vector. This simplified geometry allows us to directly com-pare the multicomponent mass diffusion model to the commonly used mixture-averaged diffusion model.The simulations used an unburnt temperature of 298 K and pressure of 1 atm,with an equivalence ratio of φ = 0 . and inlet velocity equal to the laminar flamespeed for all Cantera and NGA cases. The flame was centered in a computa-tional domain comprised of 720 grid points where ∆ x = l F = ( T max − T min ) / |∇ T | max . Schlup and Blanquart [6] used an identi-cal configuration to investigate the impact of Soret and Dufour thermal diffusioneffects.We ran the Cantera simulations similarly using both mixture-averaged andmulticomponent diffusion models with matching inlet conditions, equivalence ra-tio, and domain size. The freely-propagating adiabatic flat flame solver ( FreeFlame )was used with grid refinement criteria for both slope and curvature set to 0.1and a refinement ratio of 2.0 for 860 grid-points.
We simulated a three-dimensional, turbulent, premixed, freely propagatingflame as a test of the proposed algorithm for multicomponent mass diffusionand to assess the impact of diffusion model choice on global statistics such asthe turbulent flame speed. The computational domain consists of inflow andconvective outflow boundary conditions in the streamwise direction. The twospanwise directions use periodic boundaries. We set the inflow velocity to themean turbulent flame speed, which keeps the flame statistically stationary suchthat turbulent statistics can be collected over an arbitrarily long run time. Inthe absence of mean shear, we use a linear turbulence-forcing method [33, 51]to maintain the production of turbulent kinetic energy through the flame. Wecarefully selected the Karlovitz number to fall within the distributed reactionzone regime while avoiding the broken reaction zone regime [52]. Moreover,the computational setup for this case is similar to those of Lapointe et al. [52],Burali et al. [5], and Schlup et at. [6], who studied differential-diffusion effects,local extinction, and flame broadening using the mixture-averaged model andconstant non-unity Lewis number assumptions.Table 2 provides further details of the computational domain, unburnt mix-ture, and inlet turbulence. The unburnt temperature and pressure are 298 Kand 1 atm, respectively. The inlet equivalence ratio is φ = 0 . , with an unburntKarlovitz number Ka u = τ F /τ η = 149 , where τ F = l F /S L is the flame time scaleand τ η = ( ν u /(cid:15) ) / is the Kolmogorov time scale of the incoming turbulence withunburnt kinematic viscosity ν u and turbulent energy dissipation (cid:15) . The unburntturbulent Reynolds number is Re t = u (cid:48) l/ν u = 289 , where u (cid:48) is the fluctuation ofthe mean velocity and l is the integral length scale. The mean inflow velocity atthe inlet boundary condition approximately matches the turbulent flame speedso that the flame remains relatively centered in the domain and we can performarbitrarily long simulations. Once the turbulence has fully developed, we runthe simulations for 22 eddy turnover times, τ eddy = k/(cid:15) ≈ µ s .The domain has 1520 points in the streamwise direction and 190 points inboth spanwise directions, with a uniform grid size of ∆ x = l F / . This domainis about l F long and l F in the spanwise directions. Given the prescribedturbulence intensity, this mesh has a grid spacing equivalent to ∆ x ≈ η u , where η u is the Kolmogorov length scale for the unburnt region; this resolution im-proves in the burnt region of the flame. Lapointe et al. [52] previously confirmed16 able 2: Three-dimensional simulations parameters. ∆ x is the grid spacing, η u is Kolmogorovlength scale of the unburnt gas, ∆ t is the simulation time-step size, φ is the equivalence ratio, P is the thermodynamic pressure, T u is the temperature of the unburnt mixture, T peak is thetemperature of peak fuel consumption rate in the one-dimensional laminar flame, S L is thelaminar flame speed, l F = ( T b − T u ) / |∇ T | max is the laminar flame thickness, l = u (cid:48) /(cid:15) is theintegral length scale, u (cid:48) is the turbulence fluctuations, (cid:15) is the turbulent energy dissipationrate, Ka u is the Karlovitz number of the unburnt mixture, Re t is the turbulent Reynoldsnumber of the unburnt mixture, and ν u is the unburnt kinematic viscosity. Parameter MA MCDomain L × L × LL x Grid × × x [m] 4.24 × −5 η u [m] 2.1 × −5 ∆ t [s] 6 × −7 φ P [atm] 1 T u [K] 298 T peak [K] 1190 1180 S L [m/s] 0.230 0.223 l F [mm] 0.643 0.631 l/l F u (cid:48) /S L
18 18.6Ka u = τ F /τ η
149 151Re t = ( u (cid:48) l ) /ν u T peak defining the flamefront, where T peak is the temperature of peak fuel consumption rate in the one-dimensional laminar flame. The flame surface shows the complex behavior ofthe flame in the turbulent field. End forcing y L t Begin I n Flow �°d� ��( �sc \) Out Flow �0('�
0 0.25 X L forcing
Figure 2: Two-dimensional schematic of the three-dimensional flame configuration. Adaptedfrom Burali et al. and Schlup and Blanquart [5, 6]. The red line indicates the approximatelocation of the flame.Figure 3: Iso-surface of peak temperature colored by OH mass fraction for a three-dimensionalturbulent hydrogen/air flame with multicomponent mass diffusion.
5. Results and discussion
To start, we present a practical assessment of the method’s convergence andstability, by comparing the numerical rate of convergence to the theoretical rateof convergence. Following this demonstration of the proposed method’s stabil-ity, we verify the accuracy of the method through a posteriori assessment of18ne-dimensional, unstretched, premixed, laminar flame simulations. Finally, wepresent a preliminary evaluation of the relative differences between the mixture-averaged and multicomponent diffusion models for the three-dimensional tur-bulent premixed flame simulations.
We use the one-dimensional flame to numerically evaluate the convergencestability of the sub-iterations with respect to time-step size. The simulationsfor these tests were initialized from a mixture-averaged data file to provide aworst-case scenario for the initial iterative step in converging to the multicom-ponent solution. While the theoretical analysis was performed assuming explicittransport of the convective terms and constant density/diffusion coefficients, weperformed this test with semi-implicit transport and variable density/diffusioncoefficients. This demonstrates the stability of the proposed preconditioner forthe semi-implicit multicomponent diffusion transport in a practical numericalsimulation. Savard et al. [23] previously showed the numerical stability of thechemical and convective terms, so we do not discuss these terms in detail in thisanalysis. -10 ImplicitExplicit (a) ∆ t = 10 − s -10 k ImplicitExplicitFirst mode (b) ∆ t = 10 − s Figure 4: Convergence of the density residual as a function of sub-iteration for the proposedsemi-implicit method, for a smaller and larger time-step size. Dashed lines are the spectralradii shown in Figure 5 and are determined by numerically fitting an exponential curve to theslope of the density residual.
We focus on the maximum density residual over the whole domain, becauseits convergence is controlled by the convergence of all chemical species. Fig-ures 4a and 4b present the density residuals as a function of sub-iteration,starting from the initial time step, for a small and large time-step size, respec-tively. For the time-step sizes tested, converging (as opposed to converged)sub-iterations implies a stable simulation, which agrees with behavior shown bySavard et al. [23]. In other words, unless the sub-iterations diverge, the simula-tion remains stable. As expected, the explicit method diverges quickly even at19ery small time-step sizes (Figure 4a), while the semi-implicit method remainsstable up to a time-step size of ∆ t ≤ × −5 s (Figure 4b).The rate of convergence of the sub-iterations for each of the source termsin Figures 4a and 4b follows an exponential relationship, i.e., Res k ∼ r k , whereRes k is the residual of the k th sub-iteration and r is the convergence rate. Wecompute the numerical convergence rate r by fitting an exponential curve tothe slope of the density residuals; the convergence rate is represented by dashedlines in Figure 4. Since density is a function of the species mass fractions,its convergence rate should tend towards that of the slowest-converging speciesmass fraction. -8 -6 -4 -2 implicit theoryexplicit theoryimplicit numericalexplicit numerical Figure 5: Theoretical convergence rate determined from diagonalizing matrix B correspondingto the worst-case modified wavenumber for the one-dimensional premixed flame, comparedwith the numerical convergence rates determined by fitting an exponential curve to the slopeof the density residual. Figure 5 compares the results of the theoretical and numerical stability anal-yses, showing the spectral radius of matrix B as a function of the time-step sizefor the one-dimensional test case. For the explicit scheme, the theoretical andnumerical results agree well for the full range of time-step sizes. However, for theimplicit scheme, the predicted spectral radius is much smaller than the measuredone. The proposed implicit formulation thus yields a convergence rate that isnot limited by diffusion, but rather constrained by other processes that werenot considered in the stability analysis, such as chemistry. The predicted con-vergence rate can nonetheless be observed in the implicit case. As mentionedin Section 3.4, the eigenvector associated with the largest eigenvalue of B imp inTable 1 forms a small angle with Y ◦ . Hence, a fraction of the error, albeit tiny,is associated with this eigenvector, which will slowly converge. In Figure 4b, theconvergence rate for the last sub-iterations of the implicit case closely matches20hat eigenvalue.Overall, these results suggest that the theory well-approximates actual sta-bility and provides a practical limit for the numerical stability of the proposedalgorithm. To verify the multicomponent model, we compare a posteriori the one-dimensional unstretched species profiles and laminar flame speeds. Figure 6compares the nine species profiles for the steady-state one-dimensional flat flamesolutions relative to local mixture temperature for the multicomponent andmixture-averaged models from both NGA and Cantera. The profiles all agreewithin 1 % at all points, with the exception of N . The laminar flame speeds( S oL ) for these simulations are approximately 23.0 cm/s and 22.3 cm/s for themixture-averaged and multicomponent diffusion NGA cases, respectively; thelaminar flame speeds for both cases agree with those from Cantera within 1 %.The unstretched laminar flame speed is S oL = − (cid:82) ρ ˙ ω H dxρ u Y H ,u , (44)where ρ u is the unburnt mixture density and Y H ,u is the unburnt fuel massfraction. We attribute the larger difference in the species profile for N to thecorrection velocity term associated with the mixture-averaged diffusion model,which is weighted by mass fraction and thus can be heavily impacted by dif-ferences in N due to its high concentration throughout the flame. The minordifferences between the multicomponent species profiles are less than 1 % at allpoints. The strong agreement between the other eight species profiles for boththe NGA and Cantera results verifies the multicomponent model’s functionality. With the proposed algorithm’s stability limits and functionality verified, wenow examine the accuracy for a given stable simulation. We determine theorder of accuracy of the method based on the 1D freely propagating flame caseby examining the power-law dependence of the error as a function of the time-step size.Figure 7 shows the normalized error for the 1D freely propagating flame casefor various time steps. We initialize the simulation using an input flame profilecorresponding to a fully converged statistically stationary flame, generated witha time-step size of ∆ t = 1 × − s and seven sub-iterations. A wall is thenset at the simulation inlet, allowing the flame to propagate upstream in thedomain. We then let the reference flame propagate for two flame pass-throughtimes to ensure a fully converged freely propagating flame profile free of anyinitial transients due to the transition from the input stationary flame profile.This reference file then serves as the input for a set of freely propagating flameswith time-step sizes ranging 10 −5 –10 −7 s and for seven sub-iterations. Finally,we allow these test flames to propagate for an additional flame pass-through21
00 800 12000.7560.7580.76 400 800 1200024 10 -5
400 800 12000.150.20.25400 800 12000510 10 -4
400 800 12000510 10 -4
400 800 120000.0050.01400 800 120000.050.1 400 800 1200012 10 -4
400 800 12000510 10 -5 NGA MC NGA MA Cantera MC Cantera MA
Figure 6: A posteriori comparisons of species mass fractions relative to mixture local temper-ature in a hydrogen/air flame with φ = 0 . using NGA and Cantera. time to ensure statistical independence from the initial reference-flame inputfile.With the freely propagating flame tests completed, we interpolated thespecies and density fields to a constant temperature space corresponding tothe temperature distribution in the flame region. This interpolation ensures adirect comparison of the species error independent of variation in the temper-ature space over the range of time-step sizes. We then calculate error as the L -norm of the species and density profiles in temperature space, relative to thereference flame profile with ∆ t = 1 × − s and seven sub-iterations:error = (cid:115) (cid:82) ( Y i − Y i, ref ) dT (cid:82) Y i, ref dT (45)22 -8 -6 -4 -2 y = -x H H OOHH (a) Seven sub-iterations -6 -4 -2 y = -x (b) Four sub-iterations Figure 7: Relative accuracy of the method as a function of time step size for the one-dimensional, freely propagating flame test case with seven sub-iterations. Errors are definedas the absolute difference of their integrated value in temperature space compared with areference solution obtained for ∆ t = 1 × − s and seven sub-iterations. Black dashed linecorresponds to y = x − . and error = (cid:115) (cid:82) ( ρ − ρ ref ) dT (cid:82) ρ ref dT . (46)We selected the species H , H O, OH, and H to evaluate the accuracy of themethod because they represent the reactants, intermediate species, and prod-ucts present in hydrogen combustion. Density ( ρ ) is also included to globallyassess error, since it depends on all species. As shown in Figure 7a, all quan-tities exhibit second-order accuracy in time with seven sub-iterations. The er-rors corresponding to the L - and L ∞ -norms are similar in magnitude and alsodemonstrate second-order accuracy in time with seven sub-iterations.While the method is fully second-order accurate for seven sub-iterationsand above, the solution transitions to first-order accuracy as the number ofsub-iterations decreases. Figure 7b shows that the solution exhibits first-orderaccuracy when using four sub-iterations. Between four and seven sub-iterationsthe solution is second-order accurate for large time-step sizes but transitionsto first-order accuracy as the time step size decreases. The range of time-stepsizes that achieve second-order accuracy grows until the solution becomes fullysecond-order accurate at seven sub-iterations for all time-step sizes considered.To evaluate the of absolute magnitude of error associated with the proposedmethod, as opposed to the order of accuracy (as time step size approaches zero),Figures 8a and 8b present the temperature as a function of distance and fuelmass fraction as a function of temperature, respectively, for a range of freelypropagating flames with several time-step sizes and sub-iterations. The solutions23xhibit negligible error in both temperature and fuel mass fractions for the time-step sizes considered, and even when using as few as four sub-iterations; thesetests demonstrate the high accuracy and robustness of the proposed method. t=1e-7, 7 sub t=2e-7, 4 sub t=5e-7, 4 sub t=1e-6, 4 sub t=2e-6, 4 sub (a) Temperature vs. distance
500 1000 150000.0020.0040.0060.0080.010.012 t=1e-7, 7 sub t=2e-7, 4 sub t=5e-7, 4 sub t=1e-6, 4 sub t=2e-6, 4 sub (b) H mass fraction vs. temperature Figure 8: Impact of time-step size and number of sub-iterations on the accuracy of one-dimensional freely propagating flames.
In this section we assess a posteriori the species mass diffusion fluxes inthe doubly periodic three-dimensional flames [4–6]. Differential diffusion effectscause the instabilities found in lean hydrogen/air flames, and at high Karlovitznumbers the turbulence time scales match the order of diffusion time scales.To assess the impact of the mixture-averaged and multicomponent massdiffusion models on flame chemistry, we compare a posteriori the turbulent andchemistry statistics. We allow the flames to develop in a turbulent flow field,and compute the statistics after the transients from the initial flow and scalarfields have advected through the domain. As an initial assessment, we calculatethe effective turbulent flame propagation speeds: S T = − (cid:82) V ρ ˙ ω H dVρ u Y H ,u L . (47)Figure 9 shows the time history of the turbulent flame speed over twenty-two eddy turn-over times ( τ eddy ). The average normalized flames speeds fromthe mixture-averaged and multicomponent models differ by 15 %: S MA T /S L =29 . and S MC T /S L = 34 . , respectively. Further study is needed on whetherthe mixture-averaged diffusion model fully captures the fundamental physics ofmulticomponent diffusion.To further assess any differences between the mixture-averaged and multi-component mass diffusion models, Figure 10 presents the means of fuel mass24 MCMA
Figure 9: Turbulent flame speed history for three-dimensional, freely propagating, premixed,turbulent hydrogen/air flame with φ = 0 . . fraction and its source term conditioned on temperature for the full time do-main. The differences in the calculated conditional means are small: less than5.5 %. This agreement also extends into super-adiabatic regions for the hy-drogen/air flame; these regions, also called “hot spots”, result from differentialdiffusion and have been predicted both in theoretical studies [53] and numericalanalyses of lean hydrogen/air mixtures [54–56]. However, these small differ-ences in global flame statistics do not explain the 15 % difference observed inthe turbulent flame speeds between the mixture-averaged and multicomponentdiffusion models. These results raise questions on the appropriateness of themixture-averaged diffusion assumption for direct numerical simulation and war-rants further investigation. MCMA (a) Fuel mass fraction (b) Fuel chemical source term
Figure 10: Conditional means on temperature for the three-dimensional, freely propagating,premixed, turbulent hydrogen/air flame with φ = 0 . . .5. Computational cost This section discusses the relative cost for implementing the full multicom-ponent mass diffusion to provide context for its use. The presented timing com-parisons examine how the method scales with both number of chemical speciesand spatial dimension.We tested three chemical kinetic models (containing 9 [5], 35 [36, 37], and 172species [57–59]) in a one-dimensional flat flame simulation to determine the costof multicomponent mass diffusion over a wide range of model sizes. Figure 11shows the computational time per grid point for computing the diffusion massfluxes on a desktop workstation using an Intel Xeon-X5660 CPU with a 2.80 GHzclock speed. The presented timings include calculation of both the diffusioncoefficients and mass diffusion fluxes for all aspects of the code. y = xy = x MCMA
Figure 11: Computational time per grid point for computing diffusion coefficients and diffusionmass fluxes using kinetic models with 9, 35, and 172 species; black dashed lines correspondto linear ( y = x ) and quadratic ( y = x ) scaling trends respectively. MC and MA stand formulticomponent and mixture-averaged, respectively. For the tested chemical kinetic models, the mixture-averaged model scaleslinearly while the multicomponent model scales quadratically with the numberof species. The multicomponent model is more expensive and does take moretime per-point for all three test cases. For the largest kinetic model (with 172species) the multicomponent case is noticeably more expensive than the mixture-averaged model. The increased cost for the multicomponent simulations comesprimarily from the CHEMKIN II [28] routine used to determine the ordinarymulticomponent diffusion coefficient matrix.The relevant cost for the proposed method can be split into three primarycategories: the costs of calculating the multicomponent diffusion coefficients,26alculating the multicomponent diffusion fluxes, and the semi-implicit integra-tion scheme. Since the proposed method for implementing full multicomponentmass diffusion focuses on efficient low-memory calculation of the diffusion fluxes,rather than the multicomponent diffusion coefficients, the cost of CHEMKINshould be considered independently of the proposed algorithm. Moreover, thesemi-implicit scheme is the same for the mixture-averaged and multicomponentcases, because both cases use the mixture-averaged diffusion coefficient matrixto approximate the Jacobian for the diffusion source terms. As a result, thetwo methods have similar implementation and computational expense, with theexception of using CHEMKIN II [28]. S ca l a r C h e m i s t r y D i ff u s i on C h e m k i n V e l o c it y P r e ss u r e R e s t (a) Mixture-averaged model S ca l a r C h e m i s t r y D i ff u s i on C h e m k i n V e l o c it y P r e ss u r e R e s t (b) Multicomponent model Figure 12: Computational time per grid point for each of the three flame configurations: onedimensional (blue), two dimensional (red), and three dimensional (yellow).
ScalarChemistryDiffusionVelocityPressureRestChemkin
Figure 13: Comparison of numerical costs for the three-dimensional hydrogen flame simula-tions.
To evaluate how the multicomponent model scales with increasing spatialdimension, and evaluate the relative cost of using CHEMKIN II [28], we ac-quired timings for one- ( grid points), two- ( × grid points), and27hree-dimensional ( × × grid points) configurations covering thecases presented in this work, with the additional two-dimensional case match-ing similar timing tests by Schlup et al. [6]. These timing tests represent anaverage cost per point and are determined by averaging the timings taken forthe 20 time steps, skipping the first and last integrations. Figure 12 presentsthe computational timings for each part of the code for both diffusion models,where “Scalar” includes scalar field calculation; “Diffusion” includes the flux cal-culation and D MA calculation for the implicit solver; “Chemistry”, “Velocity”,and “Pressure” are as named; and “Rest” account for any remaining compu-tations. “Scalar” includes the semi-implicit solver for integrating the diffusionsource terms, while the semi-implicit solvers for chemistry and velocity are in-cluded in their named categories. To facilitate comparison between the twomodels, Figure 13 presents the total computational time per grid point for boththree-dimensional hydrogen simulations as a stacked bar chart broken down byeach section of code. We performed these computations on the National EnergyResearch Scientific Computing Center (NERSC) high-performance computingcluster Cori (Cray XC40) [60].While much of the code exhibits a similar cost per grid point, regardlessof the dimensionality of the problem, the chemistry is more expensive for theone- and two-dimensional cases. This cost increase is due to NGA’s structure.NGA was written and optimized for three-dimensional configurations, thus theone- and two-dimensional cases are artificially more expensive, especially inthe chemistry calculations [6]. In addition, for three dimensions the cost ofthe pressure solver increases due to using the HYPRE package [42]. The one-and two-dimensional cases both implement an exact FFT-tridiagonal solver,while HYPRE—used for the three-dimensional cases—is iterative and thus moreexpensive. Despite the minor increase in cost for the pressure solver in threedimensions, the cost is negligible when considering larger kinetic models (i.e.,more than 35 species).Consistent with Figure 11, the cost of calculating “Diffusion” increases withmodel complexity; recall that D MA is calculated for both the mixture-averagedand multicomponent solvers. However, the multicomponent diffusion mass fluxcalculation represents only 21 % of the total simulation time for the three-dimensional case. As expected, the cost of calling CHEMKIN II for the diffusioncoefficients is large and accounts for roughly 23 % of the three-dimensional sim-ulation time. Interestingly, the cost of diffusion increases only slightly movingfrom one dimension to two dimensions. This results from the high efficiency ofthe dynamic memory-allocation algorithm used to implement this model (seeSection 3.3). Moreover, the multicomponent diffusion implementation is lessexpensive than the mixture-averaged model for the one-dimensional case andequivalent in cost for the two-dimensional case. Overall, by reducing memoryrequirements and optimizing calls to memory, the memory algorithm imple-mented for the multicomponent model maintains low computational expense.These results indicate that, for hydrogen-air combustion, the multicompo-nent model is more expensive than the mixture-averaged model; however, thedifferences in “Diffusion” costs between the two models are due to the use of28HEMKIN II [28]. Thus, the slowdown could be minimized by implementinga more-efficient package for calculating the mass-diffusion coefficients such asEGLIB [9, 19, 20]; however, the total cost of computing mass diffusion fluxesremains notable, even for the mixture-averaged case.
6. Summary and future work
This article presents an efficient and stable scheme for implementing multi-component mass diffusion in reacting-flow DNS with minimal memory expense.The proposed scheme exhibits reasonable computational cost for chemical ki-netic models of up to 100 species; this performance could be further improvedby implementing a more-efficient method for calculating the multicomponentdiffusion coefficient matrix.The results presented for hydrogen flames suggest that the mixture-averagedmass diffusion model may suffice for DNS of three-dimensional, premixed turbu-lent flames in the regimes and configurations considered. However, we observeda 15 % difference in the turbulent flame speeds between the two models, thoughthe differences in the conditional means of the fuel source term and mass frac-tion were negligible. The difference observed in turbulent flame speeds raisesquestions about using the mixture-averaged model in DNS of turbulent reactingflows. Moreover, the algorithm proposed in this study provides a fast, efficient,method for implementing multicomponent mass diffusion in reacting-flow sim-ulations, which may eliminate the need for the mixture-averaged assumption.However, despite these results, we do not have sufficient data to draw firm con-clusions on the accuracy and appropriateness of mixture-averaged assumptionsfor all flames (i.e., all fuels, configurations, and regimes). Additional data areneeded from studies of different fuels—namely large hydrocarbons—and kineticmodels with more species. Thus, future work should focus on extending thesecomparisons to other fuels and flame configurations.
Acknowledgements
This material is based upon work supported by the National Science Founda-tion under Grant No. 1314109-DGE. This research used resources of the NationalEnergy Research Scientific Computing Center (NERSC), a U.S. Department ofEnergy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.
References [1] R. B. Bird, W. E. Stewart, E. N. Lightfoot, Transport Phenomena, JohnWiley & Sons, Inc., New York, 1960.[2] J. Warnatz, Calculation of the Structure of Laminar Flat Flames I: FlameVelocity of Freely Propagating Ozone Decomposition Flames, Berichte der29unsengesellschaft für physikalische Chemie 82 (2) (1978) 193–200. doi:10.1002/bbpc.197800010 .[3] M. E. Coltrin, R. J. Kee, J. A. Miller, A Mathematical Model of SiliconChemical Vapor Deposition, Journal of The Electrochemical Society 133 (6)(1986) 1206. doi:10.1149/1.2108820 .[4] S. Lapointe, G. Blanquart, Fuel and chemistry effects in high Karlovitzpremixed turbulent flames, Combustion and Flame 167 (2016) 294–307. doi:10.1016/j.combustflame.2016.01.035 .[5] N. Burali, S. Lapointe, B. Bobbitt, G. Blanquart, Y. Xuan, Assessmentof the constant non-unity Lewis number assumption in chemically-reactingflows, Combustion Theory and Modelling 20 (4) (2016) 632–657. doi:10.1080/13647830.2016.1164344 .[6] J. Schlup, G. Blanquart, Validation of a mixture-averaged thermal diffu-sion model for premixed lean hydrogen flames, Combustion Theory andModelling 22 (2) (2018) 264–290. doi:10.1080/13647830.2017.1398350 .[7] T. P. Coffee, J. M. Heimerl, Transport algorithms for premixed, laminarsteady-state flames, Combust. Flame 43 (1981) 273–289. doi:10.1016/0010-2180(81)90027-4 .[8] J. Warnatz, Influence of transport models and boundary conditions onflame structure, in: N. Peters, J. Warnatz (Eds.), Numerical Methods inLaminar Flame Propagation, Vol. 6 of Notes on Numerical Fluid Mechanics,Vieweg+Teubner Verlag, Wiesbaden, 1982, pp. 87–111.[9] A. Ern, V. Giovangigli, Impact of detailed multicomponent transport onplanar and counterflow hydrogen/air and methane/air flames, Combus-tion Science and Technology 149 (1-6) (1999) 157–181. doi:10.1080/00102209908952104 .[10] H. Bongers, L. P. H. De Goey, The effect of simplified transport modelingon the burning velocity of laminar premixed flames, Combust. Sci. Technol.175 (10) (2003) 1915–1928. doi:10.1080/713713111 .[11] Y. Xin, W. Liang, W. Liu, T. Lu, C. K. Law, A reduced multicomponentdiffusion model, Combust. Flame 162 (1) (2015) 68–74. doi:10.1016/j.combustflame.2014.07.019 .[12] M. Faghih, W. Han, Z. Chen, Effects of Soret diffusion on premixed flamepropagation under engine-relevant conditions, Combustion and Flame 194(2018) 175–179. doi:10.1016/j.combustflame.2018.04.031 .[13] J. De Charentenay, A. Ern, Multicomponent transport impact on turbulentpremixed H /O flames, Combust. Theor. Model. 6 (3) (2002) 439–462. doi:10.1088/1364-7830/6/3/304 .3014] V. Giovangigli, Multicomponent transport in laminar flames, Proceedingsof the Combustion Institute 35 (1) (2015) 625–637. doi:10.1016/j.proci.2014.08.011 .[15] S. Dworkin, M. Smooke, V. Giovangigli, The impact of detailed multi-component transport and thermal diffusion effects on soot formation inethylene/air flames, Proceedings of the Combustion Institute 32 (1) (2009)1165–1172. doi:10.1016/j.proci.2008.05.061 .[16] T. Lu, C. K. Law, Toward accommodating realistic fuel chemistry in large-scale computations, Progress in Energy and Combustion Science 35 (2)(2009) 192–215. doi:10.1016/j.pecs.2008.10.002 .[17] G. Borghesi, J. Bellan, A priori and a posteriori investigations for de-veloping large eddy simulations of multi-species turbulent mixing underhigh-pressure conditions, Physics of Fluids 27 (3) (2015) 035117. doi:10.1063/1.4916284 .[18] E. Masi, J. Bellan, K. G. Harstad, N. A. Okong’o, Multi-species turbulentmixing under supercritical-pressure conditions: modelling, direct numericalsimulation and analysis revealing species spinodal decomposition, Journalof Fluid Mechanics 721 (2013) 578–626. doi:10.1017/jfm.2013.70 .[19] A. Ern, V. Giovangigli, Fast and Accurate Multicomponent TransportProperty Evaluation, Journal of Computational Physics 120 (1) (1995) 105–116. doi:10.1006/JCPH.1995.1151 .[20] A. Ern, V. Giovangigli, Thermal diffusion effects in hydrogen-air andmethane-air flames, Combust. Theor. Model. 2 (4) (1998) 349–372. doi:10.1088/1364-7830/2/4/001 .[21] S. Ambikasaran, K. Narayanaswamy, An accurate, fast, mathematicallyrobust, universal, non-iterative algorithm for computing multi-componentdiffusion velocities, Proceedings of the Combustion Institute 36 (1) (2017)507–515. doi:10.1016/j.proci.2016.05.055 .[22] O. Desjardins, G. Blanquart, G. Balarac, H. Pitsch, High order conservativefinite difference scheme for variable density low Mach number turbulentflows, Journal of Computational Physics 227 (15) (2008) 7125–7159. doi:10.1016/j.jcp.2008.03.027 .[23] B. Savard, Y. Xuan, B. Bobbitt, G. Blanquart, A computationally-efficient,semi-implicit, iterative method for the time-integration of reacting flowswith stiff chemistry, Journal of Computational Physics 295 (2015) 740–769. doi:10.1016/j.jcp.2015.04.018 .[24] J. O. Hirschfelder, C. F. Curtiss, R. B. Bird, Molecular Theory of Gasesand Liquids, Wiley, New York, 1954.3125] C. F. Curtiss, J. O. Hirschfelder, Transport properties of multicomponentgas mixtures, The Journal of Chemical Physics 17 (6) (1949) 550–555. doi:10.1063/1.1747319 .[26] J. F. Grcar, J. B. Bell, M. S. Day, The Soret effect in naturally propagat-ing, premixed, lean, hydrogen-air flames, Proceedings of the CombustionInstitute 32 (2009) 1173–1180. doi:10.1016/j.proci.2008.06.075 .[27] J. Schlup, G. Blanquart, A reduced thermal diffusion model for H and H ,Combustion and Flame 191 (2018) 1–8. doi:10.1016/j.combustflame.2017.12.022 .[28] R. Kee, F. Rupley, J. Miller, Chemkin-II: A Fortran chemical kinetics pack-age for the analysis of gas-phase chemical kinetics, Sandia National Labo-ratories Report SAND89-8009 (1989).[29] G. Dixon-Lewis, Flame structure and flame reaction kinetics. II. transportphenomena in multicomponent systems, Proceedings of the Royal Societyof London A: Mathematical, Physical and Engineering Sciences 307 (1488)(1968) 111–135. doi:10.1098/rspa.1968.0178 .[30] Y. Xuan, G. Blanquart, M. E. Mueller, Modeling curvature effects in diffu-sion flames using a laminar flamelet model, Combustion and Flame 161 (5)(2014) 1294–1309. doi:10.1016/j.combustflame.2013.10.028 .[31] Y. Xuan, G. Blanquart, Effects of aromatic chemistry-turbulence interac-tions on soot formation in a turbulent non-premixed flame, Proceedingsof the Combustion Institute 35 (2) (2015) 1911–1919. doi:10.1016/j.proci.2014.06.138 .[32] Y. Xuan, G. Blanquart, Two-dimensional flow effects on soot formationin laminar premixed flames, Combustion and Flame 166 (2016) 113–124. doi:10.1016/j.combustflame.2016.01.007 .[33] P. L. Carroll, G. Blanquart, A proposed modification to Lundgren’s physi-cal space velocity forcing method for isotropic turbulence, Physics of Fluids25 (10) (2013) 105114. doi:10.1063/1.4826315 .[34] S. Verma, Y. Xuan, G. Blanquart, An improved bounded semi-Lagrangianscheme for the turbulent transport of passive scalars, Journal of Computa-tional Physics 272 (2014) 1–22. doi:10.1016/j.jcp.2014.03.062 .[35] M. E. Mueller, H. Pitsch, LES model for sooting turbulent nonpremixedflames, Combustion and Flame 159 (6) (2012) 2166–2180. doi:10.1016/j.combustflame.2012.02.001 .[36] F. Bisetti, G. Blanquart, M. E. Mueller, H. Pitsch, On the formation andearly evolution of soot in turbulent nonpremixed flames, Combustion andFlame 159 (1) (2012) 317–335. doi:10.1016/j.combustflame.2011.05.021 . 3237] B. Savard, G. Blanquart, Broken reaction zone and differential diffusion ef-fects in high Karlovitz n -C H premixed turbulent flames, Combustionand Flame 162 (5) (2015) 2020–2033. doi:10.1016/j.combustflame.2014.12.020 .[38] M. Herrmann, G. Blanquart, V. Raman, Flux Corrected Finite VolumeScheme for Preserving Scalar Boundedness in Reacting Large-Eddy Simu-lations, AIAA Journal 44 (12) (2006) 2879–2886. doi:10.2514/1.18235 .[39] C. Pierce, Progress-variable approach for large-eddy simulation of turbulentcombustion, Ph.D. thesis, Stanford University (2001).[40] L. Shunn, F. Ham, P. Moin, Verification of variable-density flow solversusing manufactured solutions, Journal of Computational Physics 231 (9)(2012) 3801–3827. doi:10.1016/j.jcp.2012.01.027 .[41] O. Desjardins, G. Blanquart, G. Balarac, H. Pitsch, High order conservativefinite difference scheme for variable density low Mach number turbulentflows, Journal of Computational Physics 227 (15) (2008) 7125–7159. doi:10.1016/j.jcp.2008.03.027 .[42] R. D. Falgout, U. M. Yang, hypre: A Library of High Performance Precon-ditioners, in: P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, J. J. Dongarra(Eds.), Computational Science — ICCS 2002, Springer, Berlin, Heidelberg,2002, pp. 632–641. doi:10.1007/3-540-47789-6_66 .[43] H. A. van der Vorst, Bi-CGSTAB: A fast and smoothly converging variantof bi-CG for the solution of nonsymmetric linear systems, SIAM Journalon Scientific and Statistical Computing 13 (2) (1992) 631–644. doi:10.1137/0913035 .[44] M. Frigo, S. Johnson, The design and implementation of FFTW3, Pro-ceedings of the IEEE 93 (2) (2005) 216–231. doi:10.1109/jproc.2004.840301 .[45] L. F. Richardson, The approximate arithmetical solution by finite differ-ences of physical problems involving differential equations, with an appli-cation to the stresses in a masonry dam, Philosophical Transactions of theRoyal Society of London A: Mathematical, Physical and Engineering Sci-ences 210 (459-470) (1911) 307–357. doi:10.1098/rsta.1911.0009 .[46] F. Perini, E. Galligani, R. D. Reitz, A study of direct and Krylov iterativesparse solver techniques to approach linear scaling of the integration ofchemical kinetics with detailed combustion mechanisms, Combustion andFlame 161 (5) (2014) 1180–1195. doi:10.1016/j.combustflame.2013.11.017 .[47] Z. Hong, D. F. Davidson, R. K. Hanson, An improved H /O mechanismbased on recent shock tube/laser absorption measurements, Combustion33nd Flame 158 (4) (2011) 633–644. doi:10.1016/j.combustflame.2010.10.002 .[48] K.-Y. Lam, D. F. Davidson, R. K. Hanson, A shock tube study of H + OH → H O + H using OH laser absorption, International Journal of ChemicalKinetics 45 (6) (2013) 363–373. doi:10.1002/kin.20771 .[49] Z. Hong, K.-Y. Lam, R. Sur, S. Wang, D. F. Davidson, R. K. Hanson,On the rate constants of OH + HO and HO + HO : A comprehensivestudy of H O thermal decomposition using multi-species laser absorption,Proceedings of the Combustion Institute 34 (1) (2013) 565–571. doi:10.1016/j.proci.2012.06.108 .[50] D. G. Goodwin, H. Moffat, K., R. L. Speth, Cantera: An Object-orientedSoftware Toolkit for Chemical Kinetics, Thermodynamics, and TransportProcesses, Version 2.3.0 (2017). doi:10.5281/zenodo.170284 .URL [51] C. Rosales, C. Meneveau, Linear forcing in numerical simulations ofisotropic turbulence: Physical space implementations and convergenceproperties, Physics of Fluids 17 (9) (2005) 095106. doi:10.1063/1.2047568 .[52] S. Lapointe, B. Savard, G. Blanquart, Differential diffusion effects,distributed burning, and local extinctions in high Karlovitz premixedflames, Combust. Flame 162 (9) (2015) 3341–3355. doi:10.1016/j.combustflame.2015.06.001 .[53] F. A. Williams, Combustion Theory, Benjamin/Cummings, 1985.[54] M. Day, J. Bell, P. Bremer, V. Pascucci, V. Beckner, M. Lijewski, Tur-bulence effects on cellular burning structures in lean premixed hydro-gen flames, Combust. Flame 156 (5) (2009) 1035–1045. doi:10.1016/j.combustflame.2008.10.029 .[55] A. J. Aspden, M. S. Day, J. B. Bell, Turbulence-flame interactions in leanpremixed hydrogen: Transition to the distributed burning regime, J. FluidMech. 680 (2011) 287–320. doi:10.1017/jfm.2011.164 .[56] A. Aspden, M. Day, J. Bell, Turbulence-chemistry interaction in lean pre-mixed hydrogen combustion, Proc. Combust. Inst. 35 (2) (2015) 1321–1329. doi:10.1016/j.proci.2014.08.012 .[57] G. Blanquart, H. Pitsch, Thermochemical Properties of Polycyclic Aro-matic Hydrocarbons (PAH) from G3MP2B3 Calculations, J. Phys. Chem.A 111 (28) (2007) 6510–6520. doi:10.1021/JP068579W .[58] G. Blanquart, P. Pepiot-Desjardins, H. Pitsch, Chemical mechanism forhigh temperature combustion of engine relevant fuels with emphasis onsoot precursors, Combustion and Flame 156 (3) (2009) 588–607. doi:10.1016/j.combustflame.2008.12.007 .3459] K. Narayanaswamy, G. Blanquart, H. Pitsch, A consistent chemical mecha-nism for oxidation of substituted aromatic species, Combustion and Flame157 (10) (2010) 1879–1898. doi:10.1016/j.combustflame.2010.07.009 .[60] Y. He, B. Cook, J. Deslippe, B. Friesen, R. Gerber, R. Hartman-Baker, A. Koniges, T. Kurth, S. Leak, W.-S. Yang, Z. Zhao, E. Baron,P. Hauschildt, Preparing NERSC users for Cori, a Cray XC40 system withIntel many integrated cores, Concurrency and Computation: Practice andExperience 30 (1) (2017) e4291. doi:10.1002/cpe.4291 .[61] A. J. Fillo, J. Schlup, G. Beardsell, G. Blanquart, K. E. Niemeyer, Fig-ures, plotting scripts, and data for “A fast, low-cost, and stable memoryalgorithm for implementing multicomponent transport in direct numericalsimulations” [dataset], Zenodo (2019). doi:10.5281/zenodo.3519910 .[62] A. J. Fillo, J. Schlup, G. Beardsell, G. Blanquart, K. E. Niemeyer, Di-rect numerical simulation results for turbulent hydrogen/air flame: mul-ticomponent diffusion model [dataset], Oregon State University (2019). doi:10.7267/2f75rf80s .[63] A. J. Fillo, J. Schlup, G. Beardsell, G. Blanquart, K. E. Niemeyer, Directnumerical simulation results for turbulent hydrogen/air flame: mixture-averaged diffusion model [dataset], Oregon State University (2019). doi:10.7267/c247f0072 . Appendix A. Availability of material
The figures in this article, as well as the data and plotting scripts neces-sary to reproduce them, are available openly under the CC-BY license [61].Furthermore, the full simulation results from NGA are available for the three-dimensional multicomponent [62] and mixture-averaged [63] hydrogen/air flames.
Appendix B. Method verification
To verify the method implementation, we generated an artificial species pro-file where the direction and relative magnitudes of the flux vectors could bepredicted a priori to remain independent of any differential diffusion effects thatmay exist in a physical system. Specifically, we created a two-dimensional V-shaped species profile with a central angle of 45° and projected it into threedimensions as shown in Figure B.14a.Such a profile results in flux vectors that are constant in the y -direction, areof equal magnitude and opposite sign in the z -direction reflected over the x - y -plane, and vary in magnitude but remain constant in sign matching the initialinput profile in the x -direction. These predictions should be consistent indepen-dent of chemical species or other scalar value for the artificial input profile. Weran the algorithm for one “complete” set of sub-iterations to convergence and35 a) Input species profile(b) x -component of mass flux(c) y -component of mass flux(d) z -component of mass flux Figure B.14: Normalized flux vectors resulting from an artificial species profile after one fulliteration of semi-implicit multicomponent diffusion calculation. normalized the resulting diffusion flux vectors to ensure the relative magnitudesand direction were consistent with our expectations.Figure B.14 shows the results of this artificial test case. The resulting nor-malized flux vectors agree with expectation and have equal magnitudes in the x - and zz