Dimensional Decomposition of Turbulent Reacting Flows

Abstract

A dimensional-decomposition approach, decomposing 3D into 3×1D for turbulent (reacting) flows are motivated, discussed and investigated. In the three-dimensional Linear Eddy Model (LEM3D), three orthogonally intersecting arrays of 1D domains are coupled to capture the 3D characteristics of fluid flows. The currently used recouplings for LEM3D are enlightened and thoroughly discussed. A study of the flame front of a freely propagating laminar premixed flame shows that the flame stabilizes at the upstream face of the initial solution when both the advective and auxiliary recouplings are activated. Furthermore, results from LEM3D simulations of a vitiated co-flow burner are re-visited providing a more detailed discussion of the noted early mixing and reaction of the hydrogen fuel of the burner. The main conclusion of the present work is that the auxiliary coupling based on rotations of the 3D control volumes introduces very large gradients in the near-field geometry of jets, leading to a significant amount of artificial diffusion and locally increased burning rates. This implies that applications of LEM3D should be restricted to sub-regions where high-resolution treatment of scalar mixing and reaction is of particular interest.

1 Introduction

The key limitation in simulations of turbulent reacting flows is the computational cost. Ideally, the whole range of spatial and temporal scales should be resolved and captured. The computational cost stems from solving a combination of fundamental equations involving turbulent fluid flow, heat transfer, chemical reaction, radiative heat transfer and other complicated physical and chemical processes. A number of methods for turbulent mixing and reaction have been developed and investigated in order to meet the required physical description within the computational resources available. One among these is the dimensional-decomposition approach presented by Kerstein (2009), where it was proposed to decompose a three-dimensional flow into one-dimensional domains and recouple these to describe the physical processes. The motivation behind this is to obtain a fully resolved spatial and temporal resolution at reduced cost. In the original work, a line segment was defined to correspond to an edge of a linear stack of cubic control volumes. This approach, however, was later modified for numerical implementation since it is preferable to interpret the one-dimensional evolution as occupying a volume of space, enabling a finite-volume numerical representation rather than the more conceptual understanding of a line segment.

Such an approach, a 3×1D ’DNS’, will have the clear advantage of a much lower computational cost in comparison with traditional Direct Numerical Simulations (DNS). The cost reduction is a consequence of requiring fewer computational cells, solving down to the smallest resolution, to the Kolmogorov or the Batchelor scale, in one spatial direction only. DNS in quotes is used as neither 2D DNS nor 1D DNS in principle are able to resolve all turbulent scales. Figure 1 illustrates the number of computational cells for the dimensional-decomposition approach (3×1D) and traditional DNS for Reynolds numbers of $10^4$ and $10^5$.

Fig. 1

Number of computational cells in conventional DNS versus the number of computational cells needed with the dimensional-decomposition approach

In turbulent reactive flows, one-dimensional ’DNS’ has no perception of turbulent mixing. The Linear Eddy Model (LEM) formulated by Kerstein (1988, 1991, 1992) and the One Dimensional Turbulence Model (ODT), also developed by Kerstein (1999, 2002), are one-dimensional models where turbulent eddies are represented through a special measure-preserving map called the triplet map. The discrete implementation of the triplet map is a permutation of fluid elements, also named wafers. Hence, LEM/ODT is a natural choice for the dimensional-decomposition strategy. Several models have been developed based on this strategy which all incorporate a three-array structure of LEM or ODT domains, e.g. ODTLES (Schmidt et al. 2005, 2010), LEM3D (Sannan et al. 2013; Weydahl 2010), and Lattice-Based Multiscale Simulation (LBMS) (Sutherland 2018). The focus in recent years for this strategy has been on the coupling of arrays of ODT/LEM domains, so as to obtain a self-contained 3D flow simulation with the smallest scales resolved in each 1D domain.

An alternative approach to the dimensional-decomposition strategy is the LEMLES model in which LEM is used as a sub-grid scalar closure for LES (Menon and Calhoon 1996; Chakravarthy and Menon 2000). In LEMLES, the unresolved length and time scales of LES are captured by individual LEM domains embedded in each LES cell. The LEM domains communicate with each other by means of a ’splicing’ procedure, by which portions of the domains are extracted and pasted onto domains in the adjacent LES cells. The advective transport of fluid elements across each LES cell interface by the splicing is determined by the LES-resolved mass fluxes across the cell interfaces. A survey of LEMLES applications is given by Menon and Kerstein (2011).

The spatial structure of LEM3D has been described in detail by Sannan et al. (2013) and Weydahl (2010) and is different from that of LEMLES. The structure is the same as in ODTLES and involves three mutually orthogonal arrays of 1D domains that each spans the three-dimensional flow domain in one direction. One difference from LEMLES is therefore that the 1D domains are not constrained to 3D control volumes, and that the turbulence-emulating triplet maps are not limited in size to such control volumes. The 1D domains must be coupled, however, in order to capture the 3D characteristics of fluid flows. In ODTLES, the domains are coupled through a two-way procedure in which temporal and spatial filtering enable an up-scale coupling from ODT to LES, and a subsequent down-scale procedure that reconstructs the ODT velocity field based on the 3D continuity equation (Schmidt et al. 2010). LBMS, on the other hand, is based on a lattice-like structure of 1D ODT models which are coupled through a particular reconstruction treatment of the fluxes perpendicular to the ODT lines (Sutherland 2018). In comparison, LEM3D is formulated with an advective coupling in addition to an auxiliary coupling, as explained in Sect. 2.1.

A recent study by Grøvdal et al. (2020) addressing issues connected to the dimensional recoupling for LEM3D revealed certain issues connected with the dimensional-decomposition approach. Specifically, an overestimation of the molecular diffusion was found in areas strongly dominated by flow in one direction, e.g., in regions near the nozzle exit of a turbulent reactive jet. One aim of the present paper is to gain a deeper understanding of the previously observed artifact and to assess its implications. Since the previous studies of a turbulent flame using LEM3D dealt with a composite flow configuration (Grøvdal et al. 2018, 2020), i.e., the lifted hydrogen jet flame of Berkeley’s vitiated co-flow burner (Cabra et al. 2002), we here first consider the simple configuration of a freely propagating flame to validate the diffusion and chemical kinetics implementation of the code.

