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

An Artificial Compressibility Based Approach to Simulate Inert and Reacting Flows

Keywords: Mach zero framework, finite volume, fully staggered mesh, steady flows, unsteady flows, diffusion flames

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 [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] for steady flows simulations, this method consists in that case in modifying the continuity equation by adding a nonstationary pressure term, namely:

𝛽 ̂-1 𝜕 𝜏 ̂𝑝̂= -𝛻 ⋅ 𝒖 ̂ ( 1
)
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 [START_REF] Bruel | A pseudocompressibility method for reactive flows at zero Mach number[END_REF]. 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 [START_REF] Soh | Unsteady solution of incompressible Navier-Stokes equations[END_REF][START_REF] Rogers | Numerical solution of the incompressible Navier-Stokes equations for steadystate and time-dependent problems[END_REF][START_REF] Corvellec | A timeaccurate scheme for the calculations of unsteady reactive flows at low Mach number[END_REF][START_REF] Chang | On the method of pseudocompressibility for numerically solving incompressible flows[END_REF]. 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):

𝜕 𝑡 𝒖 + 𝜕 𝜏 𝒖 = -𝛻 ⋅ (𝒖𝒖) -𝛻𝑝 + 𝑅𝑒 -1 𝛻 2 𝒖 ( 2
)
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 𝜕 𝜏 𝑝] +(𝑢 ± 𝑐)[𝜕 𝑥 𝑢 + (𝑢 ± 𝑐) -1 𝜕 𝑥 𝑝] = 0

(3)

The artificial sound speed 𝑐, and the corresponding artificial Mach number 𝑀, are related to 𝛽 by:

𝑐 = √𝑢 2 + 𝛽, 𝑀 = 𝑢 𝑐 = 𝑢 √𝑢 2 + 𝛽 < 1 (4)
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 [START_REF] Chang | On the method of pseudocompressibility for numerically solving incompressible flows[END_REF] 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:

𝑁 > √1 + 𝛽 𝛽 2𝐿 𝛥𝜏 (5)
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 [START_REF] Hodges | An artificial compressibility method for 1D simulation of open-channel and pressurized pipe flow[END_REF]), 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 [START_REF] Oran | Numerical simulation of reactive flow[END_REF]. 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. [START_REF] Bassi | Artificial compressibility Godunov fluxes for variable density incompressible flows[END_REF], Shapiro and Drikakis [START_REF] Shapiro | Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part I. Derivation of different formulations and constant density limit[END_REF]) 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. [START_REF] Bruel | A pseudocompressibility method for reactive flows at zero Mach number[END_REF]), unsteady turbulent premixed flames (Corvellec et al. [START_REF] Corvellec | A timeaccurate scheme for the calculations of unsteady reactive flows at low Mach number[END_REF], Dourado et al. [START_REF] Dourado | A timeaccurate pseudo-compressibility approach based on an unstructured hybrid finite volume technique applied to unsteady turbulent premixed flame propagation[END_REF]), 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], Donini et al. [START_REF] Donini | Buoyant Tsuji diffusion flames: global flame structure and flow field[END_REF]).

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: 𝜷 -𝟏 𝝏 𝝉 𝒑 + 𝝏 𝒕 𝝆 + 𝜵 ⋅ 𝝆𝒖 = 𝟎 (6) 𝝆𝝏 𝝉 𝒖 + 𝝆𝝏 𝒕 𝒖 + 𝝆𝒖 ⋅ 𝜵𝒖 = -𝜵𝒑 + 𝑹𝒆 -𝟏 𝜵 ⋅ 𝝉 ̿ + 𝑭𝒓 -𝟐 (𝟏 -𝝆)𝐞 𝒚 [START_REF] Dourado | A timeaccurate pseudo-compressibility approach based on an unstructured hybrid finite volume technique applied to unsteady turbulent premixed flame propagation[END_REF] 𝝆𝝏 𝝉 𝑻 + 𝝆𝝏 𝒕 𝑻 + 𝝆𝒖 ⋅ 𝜵𝑻 = 𝑷𝒆 -𝟏 𝜵 ⋅ (𝝆𝑫 𝑻 𝜵𝑻)

Where 𝒖=(𝑢, v) designates the dimensionless velocity vector, 𝜏̿ = μ(∇𝒖 + ∇𝒖 𝑻 ) -

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 [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 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: 𝛽 -1 𝜕 𝜏 𝑝 + 𝜕 𝑡 𝜌 + 𝛻 ⋅ 𝜌𝒖 = 0 (9) 𝜌𝜕 𝜏 𝒖 + 𝜌𝜕 𝑡 𝒖 + 𝜌𝒖 ⋅ 𝛻𝒖 = -𝛻𝑝 + 𝑅𝑒 -1 𝛻 ⋅ 𝜏̄̄+ 𝐹𝑟 -2 (1 -𝜌)𝒆 𝑦 [START_REF] Bassi | Artificial compressibility Godunov fluxes for variable density incompressible flows[END_REF] 𝜌𝜕 𝜏 𝑍 + 𝜌𝜕 𝑡 𝑍 + 𝜌𝒖 ⋅ 𝛻𝑍 = 𝑃𝑒 -1 𝛻 ⋅ (𝜌𝐷 𝑇 𝛻𝑍) (11) 𝜌𝜕 𝜏 𝐻 + 𝜌𝜕 𝑡 𝐻 + 𝜌𝒖 ⋅ 𝛻𝐻 = 𝑃𝑒 -1 𝛻 ⋅ (𝜌𝐷 𝑇 𝛻𝐻) (12) 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:

𝑑𝑞 (𝑎𝑐) 𝑑𝜏 + 𝑑𝑞 𝑑𝑡 = 𝑅𝐻𝑆(𝑞) 𝑖,𝑗 (13)
where 𝑞 = [𝜌, 𝑢, v, 𝑇] 𝑇 and 𝑞 (𝑎𝑐) = [𝑝/𝛽, 𝑢, v, 𝑇] 𝑇 are the vector of primitive variables for non reacting flows, and 𝑞 = [𝜌, 𝑢, v, 𝑍, 𝐻] 𝑇 and 𝑞 (𝑎𝑐) = [𝑝/𝛽, 𝑢, v, 𝑍, 𝐻] 𝑇 are those for reacting flows cases. 𝑅𝐻𝑆(𝑞) 𝑖,𝑗 is the right-hand side of the discretized equation. Then, introducing the residual value for the artificial time step as:

𝑅𝑒𝑠(𝑞 𝑛+1,ν ) 𝑖,𝑗 = ΔΩ -3𝑞 𝑛+1,ν + 4𝑞 𝑛 -𝑞 𝑛-1 2Δ𝑡 +𝑅𝐻𝑆(𝑞 𝑛+1,ν ) 𝑖,𝑗 (14)
permits to rewrite Equation 13 at times' step (𝑛 + 1, 𝜈 + 1) as:

𝑑𝑞 (𝑎𝑐) 𝑑𝜏 | 𝑛+1,𝜈+1 = 𝑅𝑒𝑠(𝑞 𝑛+1,𝜈 ) 𝑖,𝑗 (15)
The integration steps in artificial-time are finally given by: 𝑞 (1) = α 1 𝑞 𝑛+1,ν + ϕ 1 Δτ[𝑅𝑒𝑠(𝑞 𝑛+1,ν )/ΔΩ] 𝑞 (2) = α 2 𝑞 𝑛+1,ν + ϕ 2 [𝑞 (1) + Δτ𝑅𝑒𝑠(𝑞 (1) )/ΔΩ] 𝑞 𝑛+1,𝜈+1 = 𝛼 3 𝑞 𝑛+1,𝜈 + 𝜙 3 [𝑞 (2) + Δ𝜏𝑅𝑒𝑠(𝑞 (2) )/ΔΩ] 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 15is 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: 𝑦 ̂= 0: 𝑢 ̂𝑝𝑙𝑎𝑡𝑒 (0, 𝑡 ) = 𝑢 ̂𝑚𝑎𝑥 𝑐𝑜𝑠(2𝜋𝑓 ̂ 𝑡 ) [START_REF] Scha ¨fer | Benchmark computations of laminar flow around a cylinder[END_REF] 𝑓 ̂ 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 [START_REF] Sclichting | Boundary layer theory[END_REF]):

