Numerical investigation of a hypersonic flow around a capsule in CO2–N2 environment

European Journal of Mechanics / B Fluids 80 (2020) 146–156

Contents lists available at ScienceDirect

European Journal of Mechanics / B Fluids journal homepage: www.elsevier.com/locate/ejmflu

Numerical investigation of a hypersonic flow around a capsule in CO2 –N2 environment ∗

Zineddine Bouyahiaoui , Rabah Haoui, Abderrahmane Zidane Laboratoire de Mécanique Énergétique et Systèmes de Conversion, Université des Sciences et de la, Technologie Houari Boumediene, BP 32 El Alia, 16111 Bab Ezzouar, Algeria

article

info

Article history: Received 18 December 2018 Received in revised form 9 November 2019 Accepted 20 December 2019 Available online 24 December 2019 Keywords: Hypersonic flow Catalytic wall Nonequilibrium Viscous flow Finite volume method

a b s t r a c t A numerical investigation is conducted to study a viscous hypersonic flow around a capsule in chemical and vibrational nonequilibrium. A 16−reaction kinetic model was implemented in our twotemperature Hypersonic Reacting Flow (HRF) code. The vibrational relaxation times are evaluated using Millikan & White semi-empiric formula along with Park’s correction for high temperature. The utilized 16−reaction kinetic model demonstrates good modeling of the chemical nonequilibrium. The simulations are performed for three wall conditions: noncatalytic, equilibrium catalytic, and fully catalytic wall. Obtained results show a good agreement with the literature. The heat flux was found to vary from 10 to 41 W/cm2 depending on the wall condition. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction The conquest of planet Mars has become a major issue, especially during the last decade. Due to overcrowding and lack of resources, the need for vital space becomes a priority for the development and survival of humanity. Researchers began to think of a way to conquer a new planet that brings together livable conditions, namely water, air ... etc. As a result, several missions have been programmed to explore nearby planets like Venus and Mars, which turned out to be the perfect candidate. To this end, studies relative to the Mars journey had to be carried out. These studies were to focus on several aspects, such as launch, atmospheric reentry, and landing. Several launches were carried out with different probes such as Viking 1 and 2 (1975–1982) [1], Pathfinder (1996–1997) [2], MER A and B (2003–2010) [3], Phoenix (2008) [4], MSL (2012) [5] and Exomars (2016) [6]. The Mars atmosphere is composed of 95.7% carbon dioxide, 2.7% nitrogen, and 1.6% argon. It is characterized by its lowpressure density: near 7 Pa at high altitude and less than 600 Pa on the ground. This makes the braking of capsules weak and experimentation difficult and expensive. Therefore we resort to numerical calculations by developing CFD codes. Several numerical studies using different solvers were published in order to study the different aspects of the reentry, such as the study of chemical and vibrational non-equilibrium flow behind a normal shock wave [7], the flow around a blunt ∗ Corresponding author. E-mail address: [email protected] (Z. Bouyahiaoui). https://doi.org/10.1016/j.euromechflu.2019.12.009 0997-7546/© 2019 Elsevier Masson SAS. All rights reserved.

body for Martian reentry [8–12], aerothermal heating effects on capsules [13,14] and many other aspects involved in such flows. Regarding the numerical studies of aerodynamic and thermal effects during reentry, we can cite among them: Gnoffo et al. [15] investigated the influence of sonic-line location on Mars Pathfinder probe aerothermodynamics. They used two codes, the first code HALIS provides a perfect inviscid gas solution, and the second is a real gas, laminar and viscous code, LAURA [16]. The variation of angle-of-attack from 0◦ to 11◦ is studied to provide an understanding of the resulting aerodynamics and heating distribution changes. The work concluded that at small angle-of-attack (α ≤ 5◦ ), the sonic line location shifts from the shoulder to the nose cap and back again to the leeside symmetry plane because of the heat capacity ratio γ and cone half-angle of 70◦ . It is also showed that the sonic line movement affects the heating distribution. Mitcheltree et al. [17] proposed an aerothermal heating prediction for Mars Microprobe using the CFD code LAURA and experimental measurements. The heating rate at the stagnation point at 0 angle-of-attack (α = 0◦ ) was predicted to be 194 W/cm2 associated with a pressure of 0.064 atm. No significant heating augmentation due to radiation was predicted according to their model, and the forebody shock layer should remain laminar. The effect of angle-of-attack was also investigated and shows that, on the forebody frustum, the experimental data is as much as 30 percent lower than the CFD prediction. Edquist et al. [18] presented an aerothermodynamic design environment for a Mars Science Laboratory entry capsule heatshield. The design was based on Navier–Stokes flowfield simulation and entry trajectory from the 2009 launch. Baldwin–Lomax

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