It should be emphasized that LEM3D in its present formulation is used as a post-processing tool to RANS. In this modeling approach, mean-flow information from RANS (or a global flow solver) provides model input to LEM3D, which complements RANS with unsteadiness and scalar statistics needed for more accurate mixing and reaction calculations. The capability of LEM3D to provide any type of computational cost savings for equivalent accuracy (or enhanced accuracy for equivalent cost) in comparison with established turbulent combustion (or turbulence-chemistry-interaction) models is not attempted nor discussed in the present work. Furthermore, since only simple flame cases are considered here, it remains to be seen how accurate and robust LEM3D would be in realistic flows such as those in piston engines and gas turbines.

The paper is organized as follows. Section 2 gives a brief introduction of the LEM3D model, describing the recoupling mechanisms in detail. Section 3 introduces the two flow configurations considered; a numerical setup based on a freely propagating flame and the other for the near-field study of a jet nozzle, here given by the vitiated co-flow burner (Cabra et al. 2002). The results and a discussion of these are contained in Sect. 4, while the conclusions of the paper are given in Sect. 5.

2 The LEM3D Model

The LEM3D model has previously been described in detail by Sannan et al. (2013), Weydahl (2010), and Grøvdal et al. (2018, 2020). In LEM3D, three orthogonally intersecting arrays of 1D LEM domains are coupled so as to capture the 3D character of fluid trajectories. The LEM3D flow domain is a cuboid discretized in two distinct but interdependent ways, as illustrated in Fig. 2. The cubic control volumes (3DCVs) of the Cartesian mesh is formed by the intersections of sets of three orthogonal LEM domains. The LEM domains thus provide fine-scale resolution in all three spatial directions to the coarser 3D mesh.

Fig. 2

The LEM3D flow domain discretized by a coarse Cartesian mesh and a superimposed fine-scale resolution by orthogonal arrays of LEM domains in each spatial direction. The cubic control volumes (3DCVs) correspond to intersections of sets of three orthogonal LEM domains. The superimposed fine-scale resolution is illustrated by the colored LEM domains in blue, green and red. One domain is illustrated in each coordinate direction and these intersect in the top-front corner control volume. Note that the simulation domain in general is discretized by a much larger number of 3DCVs and LEM domains than shown in the figure. Also for illustrative purposes, the fine-scale resolution (number of LEM wafers within each 3DCV) is kept low and is much higher in actual simulations

In the brief introduction presented here, we will focus on the recoupling mechanism of the 1D domains. The evolution of the species mass fractions $Y_k$ on an LEM domain oriented in the Cartesian $x_j$-direction is governed by the equation

$$\begin{aligned} \frac{\partial \rho Y_k}{\partial t} + \frac{\partial \,\overline{\rho\:\! u_\alpha } Y_k}{\partial x_\alpha } = \frac{\partial }{\partial x_j}\left[ \rho D_k \frac{\partial Y_k}{\partial x_j}\right] + TM_j + \rho\;\!\omega _k, \end{aligned}$$

(1)

where $\rho$ is the density, $D_k$ is the molecular diffusivity of species k, and $\omega _k$ is the chemical reaction rate. The index j indicates that the terms are implemented on 1D LEM domains in all three spatial directions. Note that the conventional summation over the repeated indices j and k does not apply to the right-hand side of the equation. The term $\partial \,\overline{\rho\:\! u_\alpha } Y_k/\partial x_\alpha$ represents the Lagrangian advection of wafers on the domain, including the coupling to intersecting domains as described in Sect. 2.1. The repeated index $\alpha$ is summed over to reflect the 3D nature of the coupling of the LEM domains. Similarly, the governing transport equation for the specific enthalpy h is written as

$$\begin{aligned} \frac{\partial \rho h}{\partial t} + \frac{\partial \,\overline{\rho\:\! u_\alpha }h}{\partial x_\alpha } = \frac{\partial }{\partial x_j}\left[ \rho \kappa \frac{\partial h}{\partial x_j}\right] + TM_j + S_h, \end{aligned}$$

(2)

where $\kappa$ is the thermal conductivity and $S_h$ is the heat source term due to chemical reaction. The term $TM_j$ in the above equations symbolically denotes randomly occurring stirring events (triplet maps) that emulate the effects of turbulent eddies on the scalar concentration fields. The respective advection terms $\partial \,\overline{\rho\:\! u_\alpha } Y_k/\partial x_\alpha$ and $\partial \,\overline{\rho\:\! u_\alpha } h/\partial x_\alpha$ are governed by an averaged mass flux $\overline{\rho\:\! u_{\alpha }}$ which is prescribed to the model from a global flow solver or measurements.

The reactive-diffusive and turbulent transport of Eqs. (1) and (2) is time advanced for each LEM domain of the three orthogonal arrays of domains. The time advancement involves a subcycling of smaller time steps within a coarser advective time step $\varDelta t$. The subcycling is punctuated by the randomly occurring triplet maps, whose rate is determined by the locally given turbulent diffusivity $D_T$ (model input).

By the reduced description of the scalar fields to one spatial dimension, a full-scale resolution of all scales, down to the smallest turbulence scales, is computationally feasible. Hence, the scalars $Y_k$ and h are determined both spatially and temporally down to the Kolmogorov scale. The fine-scale resolution is such that the LEM wafers can be considered as homogeneous reactors, and the chemistry is therefore implemented directly without the need for any modeling. In the current formulation, detailed and finite-rate chemistry is solved using the CHEMKIN II software package. Heat release and thermal expansion, i.e., dilatation, is accounted for by increasing the wafer volume and performing a regridding prior to the transport operations of every advective time step.

Fig. 3

