An Artificial Compressibility Based Approach to Simulate Inert and Reacting Flows

An efficient methodology to simulate non-reactive and reactive flows is presented. Combining a finite-volume approach on fully staggered meshes along with the artificial compressibility method, the resulting code proves to be versatile enough to cope with flow configurations ranging from unsteady cylinder wakes, heated cylinder or steady and unsteady diffusion flames with excellent accuracy, in the limits of the underlying physical modelling.


Introduction
The artificial compressibility (AC) method is quite a well-established numerical approach for solving the constant-density incompressible Navier-Stokes equations. Originally introduced by Chorin [1] for steady flows simulations, this method consists in that case in modifying the continuity equation by adding a nonstationary pressure term, namely: where ̂ is a an artificial e.g. non-physical time and ̂ is the artificial compressibility factor which scales as a velocity square. In the following and unless stated otherwise, any quantity ̂ is dimensional whereas its dimensionless counterpart is written as .
The results are physically meaningful only when a steady-state solution in is reached so that the original continuity equation is recovered [2]. The ongoing popularity and relative success of the artificial compressibility method to deal with constant density flow simulations are mainly due to its simplicity and clear physical interpretation. To obtain a time-accurate solution, a dual time-step technique can be employed [3,4,5,6]. Accordingly, in addition to the aforementioned modification of the continuity equation, the physical time derivative of the velocity field is introduced in the momentum equation, namely (dimensionless form): where is the physical time, the Reynolds number, and are the physical and artificial-time derivatives, respectively. These equations are iteratively solved such that the velocity field approaches its new value in physical time as its divergence is driven towards zero. Thus, for each physical time step, the flow field has to go through one complete sub-iteration cycle in artificialtime.
The existence of the artificial-waves propagation phenomenon associated with the convergence towards the physically meaningful solution of the above set of equations can be evidenced by writing the momentum equation under a characteristics-like form (neglecting the viscous term and the physical time derivative and considering a one-dimensional configuration), namely: [ + ( ± ) −1 ] The artificial sound speed , and the corresponding artificial Mach number , are related to by: Thus, artificial waves of finite speed are introduced to distribute the static pressure throughout the whole computational domain. The rate of convergence of the solution during the pseudo-time integration loop heavily depends on the value of , and this can be thought of as a weak point of the method. Indeed, in order to converge to the steady-state solution during the course of each sub-iterations cycle in artificial-time, the waves associated to the hyperbolic nature of the artificial compressibility based system of equations have to undergo at least a one round-trip propagation to ensure the proper built-up of the pressure field (more precisely, of its gradient field) over the whole computational domain. Based on this representation, Chang and Kwak [6] estimated the number of artificial time-steps required to reach convergence in artificial time. By considering a characteristic length of the computational domain over which the artificial waves must travel once forth and back, they obtained a lower bound for expressed as: Compared to the numerous AC based simulations of constant density flows that can be found in the literature (quite a significant number of them have been recently recalled by Hodges [8]), much less examples of AC based simulations of non-constant density flows have been reported or are discussed in CFD textbooks, to the noticeable exception of Oran and Boris [9]. In such cases featuring a non-constant density field (in space and/or in time), one can distinguish between configurations still featuring a divergence free velocity field (Bassi et al. [10], Shapiro and Drikakis [11]) from those which did not. For the latter, they are mostly related to the simulation of Mach zero reacting flows such as steady turbulent premixed flames (Bruel et al. [2]), unsteady turbulent premixed flames (Corvellec et al. [5], Dourado et al. [7]), laminar confined and unconfined diffusion flames (Fathi et al. [12], Bianchin et al. [13], Donini et al. [14]).
Starting from the observation that the artificial compressibility is quite underused to simulate in particular reacting flows, the objective of the present work is to draw the attention of the CFD community on this valuable approach by showing how the combination of rather standard discretization tools in a finite volume framework allows to obtain an efficient and versatile simulation tool. In addition to the Introduction and Conclusion sections, the paper is organized in three main sections: Section 2 describes the different continuous systems of equations to be solved; Section 3 gives the main ingredients of the methods of solution whilst in Section 4, the results obtained for five different flow configurations are presented.

The Continuous Systems of Equations
Both inert and reacting gaseous flows are to be considered in the present study and the possible influence of the buoyancy force is accounted for. Thus, for each of the two classes of configurations, the modified (e.g. including the AC related terms) continuous system of equations dealt with is presented hereafter. 2D flow geometries are considered in this study, either in Cartesian or axisymmetric coordinates.