algebraic turbulence model was used and agreed well with fully turbulent experimental heating data. In the design conditions, a correction bias was included to account for heatshield TPS surface roughness effects on turbulent shear stress and heat flux. At the maximum heat flux location, results showed that the total bias for the heat flux plus uncertainty is 50% more than the baseline heating. The peak design conditions found were: 197 W/cm2 heat flux, 471Pa shear stress, 0.371 atm pressure and 5477 J/cm2 total heat load. McDaniel et al. [19] worked on the aerothermal analysis of the Phoenix Martian entry vehicle using two CFD codes, namely DPLR and LAURA. Both codes solve the reacting Navier–Stokes equation on a structured finite volume grid for non-equilibrium gases. The difference between the two codes is the numerical method, DPLR [20] uses a modified Steger–Warming flux splitting for Inviscid fluxes with third-order accuracy using a monotone upstream centered scheme for conservation laws (MUSCL) extrapolation with a minmod flux limiter. LAURA uses Roe’s averaging for inviscid fluxes with second-order correction via Yee’s symmetric total variation diminishing scheme. Results showed that the flow over the Phoenix vehicle could be modeled as laminar. Also, the axisymmetric CFD solution was faster and cheaper to generate than the 3D solution with conservative results. Finally, the two codes yielded consistent results within 5% on the forebody, and the maximum heating rates occurred at the stagnation point and the shoulder. Hollis and Borrelli [21] made a review article comparing computational tools, ground and flight tests data in order to illustrate the encountered challenges in the numerical modeling of these phenomena and to provide test cases to evaluate the computational fluid dynamics codes predictions. Capsules used in this article were Mars Science Laboratory (MSL), Orion Crew Exploration Vehicle (CEV), and Fire II flight investigation of the reentry environment. This review shows the effectiveness of computations at predicting the effects of physical phenomena such as turbulence, chemical and vibrational non-equilibrium, rarefied flow, and radiation transport. It was found that there is still the need to validate computational models at non-equilibrium conditions for both Earth and Mars atmospheric environments. The accuracy of chemical and vibrational models seemed to be not fully validated, especially for the Martian CO2 environment. Bansal et al. [22] developed a new flow solver to simulate coupled hypersonic flow-radiation over a reentry vehicle using existing solvers and tools available in OpenFOAM. The hypersonicFoam solver developed combines features from two solvers, the first is rhoCentralFoam, and the second is reactingFoam. An excellent agreement was observed in flowfields obtained by the new solver. The heat transfer rates for the coupled case were observed to be significantly reduced compared to the uncoupled rates. Preci and Auweter-Kurtz [23] presented a numerical simulation using the URANUS non-equilibrium code. The flow field was calculated for several chemical and thermal models of CO2 –N2 inflow mixture. A sensitivity analysis of the radiative heat flux has also been performed by employing the PRADE radiation database and the HERTA algorithm. Results showed that the calculated flow in front of the vehicle was in good agreement with the results of other authors, contrary to those at the wall, which showed significant discrepancy due to simple wall boundary conditions. Armenise et al. [24] did a comparative study for Mars free stream conditions for two models: the state to state and multitemperature models. Calculations concluded that the multitemperature model achieved a better agreement for the heating predictions than the state-to-state. There were some discrepancies for the prediction of the species mole fractions, especially close to the wall. For a noncatalytic wall, the heating values

147

obtained for Mars Pathfinder with the multi-temperature model were in better agreement with postflight calculations. The multitemperature approach is much easier to implement and costs less compared to the state-to-state model. There are also recent works concerning Martian reentry flows, among them the works of Yatsukhno et al. [25],Hollis et al. [26] and [27], as well as other experimental works allowing the derivation of new thermochemical models such as [28–33]. All the studies cited before used different numerical approaches; however as stated by Hollis and Borrelli [21], the numerical models still need to be validated, especially for the Martian atmospheric environment. Thus, the present work is another contribution to the investigation of aerothermodynamics and thermal heating for Martian reentries of the MER capsule. It is focused on the flow simulation around a space capsule, taking into consideration three wall conditions in order to predict the heat flux range that the heat shield is exposed to during a reentry. 2. Mathematical model

Fig. 1. MER entry capsule configuration. Dimensions in meter.

The MER geometry is illustrated in Fig. 1. The outer mold line (OML) consists of a 70◦ sphere–cone heatshield and a 43◦ conic backshell [34]. The Navier–Stokes equations in a flux-vector formulation in a Cartesian coordinate system for a compressible viscous flow expressed in 2D-axisymmetric form is given as:

(

mes Ci,j

) ∂ Wi,j + ∂t

∑

(

− →

( ) − →) − .→ ηa −H .aire Ci,j = Ω

Fi,j i + Gi,j j

a∈{x,x′ ,y,y′ }

(1) mes(Ci,j ) is the measurement (in m3 ) of a infinitely small volume of center (i, j) and aire(Ci,j ) is the surface of the symmetry plane passing by the center of elementary volume with ηa the integrated normal. Where vectors W , F , G, H and Ω are given by:

⎛ ⎜ ⎜ ⎜ W =⎜ ⎜ ⎝

ρ ρu ρv ρe ρs ρm ev,m

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(2)

148

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

with ρ the mass density of the fluid, ρs and ρm are the partial densities of species s and molecules m. u and v are the components of the velocity vector.

⎛ ⎜ ⎜ ⎜ F =⎜ ⎜ ⎝ ⎛ ⎜ ⎜ ⎜ G=⎜ ⎜ ⎝ ⎛ ⎜ ⎜ ⎜ H=⎜ ⎜ ⎝

ρu ρ u2 + p − τxx ρ uv − τxy (ρ e + p) u − uτxx − vτxy + qx ρs u ρm ev,m u ρv ρv u − τxy ρv 2 + p − τyy (ρ e + p) v − uτxy − vτyy + qy ρs v ρm ev,m v

⎟ ⎟ ⎟ ⎟, ⎟ ⎠

(3)

⎞ ⎟ ⎟ ⎟ ⎟, ⎟ ⎠

(4)

⎞

−2τxy 2p − 2τyy −2uτxy − 2vτyy + 2vτzz + 2qy 0 0

⎟ ⎟ ⎟ ⎟, ⎟ ⎠

(5)

τm

ωc ,s ωv,m

⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(6)

e=

h0f , kJ/mol

CO2

1890 945 (2) 3360

−393.510

N2 O2 C2 CO NO CN O C N

3395 2239 2669 3074 2817 3074 – – –

0.0 0.0 830.457 −110.535 91.271 438.684 249.175 716.680 472.680

=

∑ Xs τm,s s

(11)

where Xs is the molar fraction of the species s, and τm,s is the characteristic time of vibration–translation exchanges of the molecule m with species s, defined by the semi-empirical formula of Millikan and White [36] as:

(

(

1

)

)

pτm,s = exp am,s T − 3 − bm,s − 18.42 [atm s]

The total energy per unit of mass e is written as: ns ∑

θv,s , K

The vibrational characteristic time τm is defined from the characteristic times related to the collisions between the molecule m and the species s: 1

⎞

0 0 0 0

⎜ ⎜ ⎜ Ω=⎜ ⎜ ⎝

Molecule

⎞

0

⎛

Table 1 Vibrational characteristic temperatures and standard formation enthalpies.

cv,s T +

nm ∑

s=1