Illustrative example of the advective coupling. The blue, red and green fluid elements illustrate the wafers in the three spatial directions, while the arrows indicates the fluxes into/out of the cubic control volumes (3DCVs). Based on the prescribed mass fluxes, a single wafer is advected into the 3DCV from each side of the blue-domain direction, a single wafer is advected into and out of the red-domain direction, while a single wafer is advected into the green-domain direction and three are advected out. The three central figures illustrate the need for advective flipping of wafers due to the constraint of Eq. (3); the green domain lacks two wafers and the blue domain have two wafers in excess. These two excessive wafers are extracted from the center of the blue domain, flipped, and then inserted into the green domain to fill the gap created by the advection of that domain

2.1 Dimensional Recoupling

The recoupling mechanisms are associated with the advective time step $\varDelta t$ and consist of two operations:

• 1. Advective coupling provided by the terms ${\partial (\overline{\rho\:\! u_\alpha }Y_k)}/{\partial x_\alpha }$ and ${\partial (\overline{\rho\:\! u_\alpha }h)}/{\partial x_\alpha }$ of Eqs. (1) and (2), respectively. The displacements of LEM wafers are based upon the prescribed mass fluxes $\overline{\rho\:\!\varvec{u}}$ and the implementation of the continuity equation in every 3DCV, while accounting for dilatations from reactions due to the source terms $\rho\;\!\omega _k$ and $S_h$. Typically, this involves transfers of wafers among the differently oriented LEM domains intersecting each cubic control volume. A constant number of wafers is enforced by requiring the displacements of wafers into and out of each 3DCV to obey the equation

  $$\begin{aligned} \sum _{l=1}^{N_{\mathrm{faces}}} \delta _l = 0, \end{aligned}$$

  (3)

  where $\delta _l$ denotes the integer number of wafer displacements across the 3DCV faces. A 3D illustration of the advection operation is shown in Fig. 3.

• 2. An auxiliary coupling is implemented by stochastic rotations of the 3DCVs. The rotations give additional fluid exchanges between the LEM domains, and ensure that physical processes are consistently represented in all three spatial directions. The rotation mechanism is illustrated in Fig. 4. Thus, for every advective time step $\varDelta t$, the 3DCVs are rotated $\pm \, 90^\circ$ about any of the three coordinate axes with a locally defined probability

  $$\begin{aligned} p_{\mathrm{rot}}=\frac{3}{2}\,C_{\mathrm{rot}}\cdot \,{\hbox {CFL}}_{\mathrm{3DCV}}, \end{aligned}$$

  (4)

  where ${\hbox {CFL}}_{\mathrm{3DCV}}$ is the local Courant–Friedrichs–Lewy number, and $C_{\mathrm{rot}}$ is a model constant. In the present study we have set $C_{\mathrm{rot}}=1$, corresponding to that, on average, all wafers will be swapped to differently oriented LEM domains once per advective residence time within the 3DCVs (Sannan et al. 2013). Note that if a rotation occurs, the axis of rotation is selected with an equal probability among the three coordinate axes, and there is an equal probability for rotating clockwise or counterclockwise. Thus, there are six different options for performing a rotation, all with an equal probability. Since the rotations introduce diffusive transport, the induced diffusivity $D_T^{\mathrm{rot}}$ is deducted from the turbulent diffusivity $D_T$. $D_T$ is prescribed from experiments or calculated from the $k\text{-}\varepsilon$ turbulence model by $D_T=C_\mu k^2/\sigma _{T}\,\varepsilon$, where $C_\mu =0.09$ and $\sigma _T$ is the turbulent Schmidt number. The effective triplet map diffusivity $D_T^{TM}$ is then given by

  $$\begin{aligned} D_T^{TM}=D_T-D_T^{\mathrm{rot}}, \end{aligned}$$

  (5)

  and is used as input to the temporal Poisson distribution f(t) of triplet maps by setting the mean number of events per LEM domain (Weydahl 2010).

Fig. 4

The auxiliary coupling between the LEM domains illustrated by the rotation between a blue domain and a green domain of wafers. As observed, the red domain wafers are not affected by this particular rotation

Due to the increase of volume by heat release from reactions, following the ideal gas law for constant pressure, there is an associated deviation $\varXi _{\mathrm{flux}}$ for each 3DCV between the prescribed mass flux and the actual mass fluxed by LEM3D given by

$$\begin{aligned} \varXi _{\mathrm{flux}}=\sum _{l=1}^{N_{\mathrm{faces}}}\left( \varDelta x_w\sum _{j=1}^{N_l}\rho _{lj}-\overline{\rho \,u}_l\varDelta t\right) , \end{aligned}$$

(6)

where $\varDelta x_w$ is the wafer size, $N_l$ is the number of wafers fluxed across face l, and $\rho _{lj}$ denotes the density of the wafer j fluxed into or out of the 3DCV at face l. The given deviation is an integral part of the current implementation in the variable density framework and in the presence of chemical reactions. This deviation was not present in previous formulations of the model due to the assumption of constant density (Sannan et al. 2013). In order to adapt the deviation $\varXi _{\mathrm{flux}}$ into the algorithm, an investigation of various iteration procedures, in combination with least-squares methods to minimize the deviations at each 3DCV cell face, was previously performed by Grøvdal et al. (2020). It was shown that a least-square minimization based on a breadth-first search, with the root given where the reactions initially occur, resulted in the most accurate and comparable flame shape to experiments for the co-flow burner. In the current paper we provide a more detailed discussion and understanding of the effects of the dimensional recoupling, based on simulations using the breadth-first iteration procedure.

3 Numerical Setup

This section describes the numerical setup for the flow configurations to be investigated. For both the in-house LEM3D code and the commercially available LOGEresearch^(TM) code (Loge), the chemical reaction mechanism employed is the detailed ${\mathrm{H}}_2/{\mathrm{O}}_2$ mechanism of Li et al. (2004).

3.1 Freely Propagating Hydrogen Flame

A simple laminar freely propagating ${\mathrm{H}}_2$ flame is set up by use of LOGEresearch, with temperature, density, and mole fractions as given in Table 1. The resulting flame speed is 2.265 m/s and given as velocity in Table 1. The velocity and mass flux is subsequently prescribed to the LEM3D code together with species profiles for initialization. All turbulent parameters in LEM3D are set to zero in these simulations since the flow is laminar.

