Mariovane Donini
email: mariovane.donini@inpe.br

Fernando Fachini
email: fernando.fachini@inpe.br

Cesar Cristaldo
email: cesarcristaldo@unipampa.edu.br

Pascal Bruel
email: pascal.bruel@univ-pau.fr

Some Specific Aspects Related To The Use Of The Artificial Compressibility Methods To Simulate Unsteady Flows

Keywords: Inert wake, buoyant diffusion flame, Tsuji configuration, finite volume method

The artificial compressibility method is applied to simulate non-reactive and reactive flows with incompressible and infiniterate chemistry assumptions. To validate the implemented time-accurate approach, the resulting code is used to calculate three different flow regimes: an unsteady non-reactive flow past a circular cylinder, a steady diffusion Tsuji flame and an unsteady buoyant diffusion flame. As expected, the results evidence the strong dependence of the convergence rate to the artificial compressibility factor. They also show the ability of this approach to describe the circular cylinder wake unsteadiness as well as the flame shape of the steady diffusion flame. In return, the correct prediction of the low-frequency instability known to be featured by unsteady buoyant diffusion flames is presently out of reach of the present version of the methodology.

Introduction

The artificial compressibility method is a well-established numerical approach introduced by Chorin [1] for solving the steady incompressible Navier-Stokes equations. In this method, the continuity equation is modified by the addition of a nonstationary pressure term along with an artificial compressibility factor. With this, waves of finite speed are introduced to distribute the static pressure throughout the whole computational domain. The results are physically meaningful only when a steady-state solution is reached and the original continuity equation is recovered [START_REF] Bruel | A pseudo-compressibility method for reactive flows at zero mach number[END_REF]. The ongoing popularity and success of the artificial compressibility method are mainly due to its simplicity and clear physical interpretation.

To obtain a time-accurate solution, a dual time-step artificial compressibility technique can be employed [START_REF] Soh | Unsteady solution of incompressible Navier-Stokes equations[END_REF][START_REF] Rogers | Numerical solution of the incompressible Navier-Stokes equations for steady-state and timedependent problems[END_REF][START_REF] Corvellec | A time-accurate scheme for the calculations of unsteady reactive flows at low Mach number[END_REF][START_REF] Dourado | A time-accurate pseudo-compressibility approach based on an unstructured hybrid finite volume technique applied to unsteady turbulent premixed flame propagation[END_REF]. The equations are iteratively solved such that the velocity approaches the new value in time as its divergence approaches zero. Then, for each physical time step, the pressure field has to go through one complete steady-state iteration cycle. To satisfy the divergence-free constraint on the velocity, the same classical artificial compressibility relation in the continuity equation remains unchanged, but the physical time is introduced in the momentum equation.

As far as reacting flows are concerned, the artificial compressibility method was extended so as to deal with steady (Bruel et al. [START_REF] Bruel | A pseudo-compressibility method for reactive flows at zero mach number[END_REF]) and unsteady turbulent premixed flames (Corvellec et al. [START_REF] Corvellec | A time-accurate scheme for the calculations of unsteady reactive flows at low Mach number[END_REF], Dourado et al. [START_REF] Dourado | A time-accurate pseudo-compressibility approach based on an unstructured hybrid finite volume technique applied to unsteady turbulent premixed flame propagation[END_REF]) as well as laminar confined and unconfined diffusion flames (Fathi et al. [START_REF] Fathi Azarkhavarani | Artificial compressibility method on half-staggered grid for laminar radiative diffusion flames in axisymmetric coordinates[END_REF], Bianchin et al. [START_REF] Bianchin | On the global structure and asymptotic stability of low-stretch diffusion flame: forced convection[END_REF]).