1. Inert Flows
A Newtonian fluid is considered for which the modified conservation equations of mass, momentum and energy are given by: Where =( , v) designates the dimensionless velocity vector, ̿ = μ(∇ + ∇ ) − 2 3 μ ( ⋅ ) is the viscous stress tensor and the power-law expressions μ = ρ = σ (with = 0.7) are employed to account for the temperature dependence of the shear viscosity μ and thermal diffusivity . Vector stands for the base vector aligned with the gravity acceleration vector ̂. The dimensionless variables and differential operator are defined as: the thermal conductivity and ̂ the reference specific heat at constant pressure) will be indicated later on, on a case by case basis. The above system is supplemented by the Mach zero equation of state, namely = 1/ . When the flow configuration dealt with is isothermal, the energy equation simply reduces to = .

Reacting Flows
Two different diffusion flames in ambient atmosphere are simulated to illustrate the capability of the AC method to deal with reacting flows with a strongly varying density. They are both described through a simple thermochemical model. It assumes an infinitely fast onestep chemical reaction i.e. + 2 → (1 + ) (At stoichiometry, s mass of oxygen is consumed for each unit mass of fuel F resulting in 1 + mass of products P). Thus, the temperature and the mass fraction fields ( = 2 , ) are determined by two conserved scalars [19,20,21]: the mixture fraction ≡ − 2 + 1 and the excess enthalpy ≡ ( + 1) / + + 2 in which ≡̂/̂2 and ≡̂̂/̂̂ where ̂ and ̂2 are the reference mass fraction of fuel and mass fraction of oxidant, respectively, which supplement the reference quantities already introduced in the previous subsection for inert flows. is the heat of combustion. Both fuel and oxidant Lewis numbers are taken equal to unity. Thus, the generic system of governing equations that describes the evolution of the reacting flows at hand is given by: The above system is supplemented by the Mach zero equation of state, namely = 1/ .

Method of Solution
The numerical solution of the Mach zero system of equations relies on a dual-step time-accurate artificial compressibility method. An explicit second-order Runge-Kutta Ralston's method was adopted for the artificial-time integration and the second-order Euler method was selected for the physical-time integration due to its simplicity of implementation in a dual-time step approach.
By replacing the derivative terms with their numerical approximations, the resulting set of equations can be written in compact form as: ( ) , is the right-hand side of the discretized equation. Then, introducing the residual value for the artificial time step as: permits to rewrite Equation 13 at times' step ( + 1, + 1) as: The integration steps in artificial-time are finally given by: where ( 1 , 2 , 3 ) = (1,3/4,1/3), ( 1 , 2 , 3 ) = (1,1/4,1/3) and is the cell volume. To advance the solution by one physical time-step, the equations are iteratively solved in a segregated way such that +1, +1 approaches the new value +1 as the artificial time derivative approaches zero. For satisfying this constraint, the residual value ( ) , of Equation 15 is set to reach values below ε = 10 −5 . The flow chart of the algorithm to solve the system is presented in Figure 1. A finite volume framework is adopted. The governing equations are discretized on a staggered grid by a cell-vertex finite-volume formulation using the quadratic upstream interpolation for convective kinetics (QUICK) scheme to guarantee stability, sensitivity to the flow direction and third-order truncation error. There is no need to prescribe boundary conditions for the pressure, thanks to the recourse of the fully staggered grid, as shown in Figure 2. The non-uniform structured grid is refined (hyperbolic tangent distribution), on a case by case basis, in zones where high gradients are expected.

Numerical Experiments
A total of five flow configurations were simulated broken down into i) three inert flow configurations (two isothermal and one featuring a strongly space varying density) and ii) two reacting flow configurations (one steady laminar diffusion flame and one flickering laminar diffusion flame).

1. Inert Flows
Three different inert flow configurations are considered. The first one (the oscillating plate) is a verification test aimed at assessing the effective accuracy of the method of solution. The two others are validation tests that demonstrate the capacity of the method to cope with flow featuring unsteady vortices (unsteady cylinder wake) or strong density variation (a heated cylinder placed in a square enclosure).