Table 1 Flow boundary conditions prescribed to LEM3D for the freely propagating hydrogen flame

3.2 Vitiated Co-flow Burner

The vitiated co-flow burner, developed at UC Berkeley by Cabra et al. (2002), is a lifted turbulent H₂/N₂ jet flame with a co-axial flow of hot combustion products from lean premixed H₂/air flames. The co-flow flames are stabilized on a perforated disk with 87% blockage and an outer diameter of 210 mm. The central jet exit diameter is 4.57 mm and extends 70 mm above the surface of the perforated disk.

The burner was simulated in a previous LEM3D study for various configurations by Grøvdal et al. (2020). The input parameters for the initial RANS simulation, which provides input to LEM3D, is given by Cabra et al. (2002) and Myhrvold et al. (2006). The RANS simulation was conducted with a modified $k\text{-}\varepsilon$ model using the academic ANSYS Fluent package. A full description of the set-up can be found in Grøvdal et al. (2020).

4 Results and Discussion

4.1 Freely Propagating Hydrogen Flame

As a validation of the diffusive-reactive implementation of the LEM3D code, a simulation of a freely propagating ${\mathrm{H}}_2$ flame without any dimensional recoupling was conducted, considering a single LEM domain oriented in the flow direction. It should be emphasized that since the investigated flow is laminar, a key feature of the Linear Eddy Model, namely the triplet map, is disabled for this study. The flame position of the freely propagating premixed laminar flame is determined by the initialization alone if the inlet velocity is equal to the laminar flame speed. LEM3D is typically initialized using mean-flow information either from a preceding RANS simulation or from measurements. Therefore, as a comparison, we first initialize the flame at the coarse control volume level for a simulation. Then, using the LOGEresearch simulation tool, a solution for the simple freely propagating flame is obtained on a much finer scale, which is used to initialize LEM3D at the wafer level. Figure 5 shows the resulting density profiles, indicating the flame location for the two different initialization procedures. The black and blue-dotted curves in the sub-plots show the given initialization and the LEM3D solution, respectively, while the dashed light gray lines represent the control volume boundaries. For the left sub-plot, the wafers are initialized at the control volume level, i.e., by the step function plotted in black. The size of the control volume is here $\varDelta x=0.02$ m. In this case, since the gradients are too sharp for an immediate stabilization, the flame is shifted slightly downstream as the preheat zone is building up and reaches steady state once this zone is established. For the right sub-plot, the wafers are initialized on the wafer level based on the solution of a LOGEresearch simulation. The fine-scale resolution is defined by 160 wafers within each control volume, giving a wafer size of $\varDelta x_w=1.25\times 10^{-4}$ m. The LEM3D profile now stabilizes at the same flame-front location as given by the LOGEresearch code. Thus, while the initialization for the left-plot simulation corresponds to a typical input from RANS, we observe that the individual fine-scale initialization of the wafers results in a profile in full agreement with the LOGEresearch simulation. The actual fine-scale resolution used is somewhat arbitrary, but testing has shown that higher resolutions give no difference in the results for this laminar flame application. Hence, since the initial conditions, chemical mechanism, and transport equations are identical, the simulation demonstrates that the diffusion and chemical kinetics implementation of LEM3D is in agreement with that of the LOGEresearch tool.

Fig. 5

Freely propagating laminar flame position for two different initializations. The black and blue-dotted curves show the initial and LEM3D density profiles, respectively, while the dashed light gray lines indicate the control volume boundaries. For the left plot the flame is initialized at the level of the coarse control volume of size $\varDelta x=0.02$ m, while for the right plot the individual wafers of size are initialized by the fine-scale solution of a LOGEresearch simulation of the laminar freely propagating flame

Next, four simulations using both the dimensional recouplings of LEM3D are conducted. The setups for these are given in Table 2, where $\varDelta x$ denotes the size of the control volume and $N_i$ is the number of control volumes in the streamwise direction (denoted the I-direction). The simulation denoted Coarse refers to the use of coarse 3DCVs with step initialization, and the computational domain is given by the interval [0.020, 0.100] m. The second simulation, denoted Shifted, makes use of the same $\varDelta x$ and initialization as Coarse, but the domain is shifted 0.005 m downstream. This is done to investigate the effect on the flame-front location with respect to cell faces. The third simulation, denoted Individual, is similar to Coarse with respect to $\varDelta x$ and the domain, but here the wafers are initialized individually from the results given by the LOGEresearch tool. Finally, the fourth simulation, denoted Fine, makes use of the same initialization and domain as Coarse, but the 3DCVs are refined 20 times. The wafer size is kept constant for all four simulations, with $\varDelta x_w=1.25\times 10^{-4}$ m and a total of 640 wafers in the streamwise direction of the computational domain. This corresponds to a resolution of 160 wafers in each coordinate direction within each 3DCV for the Coarse, Shifted, and Individual simulation, and in the same manner 8 wafers for the Fine simulation.

Table 2 Simulation cases where either both the LEM3D dimensional recouplings are used (results shown in Fig. 6) or just the advective coupling (results in Fig. 7)

Fig. 6

Time-averaged density profiles of LEM3D with both the advective and the auxiliary coupling for the four cases of Table 2. The black and blue-dotted curves show the given initialization and the LEM3D solution, respectively, while the dashed light gray lines indicate the control volume boundaries. The coarse mesh of size $\varDelta x=0.02$ m is resolved by 160 wafers in each coordinate direction, while the fine-mesh control volume has a width of $\varDelta x=0.001$ m and is resolved by 8 wafers in each coordinate direction