On the basis of the code developed by Bianchin et al. [START_REF] Bianchin | On the global structure and asymptotic stability of low-stretch diffusion flame: forced convection[END_REF] to investigate steady Tsuji diffusion flames, the purpose of the present study is to develop a time-accurate version and test it on three different flow configurations: an unsteady non-reactive flow past a circular cylinder in a channel (Case A), steady Tsuji diffusion flames (Case B), and an unsteady buoyant diffusion flame (Case C). To the best of the authors' knowledge, this is the first time the artificial compressibility method is applied to compute the whole buoyant Tsuji burner geometry (Case C). As a first step towards the full and more complex implementation of the strong density variations associated with reactive flows, the effects of the thermal expansion on the CSP 107-2 flow field are only taken into account, whenever relevant, to express the buoyant force present in Case C. In that case, an isobaric approximation for the equation of state is used to implement the temperature differential in the buoyant term, which is not restricted to small density differences. Since the flow is considered bi-dimensional and the effect of the thermal expansion on the flow field is not considered ( ˆˆ0    u ), where

 

ˆˆ, uv  u is the flow velocity vector, the sources o f vorticity ω in the vorticity equation are absent i.e. ˆ(

) 0    ωu, ˆ( ) 0    ωu , ˆˆ( ) 0     g and ˆˆ( ) 0 p     ,
where p is the pressure, ĝ is the gravitational acceleration and  is the density. Thus, it may be anticipated that the calculated vortices generated by the instability of the buoyant diffusion flame (Case C) will be simply convected to the flame tip and damped there by the viscous effects [START_REF] Xia | A vortex-dynamical scaling theory for flickering buoyant diffusion flames[END_REF]. In the previous expressions and in the following sections, any quantity  is dimensional whereas its dimensionless counterpart is written as . 

Numerical Treatment and the Artificial Compressibility Scheme

The introduction of a non-stationary term in the continuity equation leads to the following generic form

1 p        u (1)
where  is the pseudo-time and  is the artificial compressibility factor. The dimensionless momentum equations are given by 12 ()

t p Re             u u uu u ( 2
)
where t is the physical time, and Re the Reynolds number, t  and   are the physical and pseudo-time derivatives. The existence of the pseudo-wave propagation phenomenon associated with the solution of the above set of equations can be evidenced by writing them under a characteristics-like form (neglecting the viscous term and the physical time derivative and considering a one-dimensional configuration), namely

u 1 u 1 (u ) 0 (u ) (u ) pp c c x c x                             (3)
The artificial sound speed , c and the corresponding artificial Mach number , M are related to  by

2 2 uu u , 1 u cM c         (4)
By the above equation, it is clear that M is always less than 1 for all 0   . The discretization of the time derivatives in the conservation equations, Eqs. (1) and (2), is expressed as

1, 1 1, 1, , n n n ij p p RHS            (5) 1, 1 1, 1, 1, , n n n n n ij RHS t                   u u u u (6)
where superscripts n and  refer to the iteration cycle in physical and pseudo-time, respectively.   is the pseudo-time step, t  is the physical time-step and RHS refers to the right-hand-side of Eqs. (1) and [START_REF] Bruel | A pseudo-compressibility method for reactive flows at zero mach number[END_REF]. The flow chart of the algorithm to solve Eqs. ( 5) and ( 6) is presented in Figure 1. To advance the solution by one physical time-step, the equations are iteratively solved such that

n n n i j i j Res RHS t       u u u . A lower limit is given in CSP 107-3
the literature for the value of   , but numerical experiments are required to obtain an optimal value. The choice of the physical time step t  is specific to the problem under consideration (steady flow, forced excitation or large vortical structures), but is generally much lower than   . The rate of convergence of the solution in the pseudo-time integration depends on  as well as the streamwise length L , of the flow geometry. 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 distribute the pressure in order to converge to the steady-state solution in pseudo-time. Analyzing the pseudo-time needed for the pseudo-waves to travel downstream and back upstream, Chang and Kwak [START_REF] Chang | On the method of pseudo compressibility for numerically solving incompressible flows[END_REF] introduced the number N of computational time-steps required to reach a converged solution in pseudo-time as

1 2L N      (7)
In the present work, an explicit second-order Runge-Kutta Ralston's method was adopted for the pseudo-time time integration due to its simplicity, robustness and low computational cost. 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.

Fig. 1: Flow chart of the algorithm where max t is the maximum physical time and q is the vector of variables u and p .

Tests Cases

1.An Unsteady Flow around Circular Cylinder (Case A) 1. 1. Geometry