𝑢 ̂𝑓𝑙𝑢𝑖𝑑 (𝑦, 𝑡) = 𝑢 ̂max 𝑒 -𝑘 ̂𝑦 ̂𝑐𝑜𝑠(2𝜋𝑓 ̂ 𝑡 ̂-𝑘 ̂𝑦 ̂) ( 17
)
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:

‖𝑬 𝒏 ‖ 𝑝 = ((Δ𝑡, Δ𝑦) ∑ |𝐸 𝑛 𝑖,𝑗 | 𝑝 𝑖=𝑖𝑚𝑎𝑥,𝑗=𝑗𝑚𝑎𝑥 𝑖=1,𝑗=1 ) 1/𝑝 (18)
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 Scha ¨fer and Turek [START_REF] Scha ¨fer | Benchmark computations of laminar flow around a cylinder[END_REF] 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 presents the time evolution of the non-dimensional normal component of the velocity at (𝑥, 𝑦) = (-0.5, 0.5) observed when the periodic regime of s heddi ng is established. For 𝑅𝑒 = 100, the experimentally obtained Strouhal number is 𝑆𝑡 = 0.287 ± 0.003 [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 computed velocity is shown in 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 [START_REF] Bruel | A pseudocompressibility method for reactive flows at zero Mach number[END_REF][START_REF] Shapiro | Artificial compressibility, characteristics-based schemes for variable density, incompressible, multi-species flows. Part I. Derivation of different formulations and constant density limit[END_REF]. 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. [START_REF] Kim | A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations[END_REF] 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)/(𝜈α ̂)

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. [START_REF] Kim | A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations[END_REF] 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. [START_REF] Kim | A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations[END_REF] e.g. 𝑁𝑢 𝑇 = 1.67. Comparison with the numerical results of Kim et al. [START_REF] Kim | A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations[END_REF] obtained with a grid resolution of 201 × 201.

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 [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 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. and -𝑓𝑤 = 0.5. 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].

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 [START_REF] Tsa | Flame stabilization over a Tsuji burner by four-step chemical reaction[END_REF], 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%, 𝑌 ̂𝐹𝑟𝑒𝑓 ≡ 𝑌 ̂𝐹𝑏 = 0.5. The fuel and air burner inlet velocities are v ̂𝐹 = 10 × 10 -2 𝑚/𝑠 and v ̂𝑟𝑒𝑓 ≡ v ̂𝑎𝑖𝑟 = 15 × 10 -2 𝑚/𝑠, respectively. The resulted Reynolds number, Péclet number and Froude number for this case are 𝑅𝑒 = 122, 𝑃𝑒 = 86, and 𝐹𝑟 = 0.42. Also, the Figure 12 shows the following boundary conditions: on the top (𝑦 = 𝑦 ̂/𝑟̂= 40) and right (𝑥 = 𝑥 ̂/𝑟̂= 10) boundaries, free-slip condition is imposed for velocities along with zero fluxes for the temperature, mixture fraction and excess enthalpy. At the inlet section located at (𝑦 = 𝑦 ̂/𝑟̂= 0), the velocity streamwise component is prescribed for fuel and air, i.e. v(𝑥 ≤ 1) = v ̂𝐹/v ̂𝑟𝑒𝑓 = 0.66 and v(𝑥 > 1) = v ̂𝑎𝑖𝑟 /v ̂𝑟𝑒𝑓 = 1, respectively. The mesh dependency tests showed no significant differences in the large-scale flame features, such as flicker frequency and points of fuel pocket detachment. The Courant-Friedrichs-Lewy condition was the base reference for the choice of the numerical parameters. For the present case, the CFL numbers for the physical and artificial-time integration, were chosen as 𝐶𝐹𝐿 𝑡 = 0.27 and 𝐶𝐹𝐿 𝜏 = 0.2, respectively. Another key parameter of choice was the ability of the physical time-step to describe sufficiently well one cycle of flame flickering. 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. [START_REF] Davis | A numerical/experimental study of the dynamic structure of a buoyant jet diffusion flame[END_REF] who also used the flame sheet model, unit Lewis number hypothesis, and validated their results against the experimental data by Chen et al. [START_REF] Chen | Buoyant diffusion flames[END_REF]. 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 [START_REF] Toong | Mechanisms of combustion instability[END_REF]. 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 [START_REF] Chen | Buoyant diffusion flames[END_REF].The frequency analysis of this behaviour of the flickering jet diffusion flame is shown in Figure 13. 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 and14f. 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).

𝜔 ̂= 𝜕v ∂ 𝑥 ̂-𝜕𝑢 ∂ 𝑦 ̂ (20)

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 [START_REF] Chang | On the method of pseudocompressibility for numerically solving incompressible flows[END_REF]. 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.

Nomenclature

Figure 1 .

 1 Figure 1. Flow chart of the algorithm where 𝑡 𝑚𝑎𝑥 is the total physical time evolution simulated and 𝑞 is the vector of variables.

Figure 2 .

 2 Figure 2. Sketch of the fully staggered mesh used for the space discretization.

Figure 3 .

 3 Figure 3. Physical time and spatial order of accuracy by the log-log scale analysis using 1-norm.