Figure 6 shows the initial (in black) and the time-averaged LEM3D density profiles (blue-dotted curve) for each of the cases given in Table 2. The control volume boundaries are indicated by the dashed light gray lines. The time averaging is taken over a period of 5 sec, subsequent to an initial relaxation time and corresponding to about 140 flow-through times for the simulation domain. We observe that the flame stabilizes close to the cell face of initialization for all four cases of Fig. 6, and that there are some ripples downstream of the flame front for the time-averaged profiles. The stabilization of the flame front is influenced by the rotations, which represent the dominating recoupling mechanism between the I-oriented domain (streamwise direction), and the J- and K-oriented domains (lateral directions). Dissimilar fluid states are regularly brought into contact by the rotations, which create large gradients and pockets within the reaction zone. In some cases hot and burning wafers are rotated to the upstream cell face, with the result that the flame front is shifted towards this location. The dominant role of the rotations is especially supported by the third sub-plot of Fig. 6, which clearly shows that the flame front has moved upstream from the point of local initialization and stabilizes near the closest up-stream cell face. In other events, unburned mixtures are moved to the upstream face, and the created pockets lead to the ripples seen in all the sub-plots. This is most clearly observed in the rightmost sub-plot of Fig. 6 since the resolution here is higher.

From the Coarse and Individual simulation plots of Fig. 6, we conclude that the flame front is independent of the fine details of the initialization. We also note, based on all four sub-plots, that the flame stabilizes within the control volume where it is initialized and close to the upstream face of that volume. The fact that the flame front stabilizes within the control volume of initialization is as expected in this case since the flame speed is balanced by the flow velocity. Also, there is no action of turbulent eddies (triplet maps) in the laminar freely propagating flame configuration. The main observation from this laminar flame study, however, is that the stochastic rotations cause large gradients and pockets of unburned wafers in the reaction zone, leading to a higher consumption rate and an eventual stabilization of the flame front at the upstream cell face of the control volume of initialization.

Fig. 7

LEM3D density profiles with the auxiliary coupling switched off for the four cases in Table 2. The black and blue-dotted curves show the given initialization and the time-averaged LEM3D solution, respectively, while the red curve shows an instantaneous LEM3D solution. The dashed light gray lines indicate the control volume boundaries. The coarse mesh contains 160 wafers and has a control volume width of $\varDelta x=0.02$ m, while the fine-mesh control volume has a width of $\varDelta x=0.001$ m and contains 8 wafers

We finally investigate the effect of switching off the auxiliary coupling of rotations to consider the LEM3D advective coupling by flipping of wafers only. To this end, four additional simulations are performed with setups as given for the cases of Table 2. The results of these simulations are shown in Fig. 7, with the initial and time-averaged density profiles shown by the black and blue-dotted curve, respectively, and an additional red curve (compared to Fig. 6) showing an instantaneous density profile of the freely propagating flame configuration. We now observe for all four simulations that the flame is shifted to its closest control-volume center, based on the initialization. Thus, the stabilization of the flame front is again independent of the initialization procedure. Furthermore, the instantaneous profiles of the four sub-plots contain spiky fluctuations within the flame front, relative to the mean profile (blue-dotted curve). These spikes result from the instantaneous flipping of wafers by which these are inserted into the middle of the flame location, splitting the flame front. Since the wafer contents of the the J- and K-oriented domains generally will differ, the instantaneous profiles can vary greatly when different wafers from these domains are flipped into the I-oriented domain. From this we deduce that the advective flipping of wafers also cause large gradients which, in the absence of the rotational recoupling, causes the flame front to stabilize at the center of the control volume of initialization.

Fig. 8

Illustration of the rotational and advective coupling for the freely propagating flame in the control volume of the flame front. Individual wafer values are shown, ranging from blue (unburned) to yellow (fully combusted), for the three spatial orientations at three consecutive advective time steps. To the left is shown the instantaneous solution at some initial time $t_0$. At $t_0+\varDelta t$, while some of the wafers in the I-direction have diffused and reacted, a few unreacted wafers from the center of the J-direction are inserted by the advective flipping into the center of the I-direction. At $t_0+2\varDelta t$, the J-oriented wafers are rotated into the K-direction and, oppositely, the K-oriented wafers are rotated into the J-direction. Note that the K-oriented domain in reality is oriented perpendicular to the paper plane but is here oriented in the same direction as the J-oriented domain for illustrative purposes

As an illustration of the LEM3D processes, the instantaneous distribution of unburned and burned wafers for all three orientations in the flame-front control volume are shown for three sequential advective time steps in Fig. 8. The K-oriented domain is oriented perpendicular to the plane of paper but for illustration purposes the wafers are here arranged in the same direction as the J-oriented domain since they both are lateral. The flame front is located in the I-oriented domain as this is the flow direction. Note that physically we are considering a one-dimensional flow, and hence, there is no massflow in the J- and K-oriented domains. The flame fronts in the J- and K-direction reside within the respective control volumes but are not resolved in the flow direction. Between $t_0$ and $t_0+\varDelta t$, the burning wafers progress towards fully combusted states and gradients are smeared out by molecular diffusion. At $t_0+\varDelta t$, a few wafers are flipped from the center of the J-oriented domain into the center of the I-oriented domain, as indicated by the dashed arrow in Fig. 8. This is due to heat release and a surplus volume in the J-direction, which in the current case implies that additional fluid elements are pushed downstream in the I-direction after the flipping of the wafers. Between $t_0+\varDelta t$ and $t_0+2\varDelta t$, the wafers again react and diffuse. At $t_0+2\varDelta t$, the wafers oriented in the J- and K-direction swap orientation by the auxiliary rotational coupling. This is indicated by the crossed arrows between $t_0+\varDelta t$ and $t_0+2\varDelta t$ for the J- and K-oriented domains. Since adjacent control volumes are not shown in Fig. 8, the effect of the rotation creating new and sharp gradients at the cell faces are strictly not highlighted in this illustration. One may imagine, however, that the wafers from either of the domains are rotated into the I-oriented domain with fully burnt wafers upstream and unburned wafers downstream. This would clearly create the large gradients stabilizing the flame at the upstream cell face.