1. 1. The Oscillating Plate (Stokes' Second Problem)
The flow motion over an infinite flat plate that oscillates parallel to itself is investigated. The coordinate system is 2D Cartesian, with ̂ being the coordinate along the plate and ̂ the coordinate normal to it. The plate oscillates in the ̂= 0 plane with a velocity given by: ̂ designates the plate oscillation frequency and ̂ the time. The fluid is air and the reference kinematic viscosity is taken at ambient temperature e.g. ̂= 1.55 × 10 −3 2 / . The maximum plate velocity is taken as the reference velocity e.g. ̂≡̂= 2 × 10 −2 / and the reference length is chosen equal to eight times the depth of penetration of the viscous wave e.g. ̂= 8 × 6.25 × 10 −3 = 5 × 10 −2 . Thus, the Reynolds number is such that = 64. The fluid velocity is given by (Schlichting [18]): where ̂= √̂ . The global error at time level will be denoted by = ℚ − , where ℚ , is the computed value at each point of the grid and , represents the exact value.
To quantify the error, the commonly used p-norm is chosen in order to estimate the error at a given physical time level , namely: Let's assume the method has an order of accuracy , then the error is expected to behave like ‖ ‖ = (Δ , Δ ) + high-order-terms, as the time step or grid spacing are decreased and (Δ , Δ ) → 0. Figure 3 shows the results of the time step and mesh refinement influence study for the Stokes' second problem using the p-norm with = 1. The exact and computed solutions are compared on a sequence of time steps and grid spacing, and the norm of the error is plotted as a function of Δ and Δ . These are shown on a log-log scale, e.g. log‖ ‖ 1 ≈ log| | + log|(Δ , Δ )| so that a linear behaviour is expected in this plot, with the slope providing the effective order of accuracy . When decreasing Δ , the mesh refinement is adapted in order to keep constant the Courant-Friedrichs-Lewy number at the value = 0.5

1. 2. The Unsteady Wake of a Flow past a Circular Cylinder
The incompressible unsteady laminar flow around a cylinder with a circular cross section of diameter ̂≡̂ placed eccentrically in a channel of height ĥ = 4.1 ̂ is considered (see Figure 4). The coordinate system is 2D Cartesian, with ̂ being the streamwise coordinate and ̂ the coordinate normal to the channel walls. This configuration corresponds to one of those used by Schäfer and Turek [16] for a benchmark of different solution approaches for solving the incompressible Navier-Stokes equations. The distances between the cylinder centre and the bottom and top walls are 2.1 ̂ and 2̂, respectively. The Reynolds number is defined by =̂̂/̂, where ̂≡̂ is the kinematic viscosity and v ≡̂ denotes the bulk velocity. The case selected here corresponds to = 100. As illustrated in Figure 4  Among all the characteristics of this type of problem (lift, drag, and pressure coefficients), the correct prediction of the periodic vortex-shedding, illustrated by the isocontour of velocity in Figure 5c was the target chosen to validate the present solution approach. In particular, the Strouhal number is computed to measure the ability of the method to produce quantitatively accurate unsteady results. Figure 5a Figure 5b. The numerically obtained Strouhal number value is St = 0.289 which agrees well (relative error of 1.35%) with its experimentally obtained counterpart. For the present case, Figure 6 shows, for a given physical time-step, an evolution of the maxima of the residuals ( ) , (Equation 15) during the artificial-time iteration cycle. It can be seen that calculations for = 40 become unstable and within 40 steps starts to diverge, whereas other cases converge to a stable solution. For low values of , it takes only 10 iterations for the solution to reach almost constant residual values but then, further iterations do not improve the solution and so the accuracy constraint cannot be met. The effect of the artificial compressibility factor on the number of artificial-time steps necessary to achieve convergence to the physical time step following that corresponding to the snapshot of Figure 5c is illustrated in Figure 7. is defined as the channel length and = 0.14. The optimum value of ≈ 34 is higher than the expected value of ≈ 8 reported in the literature [2,11]. By using Equation 5, the dashed line in Figure 7 represents the number of minimum time-steps to achieve convergence of the artificial-time integration for each value of . Considering the convergence criteria of ( ) < , it is possible to see a good agreement between the computation iterations and the iterations described by Equation 5 until .  For > , the convergence rate begins to degrade gradually until ≈ 38, above which convergence is lost.

1. 3. A Square Enclosure with a Heated Circular Cylinder
The configuration is 2D Cartesian. It consists of a cooled square enclosure kept at a constant temperature ̂, with sides of length ̂, within which a heated circular cylinder at constant temperature ̂ℎ , with a radius ̂( = 0.2̂), is located in the center of the enclosure. This configuration was studied by Kim et al. [22] to examine the natural convection phenomena by changing the location of the circular cylinder. The case selected here corresponds to = 0.2, = 6.7 and = 4.9 so that the Rayleigh number defined as = (̂̂̂3)/(̂̂) with ̂= 2(̂ℎ −̂)/(̂ℎ +̂) is such that = 10 3 , leading to reference values ̂=̂= 1.4 × 10 −2 and ̂= 7.49 × 10 −2 / . The hot cylinder and cold wall temperatures are such that ̂ℎ = 1200 and ̂= 300 , respectively.
Once the velocity and temperature fields are obtained, the local Nusselt number and the surfaceaveraged Nusselt number are defined as: Where = (̂−̂)/(̂ℎ −̂), is the normal direction with respect to the walls and is the wall surface area over which the integration is performed. The Nusselt numbers results presented by Kim et al. [22] were used to evaluate the efficacy of the present solution approach. Figure 8 present the predicted flow topology streamlines. This topology highlights the buoyant induced flow, caused by the upward motion of the hot fluid that moves along the surface of the heated cylinder. The hot fluid then encounters the cold top wall becoming gradually colder and denser while it moves horizontally outward of the vertical center line. At this stage, the fluid is cold enough to descend along the cold side walls closing the recirculation zone. Figure 9 displays the distribution of the local Nusselt number along half of the wall of the enclosure. Because the problem presented symmetry about the vertical center line at = 0, only the right half of the enclosure is considered (see Figure 9) for calculating . The present simulation yields a value = 1.73 which is within 3.59 %. of that reported by Kim et al. [22] e.g. = 1.67.

Reacting Flows 4.2.1 Tsuji Burner Flame
This configuration is concerned with the flame stabilization over a porous cylindrical burner with radius ̂≡̂ inside a channel. The geometry is 2D Cartesian. As shown in Figure 10, the gaseous fuel is injected from the forward half part of the burner with velocity ̂≡ into the incoming airflow of velocity ̂. This configuration reproduces the experimental set-up of the Tsuji Burner, where the rear side of the burner surface was coated to avoid the ejection of fuel into the wake region [17]. In such a configuration characterized by a low incoming flow velocity, an envelope steady flame is found. The flame is described by the set of Equations 9-12 where the other reference quantities are chosen as ̂≡̂, ̂2 ≡̂2 ∞ , ̂≡̂∞, ̂≡̂∞ and ̂≡̂∞ where the index and ∞ denote quantities taken at the burner exit and in the ambient atmosphere, respectively.  Figure 11 for different values of fuelejection rate − and ̂, in which − = (̂/ )( /2) 0.5 . This figure depicts the temperature profile along the axis of symmetry at the forward part of the cylinder. Figure 11a compares the predictions obtained in this study to the numerical finite-rate chemistry and experimental results of Tsa and Chen [23] and Dreier et al. [24], respectively.  [23], and the dash-dot line and its corresponding triangles are the experimental measurements of [24].
The presented infinite-rate combustion model reproduces the data reported in both numerical and Symmetry condition experimental studies, except in a small, but important, region around the maximum temperature. The profiles show that, for this study, the maximum temperature is approximately 2200 K (adiabatic flame temperature for methane) with a sharp temperature profile, while it is about 1900 K in the experimental study with a rounded distribution. This is due to the limitation of infinite-rate chemistry to describe the coexistence of reactants in the reaction layer that is approximated as a flame-sheet, i.e., the reactants must reach the flame in stoichiometric proportions. Figure 11b directly compares the flamesheet obtained in this study with the flame boundary computed from fuel reaction-rate contours by Tsa and Chen [23], represented as the dashed-line. The flamesheet shape obtained (solid-line) is similar to that given by the reaction-rate contours of the finite-rate computation, except in the wake distant from the cylindrical burner, at which the recirculation zone is affected by the thermal expansion.

A Flickering Diffusion Flame
The unconfined flickering jet diffusion flame case was chosen to validate the full implementation of the present numerical code. The geometry of this case is 2D axisymmetric. As shown in Figure 12, the burner is composed by a fuel jet with a radius ̂≡̂= 1.3 × 10 −2 surrounded by a annular air stream with radius ̂= 13 × 10 −2 . The fuel is methane diluted by 50%,  The resulting times steps were Δ̂= 4.3 × 10 −3 and Δ̂= 3.9 × 10 −4 , for physical and artificial times, respectively. The value of β = 8 led to a suitable convergence rate for this set of time steps. The case introduced above reproduces the computations done by Davis et al. [28] who also used the flame sheet model, unit Lewis number hypothesis, and validated their results against the experimental data by Chen et al. [26]. One strong motivation behind the choice of this case was to assess the ability of the present code to describe the temporal behavior and the formation of vortical structures due to large density gradients and buoyancy effects. The low frequencies of flame oscillations (flame flickering) were in the range 5 − 15 and independent of the fuel type, the geometry of the source of fuel and the flow field in the wake [25]. The coupling within the flow field between the accelerations around the flame and decelerations in the plume above it due to the buoyant force dramatically impacts the temperature and species field dynamics and is at the origin of the formation of large vortices outside of the flame. As a vortex ascends along the flame in direction of the tip, it is forced against the flame. Close to the flame tip, the vortex strangles the flame, a bottleneck appears featuring a large strain rate which leads to the local extinction of the flame ending up in the separation of part of the flame tip (fuel pockets) which is carried away by the flow [26].The frequency analysis of this behaviour of the flickering jet diffusion flame is shown in Figure 13. agrees well (relative error of 5.5%) with the value predicted by the correlation − , ̂= 0.29 −1 = 7.96 , suggested by Sato et al. [27]. Since the large-scale instability is produced mainly by buoyancy, its frequency is an increasing function of the buoyancy strength as observed when plotting the Strouhal number as a function of the Froude number. The weaker secondary frequency peak visible in Figure  13b probably results from the preceding flame bulge interaction with the trailing one and is categorized here as a sub-harmonic. The flame evolution between ̂= 1.3 and ̂= 1.52 is illustrated by Figure 14 which displays the temperatures isocontours and vorticity contours where the dimensional vorticity is defined by: The bluish regions indicate the clockwise vortex structures and the reddish regions the counter clockwise vortex ones. The flame bulge is prominent by the isocontour of temperature. Meanwhile, the red isoline represents the flame sheet and evidences the detachment of fuel pockets. As previously described, these pockets are regions of hot gas enclosed by a flame that travels upward. The effective entrainment of oxidizer in the region of the axis of symmetry is found to be the main mechanism of the flame local extinction and release of this secondary flame. Detached fuel pockets can be observed in Figures 14c and 14f. Also, the Figure  14 displays the probe location where the data from Figure 13 were acquired. This location was chosen for illustrating the passage of the large vortical structure formed by the buoyant effect, a key feature for the correct flame flickering frequency analysis. Beyond the correct prediction of the flickering frequency, the unsteady flame structure appeared to be in line with the results presented in the literature (not shown here).

Conclusions
This work investigates the implementation of the unsteady artificial compressibility approach to simulate non-reacting and reacting flow fields in a zero Mach number framework. The resulting time-accurate scheme was tested in five different situations cases: the Stokes' second problem, the unsteady wake flow past a circular cylinder, the flow around a heated cylinder placed in a square enclosure, a steady Tsuji diffusion flame, and a flickering buoyant diffusion flame. The flow over the oscillating plate and past a circular cylinder case were chosen to put into evidence the basic properties of the implemented time-accurate approach, such as the ability to describe unsteady non-reactive flows and its convergence rate. An analysis of the influence of the artificial compressibility factor on the convergence rate was carried out. As expected, the results showed that the stability of the numerical code is highly dependent on the value of the artificial compressibility factor and the computed number of time-steps required to reach convergence in artificial-time agreed very well with the expression presented by Chang and Kwak [6]. For the flow in the enclosure hosting a heated cylinder, a configuration that featured quite high density gradients, the predicted flow topology as well as the heat flux distribution along the enclosure walls matched very well those reported in the literature. Then, a comparison with experimental and numerical results was performed for a steady diffusion flame case. The simulated temperature profile and the flame shape were in good agreement with those reported in the literature, except in the region around the maximum temperature, where the reaction layer was not so well described, certainly because of the infinite-rate chemistry assumption considered here. The fifth case investigated the capability of the present approach to predict correctly the large scale unsteadiness of a flickering diffusion flame. An excellent agreement with the experimental results was observed in term of the fundamental frequency of the flame flickering. Qualitatively speaking, the predicted flame topology appeared to be in line with that observed experimentally.