The two-dimensional incompressible unsteady laminar flow around a cylinder with a circular cross section of diameter d placed eccentrically in a channel of height ˆ4.1 hd  is considered. This configuration corresponds to one of those used by Schafer and Turek [START_REF] Scha ¨fer | Benchmark computations of laminar flow around a cylinder[END_REF] for a benchmark of different solution approaches for solving the incompressible Navier-Stokes equations. As shown in Fig. 2c., the distances between the cylinder center and the bottom and top walls are ) and the normal component is set to 0 . Finally, a zero pressure gradient is imposed normal to all boundaries.

1. Results

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 Fig. 2c, 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 2a presents the time evolution of the non-dimensional normal component of the velocity at ( , ) ( 0.5,0.5)

xy  observed when the periodic regime is established. For 100 Re 

, the experimentally obtained Strouhal number is 0.287 0.003 St  [START_REF] Scha ¨fer | Benchmark computations of laminar flow around a cylinder[END_REF]. The power spectrum of the fluctuations of the streamwise component of the velocity is shown in Fig. 2b. The numerically computed Strouhal number of St = 0.289 agrees well (error of 1.35%) with the experimentally obtained value. For the present case, Figure 3a shows, for a given physical time-step, an evolution of the maxima of the residuals of the RHS (Eq. ( 5)) during the pseudo-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 pseudo-time steps necessary to achieve convergence to the new physical time step after the snapshot (Fig. 2c) is illustrated in Fig. 3b where L is defined as the channel length and 0.14 in the literature [START_REF] Chang | On the method of pseudo compressibility for numerically solving incompressible flows[END_REF][START_REF] Bruel | A pseudo-compressibility method for reactive flows at zero mach number[END_REF]. By using Eq. ( 7), the dashed line in Fig. 3b represents the number of minimum time-steps to achieve convergence of the pseudo-time integration for each value of .

 Considering the convergence criteria of max RHS   , it is possible to see a good agreement between the computation iterations and the iterations described by Eq. ( 7) until

2. Steady (Case B) and Unsteady (Case C) Diffusion Flames 2. 1 Geometry

Two different diffusion flames regimes are considered.  The first regime (Case B) is concerned with the flame stabilization over a cylindrical porous burner with radius r inside a channel. The gaseous fuel is injected from the forward half part of the burner with velocity b û into the incoming airflow of velocity air û . 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 [START_REF] Tsuji | The counterflow diffusion flame in the forward stagnation region of a porous cylinder[END_REF]. In such a configuration characterized by a low incoming flow velocity, an envelope steady flame is found.  The second diffusion flame test case (Case C) introduces some modifications to the Tsuji burner configuration, i.e., fuel injection from the whole burner surface, buoyant force, no incoming airflow and no channel walls. These modifications enable the appearance of a long buoyant diffusion flame in the open air which exhibits a self-excited periodic movement that is often referred to as 'puffing'. The low-frequency flame oscillation is used as a parameter to test the ability of the method to describe the flame unsteadiness [START_REF] Xia | A vortex-dynamical scaling theory for flickering buoyant diffusion flames[END_REF].

2. 2. Governing Equations

Thermodynamic and transport coefficients, i.e. the specific heat at constant pressure ˆp c  , the viscosity  , the thermal conductivity k and the diffusion coefficient D , are considered constant. In this analysis, both fuel and oxidant Lewis numbers are equal to unity. Also, an infinitely fast one-step chemical reaction is assumed i.e.

F

O (1 ) P ss    (At stoichiometry, s mass of oxygen is consumed for each unit mass of fuel 𝐹 resulting in 1 + 𝑠 mass of products 𝑃). Thus, the temperature T and the mass fraction

fields i Y ( , i O F 
) are determined by two conserved scalars: the mixture fraction Z and the excess enthalpy H [START_REF] Lin | Fundamental aspects of combustion[END_REF][START_REF] Fachini | Theory of flame histories in droplet combustion at small stoichiometric fuel-air ratios[END_REF][START_REF] Fachini | Extended Shvab-Zel'dovich formulation for multicomponent-fuel diffusion flames[END_REF]. The continuity equation is defined as Eq. ( 1) and the dimensionless conservation equations are