From the current study, it is noted that the base of the freely propagating flame does not move out of the control volume from which it is initialized. As observed, the flame stabilizes either at the center of the control volume, or close to or at the upstream cell face of the control volume, depending on whether the auxiliary coupling has been switched off or not. The stabilization of the flame close to the upstream cell face is due to the effect of the rotations from the auxiliary coupling. With the auxiliary coupling switched off, however, the flame stabilizes at the center of the control volume due to the advective coupling. By the flipping operation, excessive wafers in one domain is extracted from the center of the control volume and inserted into the center of the same control volume in a domain for which there is a deficit of wafers. Thus, since dissimilar fluid states are regularly brought into contact by the flipping of wafers in this case, large gradients are created close to the center of the control volume leading to the stabilization of the flame front at this location. Both the number of flipped wafers and their contents may vary greatly, leading to the large fluctuations within the flame front as observed in the instantaneous profiles of Fig. 7. However, these fluctuations are shown to decrease with the wafer size as the variation in the number of flipped wafers is smaller in the case of higher resolution.

4.2 Vitiated Co-flow Burner

The LEM3D code has previously been applied to the Berkeley vitiated co-flow burner by Grøvdal et al. (2018, 2020). In these studies, the LEM3D simulations showed good agreement with experimental data by Cabra et al. (2002) in the far field, but the current implementation did not yield a lift-off height of the hydrogen jet flame of the burner in accordance with the experiments. The proposed explanation for this was that the auxiliary rotational coupling caused an enhanced spreading of ${\mathrm{H}}_2$ in the near field, leading to early mixing and reaction of the hydrogen fuel and a stabilized flame with a low liftoff. However, this hypothesis was not validated by a thorough analysis of the near field scalar fields. In this section, we follow up on this issue with a more detailed discussion of the effect of the auxiliary coupling on the flow field, especially in the region close to the jet nozzle of the burner. For simplicity, we here present the temperature profiles only since these represent both mixing and reaction. In comparison, the mixture fraction only reflects the mixing of the species.

In the present modeling approach, LEM3D is used as a post-processing tool for an initial RANS simulation. The RANS simulation is here performed on a cuboidal Cartesian $75\times 75\times 60$ grid employing a modified $k\text{-}\varepsilon$ model. The RANS model input to LEM3D is mean-flow information such as the average mass flux $\overline{\rho\:\! \mathbf {u}}$ and the turbulent diffusivity $D_T$ of the flow field. In the RANS simulation, the jet inlet is approximated by a single grid cell such that the area of the jet is preserved, i.e., $(\varDelta x)^2=\pi (d/2)^2$, where $d=4.57$ mm is the jet diameter and $\varDelta x=4.05$ mm is the grid size. The effect of using a square nozzle is assumed to be negligible, justified by an experimental study of a similar burner configuration which indicated no detectable impact on the lift-off height of the flame from changes in the nozzle geometry (North 2013). The LEM3D simulation domain is a cuboidal $31\times 31\times 50$ grid and a sub-domain of the RANS domain. The 3DCVs of the sub-domain are equal in size to the grid cells of the RANS domain and coincide with the corresponding RANS cells. Hence, no interpolation is needed and the values of the turbulent diffusivity and the face-normal mass fluxes from RANS are used as direct input to LEM3D. The 3DCVs are resolved by 78 LEM wafers in each coordinate direction, giving a total of 11243700 $(31\cdot 31\cdot 50\cdot 78\cdot 3)$ wafers for the simulation of the vitiated co-flow burner.

Figure 9 shows the central upstream layer of the LEM3D cuboidal mesh, with the gray area indicating the jet nozzle inlet. The thick orange and blue line represent the J- and K-oriented LEM domains through the upstream centerline grid cell, denoted as ${\mathrm{J}}_r$ and ${\mathrm{K}}_r$, respectively, to indicate that profiles along these lines are radial profiles. The I-oriented domains are perpendicular to the paper plane, intersect the blue and orange lines at the center of the grid cells, and are represented by black dots in the figure. Instantaneous temperature profiles for the ${\mathrm{J}}_r$ and ${\mathrm{K}}_r$ domains, as well as the I-oriented domains ${\mathrm{I}}_0$, ${\mathrm{I}}_1$, ${\mathrm{I}}_2$, ${\mathrm{I}}_3$, and ${\mathrm{I}}_4$, are shown in Fig. 10. The bottom plot shows the profiles for the ${\mathrm{J}}_r$ and ${\mathrm{K}}_r$ domains in the interval $r/\varDelta x \in [-\,1.5,1.5]$, where the ’radial’ coordinate r jointly denotes either the j- or k-coordinate. The upper sub-plots of Fig. 10 show the instantaneous temperature profiles for the I-oriented (axial) LEM domains for $z/\varDelta x \in [0,1]$. The profiles have been rotated to illustrate the relative orientation of the axial domains compared to the ${\mathrm{J}}_r$ and ${\mathrm{K}}_r$ domains. The I-oriented profiles are given the color of the thick-colored domain intersecting the corresponding domain (cf. Fig. 9), except for the center profile of ${\mathrm{I}}_0$ which is intersected by both the orange and blue domain and given by the average RGB color of these colors. The black lines of all the plots denote the initial RANS profile prescribed from an ANSYS Fluent simulation.

Fig. 9

A 2D illustration of the central upstream layer of the LEM3D cuboidal mesh. The middle gray area represents the control volume of the jet nozzle inlet, and the dashed light gray lines give the boundaries of the neighboring control volumes. The orange and blue lines indicate J- and K-oriented LEM domains, respectively, with the thick lines denoted by ${\mathrm{J}}_r$ and ${\mathrm{K}}_r$ representing the domains through the upstream centerline grid cell. The black dots represent the I-oriented (axial) domains which are perpendicular to the paper plane. Instantaneous temperature profiles, corresponding to the denoted domains and associated colors, are given in Fig. 10

Fig. 10