3 2

 3 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 Figure4, Equations 6 and 7 are integrated with the following boundary conditions: on the top (𝑥 = 𝑥 ̂/𝑑 ̂= -2.1) and bottom (𝑥 = 𝑥 ̂/𝑑 ̂= 2) walls and at the cylinder surface (𝑟 2 = 𝑥 2 + 𝑦 2 ) the no-slip condition is imposed for velocities. At the inlet section located at (𝑦 = 𝑦 ̂/𝑑 ̂= -4), a parabolic profile is prescribed for the velocity streamwise component (with a maximum value v ̂(𝑥 ̂= ℎ ̂/2) = v ̂𝑚𝑎𝑥 =

Figure 4 .

 4 Figure 4. Unsteady cylinder wake: flow configuration and related boundary conditions.

Figure 5 .

 5 Figure 5. Unsteady cylinder wake (Mesh size of 43 × 246): (a) Velocity history signal at (𝑥, 𝑦) = (-0.5,0.5), (b) Corresponding power spectrum of such velocity signal, and (c) Snapshot of the contour of the velocity norm (The marker × indicates the location at which the data in (a) are recorded).

Figure 6 .

 6 Figure 6. Unsteady cylinder wake (Mesh size of 43 × 246): (a) Convergence history as a function of the dimensionless artificial compressibility factor 𝛽.

Figure 7 .

 7 Figure 7. Unsteady cylinder wake (Mesh size of 43 × 246): influence of the dimensionless artificial compressibility factor 𝛽 on the number of time steps required to obtain an artificialtime converged solution. Comparison with the lower bound proposed by Chang and Kwak [6].

  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:

Figure 8 .

 8 Figure 8. Heated cylinder in a square enclosure (Grid resolution of 101 × 101): Streamlines patterns.

Figure 9 .

 9 Figure 9. Heated cylinder in a square enclosure: Local Nusselt number distribution along the walls of the enclosure (Present results obtained with a grid resolution of 101 × 101).Comparison with the numerical results of Kim et al.[START_REF] Kim | A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations[END_REF] obtained with a grid resolution of 201 × 201.

Figure 10 .

 10 Figure 10. Tsuji diffusion flame: flow configuration and boundary conditions.

1 .

 1 Equations 9 to 12 are integrated with the following boundary conditions: On the symmetry axis (𝑥 = 0), 𝜕 𝑥 𝑢 = 𝜕 𝑥 v = 𝜕 𝑥 𝑍 = 𝜕 𝑥 𝐻 = 0; at the burner boundary surface (𝑟 2 | + = 𝑥 2 + 𝑦 2 = 1 + ), 𝑢 -𝑥 = v -𝑦 = 0, 𝑍 𝑠 -𝑃𝑒 -1 𝜕 𝑛 𝑍 = 𝑆 𝑍 , 𝐻 𝑠 -𝑃𝑒 -1 𝜕 𝑛 𝐻 = 𝑆 𝐻 (Robin's like boundary type for 𝑍 and 𝐻 ) where 𝑍 𝑠 ≡ 𝑆 𝑌 𝐹 𝑠 -1, 𝐻 𝑠 ≡ (𝑆 + 1) 𝑇 𝑠 /𝑄 + 𝑌 𝐹 𝑠 and the subscript 𝑛 stands for the normal to the burner surface. The terms 𝑆 𝑍 and 𝑆 𝐻 are the 𝑍 and 𝐻 fluxes which are imposed at the burner injection surface 𝑟 2 | -= 𝑥 2 + 𝑦 2 = 1 -as function of 𝑌 ̂𝐹𝑏 , 𝑇 ̂𝑏 and 𝑢 ̂𝑏, namely 𝑆 𝑍 ≡ 𝑆 + 1 and 𝑆 𝐻 ≡ (𝑆 + 1) 𝑇 𝑏 /𝑄 + Note that 𝑌 𝐹 𝑠 and 𝑇 𝑠 are found as part of the solution of the problem and this holds only in the forward part of the cylinder. The boundary conditions at the inlet (𝑦 = -7.5) are v = 1, 𝑢 = 𝑍 = 0 and 𝐻 = (𝑆 + 1) 𝑇 ∞ /𝑄 + 𝑌 𝑂 . At the outlet (𝑦 = 13), they are 𝜕 𝑦 𝑢 = 𝜕 𝑦 v = 𝜕 𝑦 𝑍 = 𝜕 𝑦 𝐻 = 0 and at the channel wall (𝑥 = 4), they read 𝑢 = v = 𝜕 𝑥 𝑍 = 𝜕 𝑥 𝐻 = 0. According to the definition of the mixture fraction function 𝑍, the flame position (𝑥 𝑓 , 𝑦 𝑓 ) is given by the isoline 𝑍(𝑥, 𝑦) = 1 where the flame temperature 𝑇 𝑓 is determined by 𝐻(𝑥 𝑓 , 𝑦 𝑓 ) = (𝑆 + 1) 𝑇 𝑓 /𝑄. The steady diffusion flame results are presented in 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.

Figure 11 .

 11 Figure 11. Tsuji diffusion flame (Mesh size of 82 × 446): (a) Temperature distribution through the flame front of a Tsuji burner with 𝑅𝑒 = 38, 𝐹𝑟 = 2.6, 𝑟̂= 0.02𝑚, 𝑢 ̂𝑎𝑖𝑟 = 1.15𝑚/𝑠, and -𝑓𝑤 = 0.318 and (b) streamlines and flame shape with 𝑅𝑒 = 18, 𝐹𝑟 = 1.95, 𝑟̂= 0.015𝑚, 𝑢 ̂𝑎𝑖𝑟 = 0.75 𝑚/𝑠,and -𝑓𝑤 = 0.5. 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].

Figure 12 .

 12 Figure 12. Flickering diffusion flame: flow configuration and boundary conditions.

Figure 13 .

 13 Figure 13. Flickering diffusion flame (Mesh size of 41 × 101): (a) Mixture Fraction history at (𝑥 ̂, 𝑦 ̂) = (1.3 × 10 -2 𝑚, 10 -1 𝑚), (b) power spectrum for the flame fluctuation. The mixture fraction history was chosen to represent the flame fluctuations at (𝑥 ̂, 𝑦 ̂) = (1.3 × 10 -2 𝑚, 10 -1 𝑚), as shown in Figure 13a. The dashed line represents the stoichiometric value, thus, when the continuous line crosses this line, it means that the detachment of a hot pocket of fuel occurred. The power spectrum derived from the flame fluctuations are presented in Figure 13b. The numerically computed predominant frequency at the probe location of 𝑓 ̂= 𝑆𝑡/𝑡 ̂𝑟𝑒𝑓 = 8.4𝐻𝑧 agrees well (relative error of 5.5%) with the value predicted by the correlation 𝑓 -𝐹𝑟, 𝑓 ̂= 𝑡 ̂𝑟𝑒𝑓 0.29𝐹𝑟 -1 = 7.96𝐻𝑧, suggested by Sato et al.[START_REF] Sato | Diffusion flames and their flickering motions related with Froude numbers under various gravity levels[END_REF]. 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 Figure13bprobably 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 Figure14which displays the temperatures isocontours and vorticity contours where the dimensional vorticity is defined by:

Figure 14 .

 14 Figure 14. (a) 𝑡 ̂= 1.3𝑠, (b) 𝑡 ̂= 1.35𝑠, (c) 𝑡 ̂= 1.39𝑠, (d) 𝑡 ̂= 1.43𝑠, (e) 𝑡 ̂= 1.47𝑠, (f) 𝑡 ̂= 1.52𝑠. The red isoline represents the flame sheet (𝑍 = 1) and the blue × marker represents the location of the probe from Figure 13a.

  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: 𝛽 ≡ 𝛽 ̂/(𝑢 ̂𝑟𝑒𝑓 ) 2 , 𝑡 ≡ 𝑡 (𝑢 ̂𝑟𝑒𝑓 /𝐿 ̂𝑟𝑒𝑓 ), 𝜏 ≡
	𝜏 (𝑢 ̂𝑟𝑒𝑓 /𝐿 ̂𝑟𝑒𝑓 ) , 𝑢 ≡ 𝑢 ̂/𝑢 ̂𝑟𝑒𝑓 , v ≡ v ̂/𝑢 ̂𝑟𝑒𝑓 , 𝜌 ≡ 𝜌 ̂/𝜌 ̂𝑟𝑒𝑓 , 𝑇 ≡
	𝑇 ̂/𝑇 ̂𝑟𝑒𝑓 and 𝑝 ≡ 𝑝/𝜌 ̂𝑟𝑒𝑓 𝑢 ̂𝑟𝑒𝑓 2 , 𝛻 = 𝐿 𝑟𝑒𝑓 𝛻 ̂. The reference
	quantities used to define the Reynolds number
	𝑅𝑒 ≡ 𝑢 ̂𝑟𝑒𝑓 𝐿 ̂𝑟𝑒𝑓 /𝜈̂𝑟 𝑒𝑓 , the Froude number 𝐹𝑟 ≡ 𝑢 ̂𝑟𝑒𝑓 /
	√𝐿 ̂𝑟𝑒𝑓 ‖𝒈 ̂‖ and the Péclet number 𝑃𝑒 ≡ 𝐿 ̂𝑟𝑒𝑓 𝑢 ̂𝑟𝑒𝑓 /𝛼 ̂𝑟𝑒𝑓
	(with 𝜈̂𝑟 𝑒𝑓 ≡ 𝜇̂𝑟 𝑒𝑓 /𝜌 ̂𝑟𝑒𝑓 , 𝛼 ̂𝑟𝑒𝑓 ≡ 𝑘 ̂𝑟𝑒𝑓 /𝑐̂𝑝 𝑟𝑒𝑓 𝜌 ̂𝑟𝑒𝑓 , where
	𝜇̂𝑟 𝑒𝑓 is the reference dynamic viscosity, 𝛼 ̂𝑟𝑒𝑓 the
	reference thermal diffusivity 𝑘 ̂𝑟𝑒𝑓 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 𝑇 = 𝑐𝑠𝑡𝑒.

μ (𝛻 ⋅ 𝒖)𝑰 is the viscous stress tensor and the power-law expressions μ = ρ𝐷 𝑇 = 𝑇 σ (with 𝜎 = 0.7) are employed to account for the

Acknowledgments

Part of this work was done during a 1-year stay of the first author in France at Laboratoire de Mathématiques et de leurs Applications (CNRS/UPPA), Pau. Mariovane Donini gratefully acknowledges Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) for its financial contribution to his research stay in France under grants numbers 88887.363086/2019-00 and 88882.444509/2019-01.

Abbreviations

AC

artificial compressibility CFD computational fluid dynamics


Dr Michel Develoux
email: micheldeveloux@yahoo.fr

F M Amona

C Hennequin

Pr C Hennequin

Histoplasmosis caused by Histoplasma capsulatum var. duboisii: a comprehensive review of the cases from 1993 to 2019

Keywords: African histoplasmosis, Histoplasma capsulatum var. duboisii, epidemiology, HIV, diagnostic, treatment Word count: 2908

Histoplasmosis Caused by Histoplasma capsulatum var. duboisii :

Introduction

H. capsulatum, a dimorphic fungus responsible for human histoplasmosis, encompasses two varieties, namely H. capsulatum var. duboisii (Hcd) and H. capsulatum var. capsulatum (Hcc) [1]. Dubois et al first described Hcd in 1952 in the Belgian Congo, today the Republic of Congo [2]. A year later Vanbreusegem used the parasitic and saprophytic characters to distinguish between the two varieties [3]. Furthermore, epidemiological and clinical characteristics have been described for each of the infections due to the 2 varieties [4,5]. Notably since Hcd infection has only been reported from Africa, the disease is commonly referred as "African histoplasmosis". However, both varieties has been described in Africa and recent phylogenetic analyses support the monophyly of H. capsulatum on this continent, with the existence of a unique African clade including both Hcc and Hcd strains [6,7]. Thus, these results and other suggest that the classification, based only on the geographic origin of the acquisition, is obsolete [8].

In Africa, Hcc and Hcd are endemic in at least 32 countries, but Hcd infection reports overpass those on Hcc infection [5,9]. Nevertheless, Hcd histoplasmosis is considered as a rare fungal infection, occurring in Sub-Saharan African countries and Madagascar, but the disease is not included in the group of the Neglected Tropical Diseases [9]. Since the first description of the disease, several papers over the years have reviewed its main characteristics [10][11][12]. In 1993, Carme et al recorded 237 Hcd histoplasmosis cases without detailing their characteristics [13]. The aim of the present review is to stand out the possible new aspects regarding this infection by the mean of a large review of the literature from 1993.

Notably, the role of HIV co-infection on the clinical presentation and outcome is analyzed.

Methods

For the present study, we searched for new cases or series published after the publication of Carme et al in 1993 [13]. This was done using the PubMed database with the limitations of date 1993-2019 and with the keywords: "African histoplasmosis", "histoplasmosis AND Africa", "Histoplasma duboisii", "Histoplasma capsulatum var duboisii". We included all publications whatever was the language, notably French and Portuguese ones. In addition, we retrieved from other internet sources (Google scholar notably) some cases published in non-referenced journals.

We only considered well-documented cases and so excluded cases described in conference abstract books or presented in case-series with limited details. For each case, we analyzed the demographic information, the clinical presentation, results of radiological and microbiological investigations, the mode of diagnosis and the therapeutic management. We focused particularly on cases occurring in HIVinfected patients, as they emerged as an at-risk population in our analysis.

Statistics

Comparisons were done by the mean of the Fisher's exact or Chi-square test for qualitative variables and the Student t-test for quantitative variables using the Prism (v7.0a) software.

Results

We found a total of 147 cases with a clear demonstration of yeasts supporting the diagnosis of Hcd infection (Supplementary Table S1). Seventeen and 36 cases were included in the two most important case-series reported from Togo and Republic of Congo,respectively [14,15]. Those cases were excluded from the clinical and therapeutic parts of the analysis due to insufficient description of the cases. Thus, 94 well-documented cases were reviewed to depict at the best the main characteristics of Hcd infections.

Epidemiology

The sex-ratio (M/F) was at 2.4 and the mean age at age at 30.7 ± 17.5 ranging from 13 months to 68 years. Children under 18 represented 30.1% of the cases. Patients originated from Central and Western Africa in 141 out of 147 cases (Figure 1).

Clinical presentation

The main clinical presentations are summarized in Table 1 and illustrated in Figure 2. Cutaneous and subcutaneous lesions, osteoarticular infections and lymphadenopathies were the most frequent, isolated or combined, and seen in 69 (74.2%), 44 (46.8%) and 39 cases (41.9%), respectively. Twenty-seven patients were febrile at the time of diagnosis. In 37 cases (39.4%), the infection presented as a single lesion, mainly localized either to the skin (n=13), or to a bone (n=9) but 57 patients (60.6%) had disseminated infection affecting at least 2 non-contiguous body-sites. Compared to a 1972 report from Drouhet, the frequency of disseminated form increased from 16.4% of the cases (19 out of 116) [12] to 60.6% in our review (Fischer's exact test; p<0.0001). Disseminated forms were not significantly more frequent in HIV-positive (81.8% versus 62.1%) (Fisher's exact test p=0.11) (Table 2).

The cutaneous lesions could be either unique (n=31) or multiple (n=38). They involved preferentially the limbs (n=68), the trunk (n=54), and the face (n=59). Their appearance was very diverse and sometimes polymorph for a given patient. Cutaneous and subcutaneous nodules, maculo-papules, abscesses were the most common with the multiple umbilicated papulo-nodules simulating molluscum contagiosum observed in 21 patients in the case of dissemination.

Lymphadenopathies were either superficial (n=37) and/or deep-seated (n=6). Superficial lymphadenopathies could be firm or indolent, fixed or mobile sometimes painful and localized with equal frequency in the inguinal (n=23), axillary (n=24) and cervical chains (n=24). A fistula to the skin was seen in 10 cases.

Osseous infection was frequently multifocal (25 out of 44 cases). Infections of the spine (n=14), the thorax (n=15), bones of the upper (n=18) and lowers (n=18) limbs were the most common. In 11 cases spinal infection mimicked a Pott's disease [16][17][18]. Lesions were always osteolytic without reconstruction. In 4 cases, they were revealed or complicated by a pathologic fracture [19][20][21][22].

Lungs infections was clinically or radiologically suspected in 10 cases but mycologically confirmed in only 3 cases. Adrenal gland infection and nasal or buccal mucosal infection was reported in 2 and 4 cases, respectively.

Abdominal infection was diagnosed in 8 patients (8.5%) either in the form of gastro-intestinal tract infection (gastric or duodenal lesions (n=2), acute colitis (n= 1)) or as intra-abdominal infection such as mass (n=2), abscess (n=1), peritoneal nodules or peritonitis (n=2).

Gathering the data of cases reports and of the more important case-series [14,15], 25 out of 120 (20.8%) tested patients were found HIV-positive. Fever was most common in HIV+ patients (59.1 vs 22.4%;

Fisher's exact test p=0.003) and we found a lower frequency in bone infection 18.2 vs 58.6% (p=0.002; Fisher's exact test) and a trend toward a higher frequency of lymph nodes infection 63.6 vs 38.6% (Fisher's exact test; p=0.08) in those patients compared to HIV-patients.

Laboratory diagnosis

Among the 94 observations reviewed, the diagnosis relied on direct examination of different fluids or skin scrapings, or histopathological examination of a biopsy specimen in 50 and 74 cases, respectively.

A confirmation of the diagnosis using either culture or PCR was only obtained in 26 and 7 cases, respectively. Among the 9 patients tested for serum antibodies, 3 were found positive. We did not find any case mentioning Histoplasma antigen detection neither ß-D glucan assay. A single case of Hcd disseminated infection in an HIV infected woman reported with a positive galactomannan (GM) serum test, with a slow decrease in titer following antifungal treatment [23]. GM test was mentioned in another case but the assay remained negative in an immunocompetent woman with pulmonary and cutaneous lesions [24].

Treatment

The treatment regimen was specified for 90 out of 94 observations analyzed. Eight patients did not receive any treatment either because of a rapidly fatal outcome (n=5), because they denied treatment (n=1) or because of spontaneous healing or absence of relapse following surgery alone (n=2). Thus, 82 patients received antifungal drug(s), either alone (n=59) or combined with surgery (drainage, debridement of necrotic tissues, osseous curettage) (n=23). First-line therapy mainly relied on iv amphotericin B (n=28) or one of its lipid formulations (n=5), itraconazole (n=27), fluconazole (n=8) and ketoconazole (n=13), the latter being no longer recommended for systemic therapy (Table 3). In 7 cases, amphotericin B toxicity (renal failure, general symptoms) led to the discontinuation or a switch for another drug. The safety profile of itraconazole was overall good with only 2 cases of adverse events mentioned in 45 patients treated with this drug. Two patients with bone infection were treated with terbinafine, either as first-line or as maintenance therapy following iv amphotericin B, leading to remission in both cases with a follow-up of 5 and 8 months, respectively [21,25].

By the end of the first-line treatment, remission, failure and relapse were notified in 50 (53.2%), 6 (6.4%) and 9 (9.6%) cases, respectively. We failed to demonstrate any significant difference in the success rate according to the antifungal drugs, but the mortality rate was significantly higher for patients that did not receive any antifungal drug as compared to those treated with antifungal(s): 75% vs 13.4% (p=0.0004 Chi-square).

By the end of follow-up, 22 (23.4%) patients died but the final outcome was not specified for 18 cases and only 17 had a follow-up exceeding one year. We found a trend toward a better prognosis for localized infection compared to disseminated form (Fisher's exact test; p= 0.0685), which was not the case for the HIV serologic status which did not influence the outcome (p=0.76).

Discussion

While not include in the official list of the WHO, Hcd histoplasmosis can be considered a neglected tropical disease. From the first report in 1952, nearly 400 Hcd cases have been reported but the cases series are few so that the epidemiologic and clinical characteristics and the optimal therapeutic management remain incompletely known. In this review, we analyzed retrospectively 94 welldocumented cases allowing us to point out some of those characteristics.

As described in the past [10,12], the majority of documented cases arose from Western and Central Africa, between latitudes 15° North and 10° South (Figure 1). However, a recruitment bias is likely as most of the cases are recorded from countries having long-standing medical structures able to diagnose the disease. More recently, cases have been reported from other countries such as Togo or Guinea-Bissau [14,26], suggesting that the endemic area could be wider. However, given the coexistence of Hcc and Hcd in Africa, the use of the results of ecological studies (detection of the fungus in the soil or of guano specimens) or of intradermal tests of histoplasmin to specify the endemic area should be done with caution.

The hypothesis that the prevalence of the disease is underestimated is also supported by the fact that 28 cases out of 94 (29.8%) have been diagnosed in countries outside Africa, mainly in Europe. The rural or urban origin of the patients was rarely specified in the literature, but people living in the bush or the forest, with outdoor activities, seemed to be the most affected [13,27]. In a Togolese series, the majority of patients came from rural areas (15 out of 17 cases) just as in the Congolese series (10 out of 14 cases).

Most had a low socio-economic level [14,15].

The pathophysiology of Hcd infection remains largely unknown, notably its mode of acquisition.

Inhalation of fungal spores is a possible route of infection. During the investigation of a particular focus of Hcd infection in the Democratic Republic of the Congo, four interviewed patients acknowledged collecting guano from caves to fertilize gardens [15]. However, only 3 out of the 94 cases we reviewed were diagnosed with a documented pulmonary infection [23,28,29]. A breach in the skin barrier is another possible mode of acquisition with several cases supporting this hypothesis: a man developed an ulcerated tumor at the needle prick of acupuncture after mud-baths [30], and another had superinfection of cicatricial lesions following tribal scarification rituals [31]. Similarly, the infection has been sometimes described after a cutaneous traumatism with initial symptoms at the site of the trauma [19,32,33]. In some cases, an intestinal penetration could have been the portal of entry as cases with isolated digestive involvement have been reported [31,34,35].

Our review does not allow a precise estimation of the incubation time of Hcd infection. Cutaneous signs at the level of a local trauma have been reported 3 weeks to 2 months later [19,32]. In contrast, HIVnegative expatriates developed symptoms years (a minimum of 2) or decades after leaving the endemic area [31,[35][36][37][38][39]. So, whatever the date of leaving the endemic area is, the diagnosis of Hcd infection should not be ruled out before appropriate investigation.

Analysis of clinical charts highlights the very pleiomorphic presentation of the disease, so that the disease should be more commonly suspected. However, the most striking result is the increase in the incidence of disseminated infection over the 25 years of our review. The HIV co-infection has been previously suggested to explain this phenomenon [38]. However, we failed to demonstrate a significative difference in the incidence of disseminated form between HIV+ and HIV-patients. So, one can assume that this phenomenon could also be due to the more frequent use of high definition imaging, allowing the detection of lesion(s) that would have been previously ignored. In their review, Loulergue et al noted that fever (58% versus 15%) and the lack of lung infection (0% versus 75%) were more common in HIV-positive patients [38]. In addition, we were able to point out the lower frequency of bone infection and the higher frequency of lymphadenopathy in those patients.

Probably due to the limited laboratory diagnosis opportunities in low-income countries, 55% of the cases were diagnosed only on direct examination or histopathological analysis of clinical specimen(s) making some diagnosis questionable as already mentioned [40]. This is illustrated by the unique autochthonous case of Hcd infection diagnosed outside of Africa [41]. This patient from Southern India denied any travel outside this region. He presented with multiple erythematous papules on the trunk and the arms.

The authors presented histopathological findings compatible with Histoplasma yeast forms. Cultures grew positive for a fungus further identified as H. capsulatum. This case underlines the difficulty in the definitive diagnosis of Hcd infection. Indeed, Hcc being, like in Africa, endemic in India, the diagnosis of Hcd infection only relies on the size of yeasts observed on the biopsy specimen, which can be questionable in this case. Several reports support the use of PCR for the diagnosis of Hcd infection [30,42] but it should be recalled that neither the PCR-based methods nor the morphology of the colonies of obtained in culture, irrespective of the phase obtained either the mycelial or the yeast phase, do allow the differentiation between Hcc and Hcd. Thus, a careful examination to characterize the morphology of the yeasts in tissue or fluids is crucial for a definitive diagnosis [40]. Large budding yeasts of 7-15 µm in diameter, ovoid, with thick cell walls are characteristic features [4]. Chain formation in the shape of an hourglass or figure-eight, bud scar and cytoplasmic lipid inclusions also support the diagnosis, but while these features are quite easy to recognize on direct examination of fluids, pathologic examination of biopsy specimen can be more difficult to interpret. This warrants the development of new diagnostic tests reliable and easy to perform in low-income countries.

To the best of our knowledge, there is no specific guidelines or consensus statement for the therapeutic management of Hcd infection. By analogy with Hcc infection, and in accordance with the IDSA guidelines, some experts recommend the use of iv amphotericin B deoxycholate or its lipid formulations, or itraconazole, according to the severity of the disease [38,43]. However, it should be kept in mind that the patient's financial resources and the availability of drugs in the considered country probably most often drove the therapeutic option. Amphotericin B is not licensed in or unavailable in 22 (14.2%) and 42 (27.1%) out of 155 countries, respectively [44] and the daily price of deoxycholate amphotericin B ranges from <$1 to $171. Considering the lack of significative difference in the outcome according to the drug used, itraconazole, thanks to its oral formulation and its safety profile, appears as a valuable alternative. Therapeutic drug monitoring is desirable notably in the case of co-administration of enzyme inductor drugs such as rifampicin, but not widely available in low-income countries. Fluconazole has the advantage of being well-tolerated but its bone penetration is lower than itraconazole [45].

The optimal duration of therapy could neither be defined. Relapses have been reported up to 2 years after an initial improvement [23]. Bone lesions should probably be treated with a minimum course of 6 months, which represents a major drawback in the limited-resource setting. For HIV-infected patients with Hcd infection, the duration of therapy should depend obviously on the recovery of an efficient immune response. Immune reconstitution inflammatory syndrome has been reported in two patients with Hcd infection [46,47], but recent guidelines do not recommend to delay the initiation of antiretroviral therapy in the context of opportunistic histoplasmosis (due to Hcc) [48].

Conclusions

Our review suggests that the prevalence of Hcd infection is underestimated because cases can be easily missed notably in low-income countries. Indeed, the clinical presentation is very pleiomorphic so that the disease should be suspected more frequently. HIV co-infection is common (20% of the cases) and substantially modifies this clinical presentation. Diagnosis mainly relies on the direct examination which requires specific laboratory materials and microscopist skills, that are not always available in those countries. There is an urgent need for reliable and easy-to-use diagnostic tools to confirm the diagnosis and distinguish between Hcd and Hcc infections, which are both endemic on the African continent.

Finally, the review points out the absolute necessity for the availability at a reduced price of the antifungal drugs used to treat Hcd infection, mainly iv amphotericin B and itraconazole.

Cutaneous/Subcutaneous lesion(s) (n=69 cases)

Face

Table 1 :

 1 Main clinical characteristics of 94 patients with Hcd infection.

	14

Table 2 :

 2 Comparison of clinical features of 94 patients with Hcd infection according to the HIV

	serologic status
	NS: not significative
	1 HIV serologic status was not specified for 14 patients
	2 Intra-abdominal infection (n=3); Adrenal glands (n=2), Lungs (n=1), Digestive tract (n=2), mucosal
	infection (n=2)
	3 Bone marrow (n=2), Digestive tract (n=1), Intra-abdominal infection (n=1), mucosal infection (n=2),
	Lungs + digestive tract (n=1), Lungs + liver + spleen + bone marrow (n=1)


Jean-Marc Salotti
email: jean-marc.salotti@ensc.fr

Michel Lamontagne

Phobos, Possible Refueling Station for a Single Stage to Orbit Vehicle

Keywords: Mars settlement, Phobos, SSTO

Numerous architectures have been proposed for a human mission to Mars. Most of them are appropriate for the first missions but in the long term, important improvements are expected to reduce the complexity and costs of the missions. Assuming the presence of water and carbonates on Phobos, which is suggested by several experts, it might be possible to use this Moon of Mars as a refueling station for rockets coming from Earth and Mars. In the long term, an interesting perspective is to manufacture on Mars a single stage to orbit space plane, called Albatros, that could be used to transfer crews and cargo from Mars orbit to Mars surface and the other way round. Technical specifications and illustrations of this SSTO are given as well as a range of possible applications. This idea was developed in a proposal in response to the recent Mars Society City State Design Competition, where we achieved the third place, with "Foundation".

Introduction

The internal composition of Phobos is uncertain. In several papers, it is suggested that Phobos is an old moon of Mars and has plenty of interesting resources like water and carbonates, which could be exploited to produce propellant [START_REF] Pollack | Multicolor observations of Phobos with the Viking lander cameras-evidence for a carbonaceous chondritic composition[END_REF][START_REF] Rivkin | Near-Infrared spectrophotometry of Phobos and Deimos[END_REF][START_REF] Fanale | Evolution of the water regime of Phobos[END_REF][START_REF] Ronnet | Reconciling the Orbital and Physical Properties of the Martian Moons[END_REF]. If it is confirmed, Phobos could play the role of a refueling station for vehicles coming from the Earth or from the Martian surface. For the first human missions, the construction of a space station on Phobos would certainly be inappropriate, but in the long term, assuming a constant human presence on the red planet and an interest in optimizing transportation means and reducing costs, the idea deserves to be developed. It is proposed in this paper to investigate the idea and to determine what could be the role of Phobos exactly. Section 2 is a review of Phobos' physical features. Section 3 presents the concept of a space station on Phobos. It would be used as a refueling station for vehicles coming from Earth. Passengers and cargo would be transferred into reusable space vehicles that would land on the surface. Section 4 shows the design of the reusable space vehicle proposed. It is a single stage to orbit rocket called Albatross. In the long term, as such a vehicle should be built and assembled on Mars. The feasibility, industrial requirements and maintenance issues are discussed.

This paper is based on a concept that was originally developed for "Foundation", a submittal to the 2020 Mars City State Design Competition organized by the Mars Society [START_REF] Lamontagne | Foundation[END_REF].

Phobos and Deimos

Chemical composition

Phobos and Deimos are the two moons of Mars. Despite numerous spectral observations, the composition of the moons is unclear [START_REF] Pollack | Multicolor observations of Phobos with the Viking lander cameras-evidence for a carbonaceous chondritic composition[END_REF][START_REF] Rivkin | Near-Infrared spectrophotometry of Phobos and Deimos[END_REF]. In 1990, Fanale suggested that there is possibly water under the surface of Phobos [START_REF] Fanale | Evolution of the water regime of Phobos[END_REF]. According to the most recent studies, the composition is tightly linked to the origin of the moons [START_REF] Ronnet | Reconciling the Orbital and Physical Properties of the Martian Moons[END_REF][START_REF] Rosenblatt | On the formation of the martian moons from a circummartian accretion disk[END_REF]. There are two hypotheses: The first is the capture of an asteroid and the second is the accretion of debris after a giant impact on the surface of Mars that occured billions of years ago. As their orbits are close to a circle and are contained in the equatorial plane of the planet, the second hypothesis is currently favoured. In this case the composition of the moons would be a mixture of material that originated from the impacting asteroid and the surface of Mars. Depending on the nature of the asteroid (carbonaceous or not), the ratios of the different chemical elements present would differ. However, as some elements come from the Martian surface, numerous light chemical elements would be present (e.g., carbonates and water) in non-negligible amounts. An important question, however, is to know why no trace of water is currently observed on the surface. According to Rosenblatt et al, volatile elements present on the surface of the moons may have been lost [START_REF] Ronnet | Reconciling the Orbital and Physical Properties of the Martian Moons[END_REF]. For that reason, carbonates and ice might not be present directly on the surface but might be found underground, as suggested by Canale [START_REF] Fanale | Evolution of the water regime of Phobos[END_REF].

IAC-21,A5,2,6,x62904
Page 2 of 6

Physical parameters

Deimos is orbiting very far from Mars, about 23 thousand kilometers, while Phobos is orbiting much closer at about 9.3 thousand kilometers. As the objective is to study the feasibility of a stopover space station for human missions to Mars, it is preferable to minimize the time and energy for the descent to the surface of Mars, as well as the time and energy for the ascent to the moon. Phobos is therefore the best candidate for a space station. 3. Phobos as a refueling station and hub

Stopover before landing on Mars

In most human missions to Mars architectures, for example in NASA Design Reference Architecture or Salotti's version of Mars Semi-Direct, the mass of the landing vehicle is in first order approximation half the mass of the vehicle arriving in the vicinity of Mars [START_REF] Drake | Human Exploration of Mars[END_REF][START_REF] Salotti | Robust, affordable, semi-direct Mars mission[END_REF]. The surplus is the mass of all systems required for entry, descent and landing:

-Heat shields.

-Thermal control system.

-Propulsion system for the last phase of the descent, engines, tanks and propellant. -Landing legs. If it is not required to land on the surface of Mars but to go to Mars orbit and land on Phobos, some mass savings are expected. The Delta V budget to go to Phobos is presented Table 2 [START_REF] Hopkins | Comparison of Deimos and Phobos as Destinations for Human Exploration and Identification of Preferred Landing Sites[END_REF][START_REF] Rebuffat | Phobos Sample Return[END_REF]. There are two cases:

• Without aerobraking, the Delta V depends on several parameters, such as the planetary configuration, which could imply a plane change and the type of transfer (Hohmann or faster). It is about 2 km/s if no effort is made to speed up the interplanetary trip. Assuming an Isp of 350s (methane oxygen), the mass of propellant would be around 44% of the total mass (Tsiolkovsky equation) and the mass of the wet propulsion system would be around 50%. It could eventually be much higher if a fast Earth Mars transit is implemented. This option is therefore not interesting in comparison with a direct to Mars surface option. • With aerobraking, the propulsive Delta V is rather limited, on the order of 1 km/s, even if a fast Earth Mars transit is implemented. It could even be lowered if the aerocapture maneuver is stronger and the first apoapsis close to the altitude of Phobos. However, an aerocapture maneuver can also be used for a direct to Mars option. The descent is demanding in terms of heat shields and thermal control, but they are already needed for aerocapture and the resulting Delta V to reach the surface would be more or less the same, on the order of 1 km/s.

A small gain might still exist, but a stopover at Phobos for refueling or for the transfer onto another vehicle before going to the surface is therefore not very interesting.

Refueling on Phobos for the return

Another important issue is the return to Earth. If a vehicle does not stopover at a space station in Mars orbit, it must bring all the propellant required for the return, which in turn would have an important impact on the mass of the propellant and systems for Mars orbit insertion. In the NASA reference mission, the impact depends on the specific impulse of the propulsion system [Nasa]. In general, from the chosen orbit, there is 1.5 km/s for trans-Earth injection. For chemical engines, the mass of propellant represents about 30% of the total mass of the return vehicle. There is however an additional impact if there is a propulsive Mars orbit insertion. In Salotti Mars semi-direct architecture, in which aerocapture is used for all vehicles, the mass of propellant for the return also represents about 30% of the total mass of that vehicle and there would be some mass penalty for scaling up the aerocapture systems. If some propellant can be manufactured on a moon of IAC-21,A5,2,6,x62904

Page 3 of 6

Mars, it would therefore allow an important mass saving.

Refueling on Phobos for the next flight from Earth

In addition to that, other mission architectures deserve to be examined. Let us consider the architecture proposed by Space X with an interplanetary reusable vehicle called Starship. Though still in early development, according to Musk, 5 heavy launches would be required to refuel the Starship vehicle in LEO before each new mission to Mars. Using Phobos as a refueling station, it could indeed be possible to avoid two or three heavy launches to LEO by refueling on Phobos during the previous Mars-Earth flight. The tanks would be at least half full after trans-Earth injection, and this would allow huge mass savings for the preparation of the next flight. If two or three heavy launches were eliminated, the cost could probably be cut by a factor of two. Starships could still go to Mars first, but on the way back, they would stop over at Phobos for refueling. Another idea is to use other vehicles for the transitions between Phobos and Mars surface and the other way round. If there is a station on Phobos, there would certainly be a need to go to the station at any time for maintenance or for swapping working teams. The idea of a Single Stage To Orbit (SSTO) rocket is therefore interesting and further developed in the next sections.

Phobos as a spaceport hub

In the future, numerous martian bases will be built on the surface and some of them might be very distant from each other. In order to reduce the travelling time to each destination, Phobos could be used as a hub, in which new visitors or settlers could make a rocket change for a direct descent to a very close spaceport.

Interesting economic perspectives

Refueling on Phobos allows other interesting and competitive options for many destinations. Considering the long term perspective of a large settlement on Mars, the most important could be the design and production of satellites on Mars to be sent to high Earth orbits (e.g., geostationary). The transfer from Mars could be made in two steps: First sending the satellite to Phobos station using a SSTO (see next sections), then refueling and sending the satellite to Earth orbit. In terms of Delta V, it is indeed easier to send satellites from Mars to GTO than from Earth, with 8,7 km/s vs 13.1 km/s. See Table 3. And thanks to the refueling station, the maximum Delta V before refueling is on the order of 6 km/s. For such a low requirement, a reusable SSTO (Mars to GTO through Phobos and back) could be exploited and could help in reducing the total cost of the missions. This is obviously for the long term, as the production of satellites requires important industries on the surface of Mars.

Other missions could also be interesting with a stopover at Phobos: missions to asteroids, missions to the Moon, robotic and human missions to other destinations in the solar system, like the moons of Jupiter or even further. In addition, if there are several permanent bases on the surface of Mars, a stop at Phobos station could be interesting to take a space taxi to the exact destination on the surface.

Propellant production

Let us assume, as suggested by Canale, that Phobos contains carbonates and water a few meters below the surface. There are different ways to produce propellant using these resources. As proposed by several authors for a Mars ascent vehicle, it is possible to produce hydrogen, methane, or ethanol, and to combine them with oxygen to get some thrust for a rocket [START_REF] Zubrin | Integrated Mars In Situ Propellant Production System[END_REF][START_REF] Drake | Mars Architecture Steering Group, Human Exploration of Mars[END_REF]. Among the different options, the easiest is probably a rocket engine powered by cryogenic liquid methane and liquid oxygen.

• Oxygen can be obtained from water electrolysis and some energy. If large amounts of ice can be found, there is no difficulty to get that element, together with hydrogen. • Methane can be obtained in several steps. First, the idea is to extract carbonates from the Phobos underground. Second, carbon dioxide can be produced by heating. Third, a Sabatier reactor can be used to transform carbon dioxide and hydrogen into methane. In order to implement these reactions, a chemical unit is required together with a power plant and several other systems such as compressors, collectors and heaters.

IAC-21,A5,2,6,x62904 Page 4 of 6 On Earth, there have been many designs for a SSTO rocket and a single attempt to build one, the NASA X33 spaceplane. A precursor for the larger Venture Star design, it had a vertical launch and a horizontal landing. Its failure was linked to the construction of its fuel tank, and not an inherent impossibility of the design. As the atmosphere of Mars is not dense, horizontal landing could be difficult (very high speed) and the wings might not be justified. As the parameters for the entry and descent to the surface would be constrained by the orbit of Phobos, increased range using aerodynamic surfaces could nevertheless be interesting, to increase the accessibility on the Martian surface. The Albatross concept is inspired from the X33 prototype, itself a precursor for the never built Venture Star. Thanks to a much lower Delta V requirements, less than 6 km/s for the ascent to Phobos, the technical design is simpler. The lifting body form facilitates the entry and descent stage, with lower temperatures at entry, a lower deceleration peak, and the possibility of landing at higher altitudes on the surface. This vehicle can also be considered as a more specifically adapted version of Starship, the vehicle SpaceX proposes for missions to Mars.

Albatross, a SSTO rocket

Capacities and design

Maintenance issues

The SSTO and any vehicle coming from Earth is likely to be unpiloted. Maintenance of the SSTO would therefore take place on Mars. As the martian landing and aerodynamic heating are likely to be the most difficult elements of the flight, this is a logical solution. Importantly, as the proposed SSTO would be relatively small and light (compared to Musk's Starship), it would be technically easy to transport it from the landing site to the maintenance site, and also to access to the highest parts of the vehicle for inspection.

The SSTO should have the capacity of bringing up maintenance personnel for the station if needed. The development of an unpiloted craft should also be easier and eliminate any risk to life. Automated orbital fuel transfer is also an essential part of many Mars exploration plans and should be old and reliable technology by the time a Phobos base is developed.

Industrial requirements for the construction on Mars

Industrial production on Mars using In-situ resources for common materials should be favorable economically when compared to transportation costs from Earth.

The Albatross SSTO requires 30 tonnes of fabrication materials, mostly structural elements, tanks and heat shields. Aluminum is readily available on Mars and could serve as the main structural and tank material for the Albatross. Iron, silica and carbonates are also available in the Martian lithosphere (for instance in hematite, sand, and sedimentary terrain), and specific separation methods could be developed. There will be a local market on Mars for structural materials used in habitats that can also serve for the Albatross. In fact, the construction of a number of Albatross type vehicles would be a small demand compared to infrastructure requirements for a Martian settlement. Martian habitats need to be built from pressure resistant materials due to the very low martian atmospheric pressure; the same industry might serve the launch vehicle market.

Large scale in-situ production of methane and oxygen propellants is required for most Mars missions planned today and the infrastructure for this would serve for the Albatross as well. Using methane rather than hydrogen also reduces the size of the tanks significantly, leaving more space for payloads and reducing structural requirements. The lower deltaV requirements simplifies the construction of the spacecraft on Mars when compared to Earth, as the volume of the tanks is not such a large fraction of the vehicle. A vehicle the size of the X-33 may then be able to carry a payload practically the size of the one planned for the larger Venture Star.

A petrochemical industry is possible on Mars, using methane from CO2 and precursor materials produced in the life support systems of a martian settlement. Products such as PVC, EPDM and other plastics can be produced using known methods and solar or nuclear energy.

Avionics and microprocessors might come from Earth in the early stages of production, until a local industry develops.

Relativity Space has proposed in its mission statement the use of 3D printing to produce a rocket on Mars using the technology it is developing for launchers on Earth [START_REF] Ellis | Relativity Space website[END_REF].

Conclusion

The Japanese Martian Moon eXploration mission (MMX) [START_REF] Usui | The Importance of Phobos Sample Return for Understanding the Mars-Moon System[END_REF], or the European candidate mission Phootprint [START_REF] Barraclough | Phootprint: A European Phobos Sample Return Mission[END_REF], planned for 2024, are sample return missions that should provide answers to the origin of the Martian moons. It should also determine if water and carbonates are present in significant quantities. If sufficient resources are identified, then Phobos might serve as an important propellant source for the development of transportation infrastructures between the Earth and Mars and between Mars and the asteroid belt. A locally produced SSTO could serve as an efficient space plane to travel between the surface and the Phobos spaceport hub. It could also be used to move resources from Mars to orbit, and could reach Earth GEO at a significantly lower deltaV cost than from the Earth's surface.

Fig. 1 .

 1 Fig.1. An Albatross SSTO moving in front of Phobos

  Albatross characteristics: • Size: Length: 21 meters; width: 25 meters • Main propulsion system: LCH4 and LO2, specific impulse 319s. • Total propellant: 240 tonnes. • Dry mass: 30 tonnes (including shielding for descent) • Maximum payload: 20 tonnes • DeltaV: 5500 m/s

Fig. 2 .

 2 Fig.2. Albatross model and Martian landing area

Table 1 .

 1 Phobos and Deimos orbital parameters.

		Deimos	Phobos
	Periapsis	23, 546	9, 234 km
	Apoapsis	23, 471	9, 518 km
	Orb. period	1.26 days	7 h 39 min
	Radius	6.2 km	11.3 km
	Surf. gravity	0.003 m/s2	0.0057 m/s2

Table 2 .

 2 DeltaV budget from Mars transfer to Phobos.

		To Phobos,	To Phobos w/o
		with	aerocapture, on
		aerocapture	average
	Mars transfer	962 m/s	2017 m/s

Table 3 .

 3 Delta V budget from Mars to geostationary orbit through Phobos.

	Mars to	Phobos to	Earth	Total
	Phobos	Earth	transfer
		transfer	(GTO) to
		(GTO)	GEO
	5.5 km/s	1.9 km/s	1.3 km/s	8.7 km/s
	Earth to LEO LEO to GEO		Total
	9.3 km/s	3.8 km/s		13.1
				km/s

nd International Astronautical Congress (IAC), Dubai, United Arab Emirates, 25-29 October 2021. Copyright 2021 by Jean-Marc Salotti. Published by the IAF, with permission and released to the IAF to publish in all forms.