2 t           E E u E C E F (
H S T Q Y Y     , in which ˆ/ b FO S sY Y   and ˆ/ b Fp Q Y Q c T   
, where ˆb F Y and ˆO Y  are the mass fraction of fuel at the burner injection surface and oxidant at the ambient, respectively, and 𝑄 is the heat of combustion. The dimensionless variables are defined as

b ˆˆ(u / ) t t r  , ˆ/ x x r  , ˆ/ y y r  , b ˆû u/ u  , ˆv v/ u b  , ˆ/ T T T   and 2 b /u pp    .
The oxidant mass fraction is normalized according to

ˆ/ O O O Y Y Y  
and the fuel mass fraction is defined as

ˆ/ b F F F Y Y Y 
. Equations ( 1) and ( 8) are integrated with the following boundary conditions:

On the symmetry axis ( 0)

x  , x u0 x x x x v p Z H           ; at the burner boundary surface 2 2 2 ) (1 r x y      , u v 0 n x y p       , 1 s n Z Z Pe Z S     , 1 s n H H Pe H S     (
    and ( 1) / O H S T Q Y     . At outlet ( 13) y  , they are u v 0 y y y y y p Z H         
 and at the channel wall ( 4)

x  , they read u v 0

x x x p Z H         . For

Results

The steady diffusion flame (Case B) results are presented in Fig. 4 for different values of fuel-ejection rate fw  and air û , in which . The continuous line and its corresponding circles are the numerical result of the current study, the dashed line and its corresponding squares are the numerical results of [START_REF] Tsa | Flame stabilization over a Tsuji burner by four-step chemical reaction[END_REF], and the dash-dot line and its corresponding triangles are the experimental measurements of [START_REF] Dreier | Comparison of CARS measurements and calculations of the structure of laminar methane-air counterflow diffusion flames[END_REF].

Fig. 4a compares the predictions obtained in this study to the numerical finite-rate chemistry and experimental results of Tsa and Chen [START_REF] Tsa | Flame stabilization over a Tsuji burner by four-step chemical reaction[END_REF] and Dreier et al. [START_REF] Dreier | Comparison of CARS measurements and calculations of the structure of laminar methane-air counterflow diffusion flames[END_REF], respectively. The presented infinite-rate combustion model reproduces the data CSP 107-7 measured in both numerical and 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 4b directly compares the flame-sheet obtained in this study with the flame boundary computed from fuel reaction-rate contours by Tsa and Chen [START_REF] Tsa | Flame stabilization over a Tsuji burner by four-step chemical reaction[END_REF]. The flame-sheet shape obtained 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. The simulation of Case C proves to be the most challenging of all. Indeed, such a low Froude diffusion flame is known to feature a self-sustained cyclic buoyancy driven puffing regime. The results obtained are illustrated by Figure 5 which displays the mass fraction isocontours and vorticity contours. Starting from the initial conditions, the flow field is experiencing:  A transient phase ( 0 30 t  ) characterized by the generation of a finite number of vortical tubular structures which are growing while being convected upwards (See in Fig. 5b to Fig. 5j, the bulged and rolled contour regions). These vortices are generated by buoyancy-driven shear flow, close to the burner, when the flame is progressively adapting its shape from the initial condition. The instability stops as soon as the last generated vortex structures leaves the computational domain (Fig. 5j).

 A (non-physical) steady state ( 30 t  ) characterized by the absence of any puffing and related vortical structures.

Here, it can be guessed that due to the absence of density gradient related vorticity sources (constant density hypothesis) and vortex stretching mechanism (two-dimensional flow), the dynamics of the flow leading to the puffing regime cannot be sustained beyond the observed transient period.

Conclusion

This work investigates the implementation of the unsteady artificial compressibility approach to simulate non-reacting and reacting flow fields in the limit of a zero Mach number. The resulting time-accurate scheme has been tested in three cases: unsteady non-reactive flow past a circular cylinder, steady Tsuji diffusion flame, and unsteady buoyant diffusion flame. The flow past a circular cylinder case has been 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 has been carried out. 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 pseudo-time agrees very well with the expression presented by Chang CSP 107-8

and Kwak [START_REF] Chang | On the method of pseudo compressibility for numerically solving incompressible flows[END_REF]. A comparison with experimental and numerical results was carried out for the steady diffusion flame case. The simulated temperature profile and flame shapes are in good agreement with those reported in the literature, except in the region around the maximum temperature, where the reaction layer is not well described certainly because of the infinite-rate chemistry considered. The third case investigated the influence of the hypothesis of constant-density flow in a buoyant diffusion flame. The results showed that the low-frequency instability generated by the displacement of the flame's initial condition was damped and led ultimately to a non-physical stationary diffusion flame. This is tentatively attributed to the absence of vorticity source terms related to the underlying constant density assumption. Hence, the full implementation of variable density and finite-rate chemistry will be the natural next step to be addressed in order to refine the description of unsteady buoyant diffusion flames through an artificial compressibility based method of solution.

Fig. 2 :

 2 Fig. 2: Case A (Mesh size of 43 246  ): (a) Velocity history at ( , ) ( 0.5, 0.5) xy  , (b) power spectrum for the vortex-shedding, and (c) snapshot of isocontour of velocity (× marker indicates the location at which the data in (a) is computed.

  opt  , value beyond which the convergence rate begins to degrade gradually until the value 38   beyond which convergence is lost. CSP 107-5

Fig. 3 :

 3 Fig. 3: Case A (Mesh size of 43 246  ): (a) Convergence history and (b) the influence of the artificial compressibility factor β on the number of time steps required to obtain a pseudo-time converged solution.

  to the definition of the mixture fraction function Z , the flame position ( , ) ff xy is given by the isoline ( , ) 1 Z x y  where the flame temperature f T is determined by ( , ) ( 1) /

  depicts the temperature profile along the forward stagnation streamline.

Fig. 4 :

 4 Fig. 4: Case B (Mesh size of 82 446  ): (a) Temperature distribution through the flame front of aTsuji burner with 38 Re  , 2.6 Fr  , ˆ0.02 rm  , air û 1.15 / ms  , and 0.318 fw  and (b) streamlines and flame shape with 18 Re  ,

Fig. 5 :

 5 Fig. 5: Case C (Mesh size of 82 446  ): Mass fraction contours at several times with 40 Re  and 0.8 Fr  .( ) a

  𝛼 ̂ thermal diffusivity, 𝑘 ̂∞ thermal conductivity and ĝ the gravitational acceleration), respectively, based on the burner properties: the fuel injection velocity b û and the cylinder radius r .

	The mixture fraction is expressed as	1   and the excess enthalpy Z S Y Y 	( 1) /
																FO				FO
																				8)
	in which	u		(u, v)	,	E		(u, v, , ) ZH	,	C		(	1111 , , , Re Re Pe Pe 	)	and	F	(     , p p F	,0,0)	. For Case C, the buoyant
																				x	y	g
	force																	g F		T Fr	, in which	Fr		b ˆˆ( u /	r g	)	is the Froude
	number and the equation of state	 	1/T		is assumed, where 𝜌 is the density. In Eq. (8), the parameters in	C are the Reynolds
	and Péclet numbers defined as	Re r 	b û/ 	and	Pe r 	b	û/ 	(with ˆ/     	,	ˆ/ p kc    

g

F is considered to be dependent on the temperature as 2 (1 1 / ) /



, where 𝜇̂ is the dynamic CSP 107-6 viscosity,

  Robin's like boundary type for Z and H ) where are the Z and H fluxes which are imposed at the burner injection surface 2 are found as part of the solution of the problem and that for Case B, this holds only in the forward part of the cylinder. For Case B the boundary conditions at the inlet

	1  , Z S Y	H	( 1) / S T Q Y   	and the subscript n stands for the normal to the burner surface. The terms
		sF			s	s	F
		s				s
						r		2 x y   as function of ˆb F 2 1  Y , ˆb T
	and ˆb u , namely	1 SS  and	S	( 1) / S T Q  	1  . Note that
					Z	Hb
	(	7.5) y 	are v1  , u0 y pZ

Z S and H S s F Y and s T

1. 2. Governing EquationsAn incompressible Newtonian fluid is considered for which the conservation equations of mass (artificial compressibility