Instantaneous temperature profiles for axially and radially oriented LEM domains in the central upstream layer of grid cells for the vitiated co-flow burner. The bottom plot shows the profiles of the radially oriented ${\mathrm{J}}_r$ and ${\mathrm{K}}_r$ domains illustrated in Fig. 9, with coloring following the colors of the domains in that figure. The upper sub-plots show the profiles of the axially oriented domains ${\mathrm{I}}_0$, ${\mathrm{I}}_1$, ${\mathrm{I}}_2$, ${\mathrm{I}}_3$, and ${\mathrm{I}}_4$ shown in Fig. 9, with coloring according to the color of the thick-colored radial domain intersecting the corresponding domain. For the upper center sub-plot, the ${\mathrm{I}}_0$ profile is colored by the average RGB color of the orange and the blue domain intersecting the upstream centerline grid cell. The black profiles of the plots show the prescribed solution from an initial RANS simulation

We note that there are unphysically large gradients in the radial plots outside the jet nozzle (central) control volume, i.e., for $|r|/\varDelta x \ge 0.5$. Also, the temperature hits very high values of about 1500 K in parts of the control volumes adjacent to the central control volume, which is much higher than the expected co-flow temperature of 1045 K given by the black profiles of Fig. 10. Hence, chemical reactions and heat release appear to take place already in the upstream first layer of control volumes, just inside the co-flow inlet of the cuboidal flow domain. This is surprising but falls in line with the previous result of early downstream burnout of the hydrogen fuel, as pointed out in the discussion of the axial scatter plots by Grøvdal et al. (2020).

Fig. 11

Conceptual illustration describing the creation of the gradients seen in Fig. 10 in control volumes adjacent to the central control volume of the jet inlet. In the upper plots, an initial solution is smeared out by diffusion from the left to the central plot. Then, the radial profile is rotated into the streamwise direction in the upper right plot. In the lower left plot, the profile is shifted downstream by pure advection. In the lower central plot, the profile is rotated back into a radial direction. Finally, in the lower right plot, local temperature maxima and minima are developed due to molecular diffusion and reaction

To explain the early mixing and reaction of the hydrogen fuel in the simulations of the vitiated co-flow burner, a series of conceptual sketches is presented in Fig. 11. The sketches are based on results from analyzing transient scalar profiles, and illustrate the creation of the large gradients seen in Fig. 10. For simplicity, only the solution for LEM wafers initially oriented either in the ${\mathrm{J}}_r$ or the ${\mathrm{K}}_r$ domain is shown. The upper left sub-plot shows the initial solution, which corresponds to the rightmost part of the black profile in the bottom plot of Fig. 10. Evolving from this, molecular diffusion smears out the steep gradient of the sharp step function prescribed by RANS, resulting in the smooth profile in the upper middle figure. Next, as seen in the upper right sub-plot, the radially oriented ${\mathrm{J}}_r$ or ${\mathrm{K}}_r$ domain within the control volume $r/\varDelta x \in [0.5,1.5]$ is rotated into the axial orientation (I-orientation) by the auxiliary coupling. With the wafers oriented in the axial direction, these are shifted downstream by the co-flow velocity and new fresh wafers are advected into the domain, as seen in the lower left sub-plot. In the lower middle sub-plot, the axially oriented profile is rotated back into either of the radially oriented ${\mathrm{J}}_r$ or ${\mathrm{K}}_r$ domains. Finally, the profile is affected by both molecular diffusion and reaction, resulting in local temperature maxima above and a local minimum below the initially prescribed profile, as seen in the lower right sub-plot.

The sketches of Fig. 11 are conceptual only, and the actual processes tend to deviate slightly from the above, e.g. by that the diffusion and chemical reactions are simultaneous and continuous processes on the small scales followed by the advection and instantaneous rotations associated with the larger advective time step $\varDelta t$. Also, a given profile is not necessarily rotated directly back to its original orientation but may go through several diffusion-reaction and advective time steps before this takes place. Nevertheless, the conceptual illustration of the development of the radially oriented profile depicted in Fig. 11 provides a plausible explanation of the emergence of the very large gradients observed in Fig. 10 and the consequent early mixing and reaction in the vitiated co-flow burner.

4.3 Discussion

Both the study of the freely propagating flame configuration in Sect. 4.1 and that of the vitiated co-flow burner in Sect. 4.2 show that the dimensional recoupling mechanisms cause large scalar gradients in the respective flow fields. In the case of the freely propagating flame, the stochastic rotations and the advective flipping of wafers create pockets of unburned wafers in the reaction zone, leading to locally increased consumption rates and stabilization of the flame front either at the upstream cell face or near the center of the control volume of initialization. The study of the laminar flame clearly shows that the rotational coupling plays a more dominant role than the advective coupling in the stabilization of the flame. With both the dimensional recoupling mechanisms switched off, however, the flame front stabilizes exactly at the same location as obtained from a LOGEresearch simulation. This serves as a validation of the diffusive-reactive part of the LEM3D code and provides strong indication that the diffusion and chemical kinetics algorithms are correctly implemented.

The investigation of the vitiated co-flow burner revealed that chemical reactions and heat release occur already in the upstream first layer of control volumes, just inside the co-flow inlet and in the neighboring control volumes of the central fuel jet. The auxiliary coupling in the form of stochastic rotations is shown to cause very large gradients in the radial temperature profiles, as illustrated in Fig. 10. Similar large gradients occur in the profiles of the chemical species, leading to enhanced molecular diffusion and the observed early mixing and reaction in the vitiated co-flow burner.

The model artifact of large gradients created by the rotational coupling has been discussed previously, and it has been proposed to use coarser 3DCVs to remedy this issue (Sannan et al. 2013); coarsening the 3DCVs will however results in larger gradients per rotation. The coarsening approach has the advantage of being computationally less expensive since the total number of LEM wafers is reduced for a given spatial resolution of the wafers. The size of the 3DCVs must be balanced by general model performance considerations, however, and the coarse-scale 3D-resolution must be fine enough so that mean-flow resolution requirements are fulfilled. In the case of the vitiated co-flow burner, we note that a coarser mesh does not seem beneficial due to the fact that the jet inlet is approximated by a single grid cell only and a fairly coarse mesh has already been chosen.