Ym ev,m +

s=1

ns ∑ s=1

Ys h0f ,s

+

1( 2

u +v 2

2

)

(12)

With: 1

(7)

4

am,s = 0.00116µm2 ,s θv,3 m

(13)

1

h0f ,s is the formation enthalpy of the species s in J/kg given in Table 1, and cv,s is the specific heat at constant volume for each species. The pressure of the mixture is obtained by the equation of state: p

ρ

=

) ns ( ∑ R s=1

Ms

T = RT

(8)

where R is the gas constant, Ms the molar weight for species s, and R the specific gas constant for the mixture. The source term of vibration energy production ωv,m is calculated as follows:

ωv,m = ρm

ev,m (T ) − ev,m (Tv )

τm

+ ev,m (Tv ) ωc ,s

(9)

The first term of the right-hand side of Eq. (9) represents the Landau–Teller equation [35] for the vibrational–translational energy exchange mode. The second term is the loss/gain of vibrational energy due to chemical effects. This term is significant, but to avoid the use of complex models that are more computationally expensive, we opted for the mentioned above model along with Park’s CVD model for the influence of vibrational nonequilibrium on kinetics. ev,m (T ) is the equilibrium energy of vibration at the temperature of translation–rotation expressed as: R

ev,m (T ) = exp

(Mm

θv,m )

θv,m T

(10)

−1

θv,m is the characteristic temperature of vibration for any molecule given in Table 1. ev,m (Tv ) is the vibrational energy at the temperature of vibration Tv specified for each molecule.

bm,s = 0.015µm4 ,s

(14)

µm,s being the reduced molar weight in g/mole given as: µm,s =

Ms Mm

(15)

Ms + Mm

For temperatures superior to 8000 K, the equation gives relaxation times lower than those observed in experiments. To overcome this problem, C. Park [37] suggests the following relation to be added to the vibrational relaxation time:

τmp =

1

(16)

Nt cm σm

where cm is the molecular average velocity in m/s defined as follows:

√ cm =

8RT

(17)

π Mm

with σm being the effective collision cross-section to vibrational relaxation in m2 computed using the following,

σm = 10

−21

(

50000 T

)2 (18)

and Nt is the density of the number of collision particles of species s, where Na is the Avogadro number,

ρ Na

Nt = ∑ns

s=1

Xs Ms

(19)

Combining the two relations, the following expression for the vibrational relaxation time is obtained:

τm = τmMW + τmp

(20)

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

The source term of the conservation equation of chemical species ωc ,s is given by: r ∑ ) ( ′′ ωcs = Ms υs − υs′ Jr

(21)

r =1

Where Jr is the elementary sheet balance for a chemical reaction written as : ns ∑

kf

νs′ As ⇐⇒ kb

s=1

[ Jr =

kf

ns ∑

νs′′ As

∏ ( ρs )νs′ s

(22)

s=1

Ms

− kb

∏ ( ρs )νs′′ s

Ms

] (23)

νs′ and νs" are, respectively, the stoichiometric mole numbers of the reactants and products of species s for each chemical reaction r. kf and kb being the forward and the backward (or reverse) reaction rates. Bf

−

Θf

B

−

Θb

kf = Af Ta exp

kb = Ab Ta b exp

Ta

Ta

(24)

√

TTv

Table 2 Forward coefficients. N◦

Reaction

Af

Bf

Θf

1 2a 2b 3a 3b 4 5a 5b 6a 6b 7a 7b 8 9 10 11 12 13 14 15 16

C2 + M ⇐⇒ C + C + M N2 + A ⇐⇒ N + N + A N2 + Mol ⇐⇒ N + N + Mol O2 + A ⇐⇒ O + O + A O2 + Mol ⇐⇒ O + O + Mol CN + M ⇐⇒ C + N + M CO + A ⇐⇒ C + O + A CO + Mol ⇐⇒ C + O + Mol NO + X ⇐⇒ N + O + X NO + D ⇐⇒ N + O + D CO2 + A ⇐⇒ CO + O + A CO2 + Mol ⇐⇒ CO + O + Mol NO + O ⇐⇒ N + O2 N2 + O ⇐⇒ NO + N CO + O ⇐⇒ C + O2 CO + C ⇐⇒ C2 + O CO + N ⇐⇒ CN + O N2 + C ⇐⇒ CN + N CN + O ⇐⇒ NO + C CN + C ⇐⇒ C2 + N CO2 + O ⇐⇒ O2 + CO

3.7 × 10+14 3.0 × 10+22 7.0 × 10+21 1.0 × 10+22 2.0 × 10+21 2.5 × 10+14 3.4 × 10+20 2.3 × 10+20 1.1 × 10+17 5.0 × 10+15 1.4 × 10+22 6.9 × 10+21 2.8 × 10+09 2.0 × 10+12 3.9 × 10+13 2.0 × 10+17 1.0 × 10+14 1.1 × 10+14 1.6 × 10+13 5.0 × 10+13 2.1 × 10+13

0.00 −1.60 −1.60 −1.50 −1.50 0.00 −1.00 −1.00 0.00 0.00 −1.50 −1.50 1.00 0.50 −0.18 −1.00 0.00 −0.11 0.10 0.00 0.00

69900 113200 113200 59750 59750 71000 129000 129000 75500 75500 63275 63275 20000 38000 69200 58000 38600 23200 14600 13000 27800

(25)

The kinetic model used in the present work is a reduced version of the model found in [38], which is based on the Park Model [7]. It was obtained by eliminating the reactions that involve ionization and NCO. Constants Af , Ab , Bf , Bb , Θf and Θb depend on the considered reaction and are summarized in Tables 2 and 3. The reaction rates of CO2 dissociation and O2 taken from different references [39,40] are compared to those included in the present reduced model (Fig. 2). It is observed that all the reaction rates have relatively the same variations in a function of temperature. In order to assess the coherence of this model, it was tested for the prediction of equilibrium composition for a temperature range of [1000 K–4000 K] at two considered pressure values 1500 Pa and 15000 Pa. The obtained results were compared to those shown by Palmer and Cruden [41] and are represented in Figs. 3 and 4. The initial gas mixture used is composed of [96%CO2 − 4%N2 ]. We can see that the CO2 dissociates more quickly at lowpressure than at high pressure when the temperature increases. For temperatures below 2500 K, CO2 is the dominant species. Above this temperature, CO and O become the dominant species, this is the case for both values of pressure. The differences observed between our results and those of Palmer and Cruden are due to the fact that the latter use chemical equilibrium, whereas in our case, we use the nonequilibrium relaxation until the equilibrium is reached. In the case of hypersonic flows behind shock-waves, it is customary to express the reaction rate of dissociation reactions as a function of both Ttr and Tv,m [42]. We use the vibration– dissociation coupling because of the molecular vibration influence on the dissociation. Indeed, the more a molecule is excited in a vibrational plan, the more it breaks up easily. The coupling allows holding into account this influence. According to this, the average temperature of the dissociation reactions is equal to the geometric mean of the temperature of translation rotation and the vibrational temperature of each molecule: Ta =

149

(26)

The heat flux vector q is calculated for each direction x, y, and z. The heat flux for a direction j is constituted of three components: the first component is given by the Fourier law, the second term corresponds to the vibrational heat flux, and the third one

Table 3 Backward coefficients. N◦

Reaction

Ab

Bb

Θb

1 2a 2b 3a 3b 4 5a 5b 6a 6b 7a 7b 8 9 10 11 12 13 14 15 16

C2 + M ⇐⇒ C + C + M N2 + A ⇐⇒ N + N + A N2 + Mol ⇐⇒ N + N + Mol O2 + A ⇐⇒ O + O + A O2 + Mol ⇐⇒ O + O + Mol CN + M ⇐⇒ C + N + M CO + A ⇐⇒ C + O + A CO + Mol ⇐⇒ C + O + Mol NO + X ⇐⇒ N + O + X NO + D ⇐⇒ N + O + D CO2 + A ⇐⇒ CO + O + A CO2 + Mol ⇐⇒ CO + O + Mol NO + O ⇐⇒ N + O2 N2 + O ⇐⇒ NO + N CO + O ⇐⇒ C + O2 CO + C ⇐⇒ C2 + O CO + N ⇐⇒ CN + O N2 + C ⇐⇒ CN + N CN + O ⇐⇒ NO + C CN + C ⇐⇒ C2 + N CO2 + O ⇐⇒ O2 + CO

3.886 × 10+11 4.351 × 10+19 1.015 × 10+19 5.856 × 10+19 1.171 × 10+19 3.083 × 10+16 2.281 × 10+13 1.543 × 10+13 2.485 × 10+15 1.129 × 10+14 2.343 × 10+14 1.155 × 10+14 1.100 × 10+10 4.400 × 10+15 3.399 × 10+13 3.043 × 10+16 2.044 × 10+13 2.921 × 10+15 1.031 × 10+12 1.582 × 10+16 3.281 × 10+09

0.62 −1.24 −1.24 −1.19 −1.19 −0.50 0.37 0.37 0.27 0.27 0.02 0.02 1.00 0.50 −0.18 −1.00 0.00 −0.63 0.38 −0.59 0.78

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4000.0 0.0 14280.0 7351.0 5814.3 0.0 1327.1 0.0 23464.0

M: CO2 , CO, O2 , C2 , N2 , NO, CN, C, N, O Mol: CO2 , CO, O2 , C2 , N2 , NO, CN D: CO2 , NO, C, N, O X: CO, O2 , C2 , N2 , CN A: C, N, O

corresponds to the species diffusion. The heat flux in a direction j is computed as follows: qj = −keq

∑ ∂T ∂ Ys −ρ hs Ds ∂ xj ∂ xj s

(27)

keq denotes the coefficient of total thermal conductivity, i.e., it includes translational–rotational and vibrational thermal conductivities, which is a function of Prandtl number Pr and dynamic viscosity µ. Ds is the diffusion coefficient given in function of the Lewis number Le as: Ds = ks

Le

ρ cp s

(28)

Lewis number characterizes mass diffusion (ratio of energy transport due to diffusion and thermal conduction). In our case, it is

150

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

Fig. 2. Comparison of dissociation rates for CO2 (a) and O2 (b). Table 4 Blottner’s viscosity law coefficients.

Fig. 3. Equilibrium molar fractions comparison for p = 1500 Pa: solid lines is our model and dashed for Palmer and Cruden [41].

µ

µ

µ

Molecule

As

Bs

Cs

CO2 N2 O2 C2 CO NO CN O C N

−0.026654

1.107305 0.4329 −0.1455 0.8486858 0.887198 −0.0609 0.779370 0.4257 0.768602 0.5930

−14.291274 −11.8153 −8.9231 −13.182071 −13.269815 −9.4596 −12.825216 −11.5803 −12.956246 −12.3805

0.0203 0.0484 −0.011809 −0.014044 0.0452 −0.007770 0.0205 −0.007140 0.012

For a pure gas, the viscosity is given by the interpolation of Blottner et al. [44]:

[( ) ] µs = 0.1 exp Aµs ln (T ) + Bµs ln (T ) + Csµ µ

µ

(29)

µ

Coefficients As , Bs and Cs depend on the species considered and are given in Table 4. For air species, coefficients are valid for a temperature range from 1000 K to 30000 K as proposed by Gupta et al. [45]. For carbon species William [46] has built a viscosity law in the form of Blöttner’s model from the gases kinetic theory of Neufeld et al. [47] for a temperature range from 100 K to 20000 K. For multicomponent gas, such as a chemically reacting mixture, the mixture value of µmix must be calculated from the values of µs of each chemical species s by means of mixture rules. A common viscosity expression is Wilke rule, which states that:

µmix =

Xs µ s

∑ ∑

s′

s

(30)

Xs′ φss′

With: Fig. 4. Equilibrium molar fractions comparison for p = 15000 Pa: solid lines is our model and dashed for Palmer and Cruden [41].

fixed at 1.0 for Martian applications [43]. cp is the specific heat at constant volume for all energy modes.

1

φs,s′ = √

8

( 1+

Ms Ms′

) −21 (

( 1+

µs µs′

) 12 (

Ms′ Ms

) 21 ) (31)

Ms , Ms′ Xs and Xs′ are respectively the molar weights and the molar fractions of species s and s′ . For the thermal conductivity kmix of mixture, the same equation (30) can be used, replacing µ by k.

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

151

Fig. 5. Computational domain and boundary conditions. Fig. 6. Grid near the wall.

3. Numerical code

− →

The numerical study presented in this paper was performed using our in-house code HRF (Hypersonic Reacting Flow). This code was specifically designed for nonequilibrium hypersonic compressible flow studies and is based on the finite volume method in order to solve nonequilibrium Euler equations using a two-temperature model. The disturbance of domain technique [48] is used to transform the 3D model into a 2D-axisymmetric model, and the time integration is performed using an explicit iterative approach (Euler forward method). It was first used for the numerical simulation of nonequilibrium nozzle flow [42,49], and later, it was upgraded with Van-Leer’s flux vector splitting numerical scheme [50] in a first-order version, in order to be used for the simulation of hypersonic nonequilibrium compressible flows around a blunt body [51,52] or for atmospheric reentry of a space vehicle [12]. An alternative of this in-house code was written in order to solve nonequilibrium Navier–Stokes equations for the simulations of viscous hypersonic compressible flows. Details of the axisymmetric formulation, discretization, and Van-Leer’s flux vector splitting can be found in Ref. [53]. The boundary conditions used in simulations are illustrated by Fig. 5 and detailed in the following: 3.1. Freestream conditions

3.3. Axis of symmetry

Inlet velocity, pressure, and temperature are fixed since the flow is supersonic at the freestream. 3.2. Body surface Since the flow is considered viscous, the no-slip wall condition must be applied. Therefore, the velocity at the surface is set equal to zero (u = v = 0). Regarding temperature, it is maintained at a value of Tw = 1500 K. In this study, the wall shear stress is calculated by:

) ( Vt ∂ Vt =µ τw = µ ∂ n wall ∆n

Here we assume that the coordinate of the unit vector t is in − → the direction of the shear force at the wall and the unit vector n − → is normal to t Fig. 6. The total friction drag is caused by shear stress acting on the top and bottom surfaces. Concerning the catalysis model, the conditions applied at the wall must be specified in order to take into consideration these phenomena. Several definitions can be found in the literature. In this work, the catalysis modeling is based on the definition given by Anderson Jr [54]. The non-catalytic wall condition implies that there are no variations in the chemical composition near the wall, and thus we have: ∂∂Yns = 0 For a fully catalytic wall, all atoms recombine in order to form molecules. This condition is applied as follows: Yatoms = 0. Equilibrium catalytic wall as indicated by its designation supposes that local chemical equilibrium is reached near the wall at the given temperature and pressure, i.e. at Tw = 1500 and a pressure of pw ≈ 9000 Pa. For these conditions, the composition at equilibrium is the same as the freestream. Therefore the condition will be: Ys = Ys (freestream).

The conditions at the axis of symmetry boundary are no flow and no scalar flux across the boundary.

3.4. Outlet boundary conditions At the exit of the computational domain, the values of the flow parameters are extrapolated from the interior values, including in the boundary layer.

(32)

where:

− →− →

Vt = V . t

(33)

and:

→ ∆n = (∆x, ∆y) .− n

(34)

with:

∆x = x (i, j + 1) − x (i, j)

(35)

and:

∆y = y (i, j + 1) − y (i, j)

(36)

Fig. 7. Grid and computational domain (60 × 197).

152

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156 Table 5 MER free stream conditions.

Fig. 8. Mach number over MER body.

Fig. 9. Flow Streamlines.

Name

Symbol & unit

Value

Inlet velocity Pressure Temperature Wall temperature CO2 Mass fraction N2 Mass fraction

V [m/s] p [Pa] T [K] Tw [K] YCO2 YN2

5223 7.87 140 1500 0.97 0.03

4.2. Viscous real gas simulations of the front domain In order to reduce the computation time, a reduced computational domain representing the fore-body is used for the viscous non-equilibrium case. Simulations are run using the same freestream conditions as the frozen inviscid simulation. Figs. 11 and 12 represent the tangential velocity profile and tangential stress normal to the wall at x = 0.2 m respectively. Velocity goes from zero at the wall to a value of around 1000 m/s in the boundary layer with a thickness of 20 mm. Then, it goes from 1000 m/s to 1800 m/s in the flow outside of the boundary layer until reaching the shock at 80 mm from the wall. From Fig. 11, we also notice that the (80 × 808) grid (80 in the normal direction and 808 in the tangential direction) gives an almost straight shock which is expected from the use of a first-order FVS type numerical scheme. Hence, this mesh will be retained to conduct the simulations. Tangential stress shows the exact positions of the shock, the flow behind, and the boundary layer. Note that the stress decreases in the boundary layer because of the fast decrease in temperature, which is a consequence of the cold wall temperature Tw = 1500 K.

Fig. 10. Temperature over MER body.

4. Results and interpretations 4.1. Inviscid perfect gas simulation of the whole domain (front and wake) A first simulation of the flow over the whole MER capsule is conducted using our inviscid solver with non-reacting gas. The free stream conditions are listed in Table 5. The numerical grid used for our calculations is illustrated in Fig. 7, with 60 cells normal to the wall and 197 streamwise. The chosen convergence criteria for our simulations is a maximum residue of 10−5 . The calculations converged after 131098 iterations equivalent to an execution time of 7328 s. From Figs. 8 and 9, we observe that the detached shock wave is located at approximately 0.25 m in front of the body. For the after-body area, the gas cools down at the shoulder. Though, hot gases are raised to the payload by the recirculation zone developed in the after-body region (Fig. 10). The temperature distribution (Fig. 10) for the fore-body is characterized by a sudden increase in temperature to a value of 20715 K after the shock due to the flow compression. As the flow goes around the body surface, it accelerates, this causes a temperature diminution. For the after-body, flow velocity decreases with the formation of a recirculation region, which will cause an increase in temperature to a maximum value of 24457 K.

Fig. 11. Tangential velocity profile in boundary layer for chemically reacting viscous flow.

Fig. 12. Tangential stress.

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

153

Fig. 14. Temperature contour for viscous flow.

Fig. 13. Mach number contour for viscous flow.

Mach number and temperature distributions are illustrated in Figs. 13 and 14 respectively. The sonic line represented by a discontinued line, shifts from the shoulder to the nose cap and back again on the leeside as reported by Gnoffo et al. [15]. We also notice that the shock position is much closer to the body than in the frozen flow case. This is due to the chemical and vibrational processes. On the other hand, the temperature reaches a peak value of almost 16000 K. A sudden decrease after the peak is caused by the dissociation reactions which are endothermic and the vibrational relaxation. After a short distance, the temperature reaches the wall value. Fig. 15 represents kinetic and vibrational temperatures of each molecule along the stagnation line. As stated before, the temperature reaches a peak behind the shock front, which induces CO2 –N2 dissociation. Other species start appearing namely CO, O2 , NO, CN, C2 , C, N and O as illustrated in Figs. 17–19 for the noncatalytic, equilibrium catalytic and fully catalytic wall respectively. Vibrational temperatures rise just behind the shock due to the excitation of the vibrational modes, then the vibrational temperatures meet the trans-rotational temperature, showing that thermal equilibrium is reached. Thus, the region from x ≈ −0.022 m to x ≈ −0.005 m represents the thermal equilibrium region. Beyond this region, all temperatures decrease to reach wall temperature. Fig. 16 illustrates a comparison of temperature profiles along the stagnation line for the same freestream conditions. The first thing to be noticed is the differences in the temperature peak after the shock, i.e. for our results and those of Preci and AuweterKurtz [23] the peak is about 16000 K, while for Rouzaud et al. [55],Druguet [56] and Potter et al. [57] it is about 7500 K. This is due to the consideration of the vibrational nonequilibrium in our calculations an those of Preci and Auweter-Kurtz. This is reasonable since, in the case of nonequilibrium, hundred of collisions are required to allow the translational degree of freedom to pass energy to the vibrational degree. After that our temperature decrease rate is different from others, we suspect that this is due to the difference in the employed chemical and vibrational nonequilibrium models. Depending on the wall condition (Fig. 17 to 19), the mass fraction on the boundary layer differs between the three figures. For

Fig. 15. Temperature profiles along the stagnation line.

Fig. 16. Comparison between temperature variations along the stagnation line for the same freestream condition.

the noncatalytic wall, there is no variation of the mass fraction at the wall. In contrary equilibrium and fully catalytic wall depend on the composition at the wall which will induce the temperature distribution in the boundary layer as seen in Fig. 23. For Figs. 17 and 18, the dashed line represent CEA [58] values for the same conditions at the wall. Pressure, heat flux and shear stress over the body are illustrated in Figs. 20–22 respectively for the three wall conditions. Three regions can be distinguished for all figures. The first region is around the nose cap, where the pressure reaches a maximum of almost 9000 Pa and diminishes to reach a constant value. For heat

154

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

Fig. 17. Mass fraction for a noncatalytic wall along the stagnation line: present work (solid lines), CEA (dashed lines).

Fig. 20. Pressure over the body surface.

Fig. 21. Heat flux over the body surface. Fig. 18. Mass fraction for an equilibrium catalytic wall along the stagnation line: present work (solid lines), CEA (dashed lines).

Fig. 22. Shear stress over the body surface. Fig. 19. Mass fraction for a fully catalytic wall along the stagnation line.

flux, the maximum value is around 41 W/cm2 for an equilibrium catalytic wall, 17 W/cm2 for a fully catalytic wall and 10 W/cm2 for a noncatalytic wall. For the latter, the heat flux maximum value is on the same order of magnitude as with those found in literature, such as Rouzaud et al. [55] 12 W/cm2 and Druguet [56] 14 W/cm2 . The differences noticed for the heat flux are due to species conditions imposed at the wall 5. Conclusion A numerical model was developed to investigate the viscous hypersonic flow around a capsule. Chemical and vibrational

nonequilibrium were considered. Therefore, a 16-reactions kinetic model based on the Park model was implemented in our in-house code HRF. Vibrational relaxation times were evaluated using Millikan and White formula with Park correction. The flow was considered laminar, and Blottner’s model was used for viscosity computation. In order to reduce the computation time, an inviscid axisymmetric simulation of the hypersonic flow around the whole capsule was first simulated. Results from this first simulation were used to locate the shock position in order to reduce the computational domain for the simulation of the viscous hypersonic flow in chemical and vibrational nonequilibrium. The modeling of the chemical and vibrational nonequilibrium was applied to a 2D axisymmetric simulation of the ballistic

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

Fig. 23. Comparison of Temperature distributions at the stagnation line for the different wall conditions.

entry of the MER front part capsule. Three wall conditions were considered: a noncatalytic wall, an equilibrium catalytic wall, and a fully catalytic wall. We can summarize the significant conclusion of this study in two points: The obtained results showed a good agreement with literature for the aerothermodynamics parameters and the heat flux estimation should be between 10 W/cm2 and 41 W/cm2 which are the bounds to be used for heatshield design to ensure optimum and secure design of the capsule. The accuracy of these results can be improved by using a higher-order numerical scheme and the inclusion of turbulence phenomena. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. References [1] R. Ingoldby, F. Michel, T. Flaherty, M. Doryand, B. Preston, K. Villyard, R. Steele, Entry data analysis for viking landers 1 and 2 final report, NASA CR-159388, 1976. [2] D.A. Spencer, R.C. Blanchard, R.D. Braun, P.H. Kallemeyn, S.W. Thurman, Mars pathfinder entry, descent, and landing reconstruction, J. Spacecr. Rockets 36 (3) (1999) 357–366. [3] R. Roncoli, J. Ludwinski, Mission design overview for the mars exploration rover mission, in: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2002, p. 4823. [4] M.D. Garcia, K.K. Fujii, Mission design overview for the phoenix mars scout mission, 2007. [5] R. Prakash, P.D. Burkhart, A. Chen, K.A. Comeaux, C.S. Guernsey, D.M. Kipp, L.V. Lorenzoni, G.F. Mendeck, R.W. Powell, T.P. Rivellini, et al., Mars science laboratory entry, descent, and landing system overview, in: 2008 IEEE Aerospace Conference, IEEE, 2008, pp. 1–18. [6] T. Tolker-Nielsen, EXOMARS 2016-schiaparelli anomaly inquiry, 2017. [7] C. Park, J.T. Howe, R.L. Jaffe, G.V. Candler, Review of chemical-kinetic problems of future NASA missions. II-Mars entries, J. Thermophys. Heat transfer 8 (1) (1994) 9–23. [8] R. Paciorri, F. Sabetta, B. Favini, On the role of vibrational excitation in hypersonic flow computations, Meccanica 33 (4) (1998) 331–347. [9] S. Surzhikov, Radiation gas dynamics of Martian space vehicles, in: Doklady Physics, vol. 57, Springer, 2012, pp. 119–124. [10] S. Surzhikov, V. Borovoi, A. Skuratov, Study of convective heating of segmental-conical Martian descent vehicle in shock wind tunnel, in: 34th AIAA Fluid Dynamics Conference and Exhibit, 2004, p. 2634. [11] O. Rovenskaya, V. Aristov, Kinetic simulation of a supersonic compressible flow over different geometry bodies, Eur. J. Mech. B Fluids 64 (2017) 102–111. [12] Z. Bouyahiaoui, R. Haoui, Hypersonic flow of CO2-N2 mixture around a spacecraft during the atmospheric reentry, World Acad. Sci. Eng. Technol. Int. J. Mech. Aerosp. Ind. Mechatronic Manuf. Eng. 11 (12) (2017) 1865–1868.

155

[13] F.S. Milos, Y.-K. Chen, W.M. Congdon, J.M. Thornton, Mars pathfinder entry temperature data, aerothermal heating, and heatshield material response, J. Spacecr. Rockets 36 (3) (1999) 380–391. [14] X. Yang, Y. Gui, W. Tang, Y. Du, L. Liu, G. Xiao, D. Wei, Surface thermochemical effects on TPS-coupled aerothermodynamics in hypersonic martian gas flow, Acta Astronaut. 147 (2018) 445–453. [15] P.A. Gnoffo, K.J. Weilmuenster, R.D. Braun, C.I. Cruz, Influence of sonic-line location on mars pathfinder probe aerothermodynamics, J. Spacecr. Rockets 33 (2) (1996) 169–177. [16] P.A. Gnoffo, F.M. Cheatwood, User’s manual for the Langley aerothermodynamic upwind relaxation algorithm (LAURA), 1996. [17] R.A. Mitcheltree, M. Difulvio, T. Horvath, Aerothermal heating predictions for mars microprobe, 1998. [18] K. Edquist, A. Dyakonov, M. Wright, C. Tang, Aerothermodynamic design of the mars science laboratory heatshield, in: 41st AIAA Thermophysics Conference, 2009, p. 4075. [19] R.D. McDaniel, M.J. Wright, J.T. Songer, Aeroheating predictions for Phoenix entry vehicle, J. Spacecr. Rockets 48 (5) (2011) 727–745. [20] M.J. Wright, G.V. Candler, D. Bose, Data-parallel line relaxation method for the Navier-Stokes equations, AIAA J. 36 (9) (1998) 1603–1609. [21] B.R. Hollis, S. Borrelli, Aerothermodynamics of blunt body entry vehicles, Prog. Aerosp. Sci. 48 (2012) 42–56. [22] A. Bansal, A. Feldick, M. Modest, Simulation of hypersonic flow and radiation over a mars reentry vehicle using openfoam, in: 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 2012, p. 650. [23] A. Preci, M. Auweter-Kurtz, Sensitivity analysis of non-equilibrium martian entry flow to chemical and thermal modelling, in: 53rd AIAA Aerospace Sciences Meeting, 2015, p. 0982. [24] I. Armenise, P. Reynier, E. Kustova, Advanced models for vibrational and chemical kinetics applied to Mars entry aerothermodynamics, J. Thermophys. Heat Transfer 30 (4) (2015) 705–720. [25] D. Yatsukhno, S. Surzhikov, D. Andrienko, J. Annaloro, P. Omaly, Different estimations of the convective and radiative heating for the martian entry probes, in: AIAA Scitech 2019 Forum, 2019, p. 0793. [26] B.R. Hollis, M. Barnhardt, M.G. MacLean, A.T. Dufrene, T.P. Wadhams, Turbulent aeroheating measurements on a 7-deg half-angle sphere-cone in a high-enthalpy CO2 expansion tunnel, in: 2018 AIAA Aerospace Sciences Meeting, 2018, p. 0377. [27] X. Lin, L. Chen, J. Li, F. Li, X. Yu, Experimental and numerical study of carbon-dioxide dissociation for mars atmospheric entry, J. Thermophys. Heat Transfer 32 (2) (2017) 503–513. [28] D. Potter, T. Eichmann, A. Brandis, R. Morgan, P. Jacobs, T. McIntyre, Simulation of radiating CO2-N2 shock layer experiments at hyperbolic entry conditions, in: 40th Thermophysics Conference, 2008, p. 3933. [29] S. Nishimura, A. Lemal, H. Takayanagi, S. Nomura, S. Matsuyama, K. Fujita, Analysis of CO2 plasma infrared radiation measurements, in: 54th AIAA Aerospace Sciences Meeting, 2016, p. 0739. [30] H. Takayanagi, S. Nomura, A. Lemal, K. Fujita, Measurements of nonequilibrium carbon dioxide infrared radiation in an expansion tube, in: AIAA paper 1369, 2017, p. 2017. [31] B.A. Cruden, A.M. Brandis, M.E. MacDonald, Characterization of CO thermochemistry in incident shockwaves, in: 2018 Joint Thermophysics and Heat Transfer Conference, 2018, p. 3768. [32] A. Brandis, C. Johnston, B. Cruden, D. Prabhu, A. Wray, Y. Liu, D. Schwenke, D. Bose, Validation of CO 4th positive radiation for Mars entry, J. Quant. Spectrosc. Radiat. Transfer 121 (2013) 91–104. [33] A. Lemal, S. Nishimura, S. Nomura, H. Takayanagi, S. Matsuyama, K. Fujita, Analysis of VUV radiation measurements from high temperature air mixtures, in: 54th AIAA Aerospace Sciences Meeting, 2016, p. 0740. [34] M. Schoenenberger, F.M. Cheatwood, P. Desai, Static aerodynamics of the Mars Exploration Rover entry capsule, in: 43rd AIAA Aerospace Sciences Meeting and Exhibit, 2005, p. 56. [35] L. Landau, E. Teller, Contribution to the theory of sound dispersion, Phys. Z. 10 (1936) 34. [36] R.C. Millikan, D.R. White, Systematics of vibrational relaxation, J. Chem. Phys. 39 (12) (1963) 3209–3213. [37] C. Park, Assessment of two-temperature kinetic model for ionizing air, AIAA Paper 87-1574, Technical Report, AIAA, Reston, VA, USA, 1987. [38] A.S. Dikalyuk, P. Kozlov, Y. Romanenko, O. Shatalov, S. Surzhikov, Nonequilibrium spectral radiation behind the shock waves in Martian and Earth atmospheres, in: 44th AIAA Thermophysics Conference (San Diego, USA, 2013), 2013, pp. 1–27. [39] K. Fujita, S. Matsuyama, T. Suzuki, Prediction of forebody and aftbody heat transfer rate for Mars aerocapture demonstrator, in: 43rd AIAA Thermophysics Conference, 2012, p. 3001. [40] B.A. Cruden, A.M. Brandis, D.K. Prabhu, Compositional dependence of radiance in CO2/N2/Ar systems, in: 44th AIAA Thermophysics Conference, 2013, p. 2502. [41] G. Palmer, B. Cruden, Experimental validation of CO2 radiation simulations, in: 43rd AIAA Thermophysics Conference, 2012, p. 3188.

156

Z. Bouyahiaoui, R. Haoui and A. Zidane / European Journal of Mechanics / B Fluids 80 (2020) 146–156

[42] A. Zidane, R. Haoui, M. Sellam, Z. Bouyahiaoui, Numerical study of a nonequilibrium H2-O2 rocket nozzle flow, Int. J. Hydrogen Energy 44 (8) (2019) 4361–4373. [43] J. Charbonnier, J. Couzi, W. Dieudonné, J. Vérant, Definition of an axially symmetric testcase for high temperature gas radiation prediction in mars atmosphere entry, Rapport CNES NG104-07-TF-001-CNES, 2003. [44] F.G. Blottner, M. Johnson, M. Ellis, Chemically reacting viscous flow program for multi-component gas mixtures, Technical Report, Sandia Labs., Albuquerque, N. Mex., 1971. [45] R.N. Gupta, J.M. Yos, R.A. Thompson, K.-P. Lee, A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K, 1990. [46] J. William, Etude des processus physico-chimiques dans les écoulements détendus à haute enthalpie: application à la soufflerie à arc F4 (Ph.D. thesis), Aix-Marseille 1, 1999. [47] P.D. Neufeld, A. Janzen, R. Aziz, Empirical equations to calculate 16 of the transport collision integrals Ω (l, s)* for the Lennard-Jones (12–6) potential, J. Chem. Phys. 57 (3) (1972) 1100–1102. [48] J. Goudjo, Désidéri,‘‘A finite volume scheme to resolution an axisymmetric Euler equations (Un schéma de volumes finis décentré pour la résolution des équations d’Euler en axisymétrique)’’, Research report INRIA 1005, 1989. [49] R. Haoui, A. Gahmousse, D. Zeitoun, Écoulement hors d’équilibre chimique et vibrationnel dans une tuyère hypersonique axisymétrique, Int. J. Therm. Sci. 40 (8) (2001) 787–795.

[50] B. Van Leer, Flux-vector splitting for the Euler equations, in: Eighth International Conference on Numerical Methods in Fluid Dynamics, Springer, 1982, pp. 507–512. [51] R. Haoui, Finite volume analysis of a supersonic non-equilibrium flow around an axisymmetric blunt body, Int. J. Aeronaut. Space Sci. 11 (2) (2010) 59–68. [52] R. Haoui, Physico-chemical state of the air at the stagnation point during the atmospheric reentry of a spacecraft, Acta Astronaut. 68 (11) (2011) 1660–1668. [53] R. Haoui, Effect of mesh size on the supersonic viscous flow parameters around an axisymmetric blunt body, Int. J. Mech. Aerosp. Ind. Eng. 8 (7) (2014). [54] J.D. Anderson Jr, Hypersonic and High-temperature Gas Dynamics, American Institute of Aeronautics and Astronautics, 2006. [55] O. Rouzaud, L. Tessé, T. Soubrié, A. Soufiani, P. Rivière, D. Zeitoun, Influence of radiative heating on a martian orbiter, J. Thermophys. Heat Transfer 22 (1) (2008) 10–19. [56] M.-C. Druguet, Prediction of the flow field over an orbiter entering the mars atmosphere, Shock Waves 20 (3) (2010) 251–261. [57] D. Potter, S. Karl, M. Lambert, Simulation of a Mars aerocapture vehicle on flow adapted unstructured grids with the DLR TAU code, in: Proceedings of the 5th International Workshop on Radiation of High Temperature Gases in Atmospheric Entry Barcelona, Spain. European Space Agency, Noordwijk, the Netherlands, 2012. [58] B. McBride, S. Gordon, Computer program for calculation of complex chemical equilibrium compositions and applications, 1996.