Another potential method to abate the coupling artifact of large gradients is to reduce the frequency of the rotations. Various iteration procedures, coupled with different values of the model constant $C_{\mathrm{rot}}$ of Eq. (4) to vary the rotational frequency, have been investigated in previous work (Grøvdal et al. 2020). That study showed that simulations with the auxiliary coupling switched off altogether, corresponding to $C_{\mathrm{rot}}=0$, result in inaccurate property profiles for the vitiated co-flow burner. It was also shown that a reduced rotational frequency did not seem to mitigate the issue of large gradients. In the current work we have used $C_{\mathrm{rot}}=1$, corresponding to that the 3DCVs on average are rotated 1.5 times during the advective residence time of the LEM wafers within the 3DCVs (Sannan et al. 2013). As a general remark, it should be noted that the auxiliary coupling of rotations is an important part of the LEM3D model to ensure that physical processes are consistently represented in all three spatial directions and cannot be switched off completely.

As demonstrated by the study of the vitiated co-flow burner, the near-field geometry of turbulent jets seems unsuitable for accurate modeling and simulation by LEM3D, notwithstanding the physically based representation of the fundamental processes. This suggests an alternative approach in which the LEM3D simulation domain is a sub-domain of a larger, more general flow domain. In this approach, a global flow solver (RANS or LES) is applied to the larger flow domain, while LEM3D simulations are applied in the smaller sub-domain for which high-resolution treatment of scalar mixing and reaction is of particular interest. This has the benefit of reduced computational cost, but also opens up for applications to more complex geometries than defined by the LEM3D cuboidal geometry. In the case of the vitiated co-flow burner, the sub-domain of interest may be defined as a region from the flame base and downstream to the flame tip. One procedure is to apply LEM3D as a post-processing tool, as is the case in the current formulation, but with mean-flow information provided to LEM3D limited to the sub-domain.

Another viable approach is to use LEM3D as a sub-grid scalar closure to a global flow solver. The LEM3D cuboidal geometry is constructed such that the 3DCVs naturally can be identified with corresponding, say, cubical LES control volumes. As a scalar closure for LES applications, LEM3D can provide small-scale resolution at the sub-grid scale in a given region of interest. This approach is broadly analogous to LEMLES in its physical treatment, but differs in its overall structure by providing orientation-dependent resolution at the sub-grid scale. One limitation to LEM3D in this regard is the restriction to a Cartesian mesh in the current implementation. Therefore, interpolation will be needed to transfer data between the LEM3D mesh and an arbitrary meshing structure, e.g. in the case that a curvilinear or and unstructured grid is used for the LES solver. Hence, a strategy for applying LEM3D to more complex geometries, such as for a gas turbine combustor or an industrial burner, is to run LES with LEM3D in a Cartesian subvolume, and to use LES with a conventional combustion model or the splicing method of LEMLES in the domain outside of that subvolume. This strategy will be well suited for cases where there is a need for high-fidelity resolution of scalar mixing and reactions in regions far from the combustor walls.

5 Conclusions

The dimensional-decomposition approach has been investigated in this paper, with a detailed description of the currently used recoupling mechanisms in LEM3D. The studies of the freely propagating flame configuration and the vitiated co-flow burner both have illuminated certain artifacts of the approach, such as the occurrence of very large gradients in the flow fields. This is mainly due to the auxiliary rotational coupling which brings dissimilar fluid states into contact and thus causes increased molecular diffusion and fuel consumption in the reaction zone.

In the case of the freely propagating flame, the auxiliary coupling leads to the stabilization of the flame front at the upstream cell face of the 3DCV in which the flame is initialized, independent of the size of the control volumes and the type of initialization. With the auxiliary coupling turned off, the advective flipping of wafers is the only dimensional recoupling mechanism of the model and the propagating flame stabilizes close to the center of the control volume of initialization. In this case, the flame front is slightly unstable due to the flipping of wafers. When both the recoupling mechanisms are turned off, the flame front stabilizes at the exact same location as obtained by the LOGEresearch software using the same chemical kinetics. This provides a validation of the diffusive-reactive part of the LEM3D code and gives further evidence that the demonstrated artifacts are not due to issues connected to the chemistry implementation of the code.

The investigation of the vitiated co-flow burner showed that chemical reactions and heat release take place already in the upstream first layer of control volumes. Due to the auxiliary rotations, large scalar gradients are created in the near field of the burner jet nozzle, leading to early mixing and reaction of the hydrogen fuel. This gives a stabilized flame with a low lift-off height in comparison to the experimental measurements.

The two flame studies presented demonstrate that the dimensional-recoupling mechanisms cause very large scalar gradients, leading to early mixing and locally increased consumption rates in the flow configurations. Since each 1D LEM domain represents species transport and chemical reaction in one direction only, a coupling of the domains is a necessity to ensure a physically consistent formulation in three dimensions (Sannan et al. 2013; Weydahl 2010). Thus, in order to limit model artifacts, a viable approach is to restrict LEM3D to a sub-region and rely on traditional mixing and chemistry solvers where creation of property discontinuities and other LEM limitations are significant. Therefore, in the present formulation of LEM3D, applications of the model should be restricted to sub-regions where high-resolution treatment of scalar mixing and reaction is of particular interest.

In the current formulation, LEM3D is used as a post-processing tool to RANS. In this approach, mean-flow information from RANS provides model input to LEM3D, which complements RANS with unsteadiness and scalar statistics needed for more accurate mixing and reaction calculations. A strategy for future work is to use LEM3D as a sub-grid scalar closure to LES. In applications to complex combustor geometries, LEM3D can then be used to solve fine-scale mixing and reactions in a Cartesian subvolume far from the combustor walls, while LES with conventional combustion closure is used in the region outside of the LEM3D sub-domain.

The capability of LEM3D to provide any type of computational cost savings for equivalent or enhanced accuracy for equivalent cost in comparison with established turbulent combustion (or turbulence-chemistry-interaction) models has not been attempted in the present work. Since only simple flame cases have been considered, it remains to be seen how accurate and robust LEM3D will be in realistic flows such as those in piston engines and gas turbines.