During hypersonic reentry, a spacecraft experiences several different fluid flow regimes, which usually require the application of different software frameworks to simulate the respective regimes. This study aims to evaluate the hyStrath library for predicting aerodynamic lift and drag coefficients of complex three-dimensional (3D) geometries during hypersonic Mars reentry, using flight data of the Viking 1 mission as reference. A range of altitudes ( h = 30 140 km) and Mach numbers ( M = 13 24) where flight data is available is considered, covering the rarefied, transitional, and continuum fluid flow regimes. The hyStrath library contains a set of modified solvers and state-of-the-art thermophysical and chemistry models within the framework of OpenFOAM, dedicated to modeling high-enthalpy hypersonic flow problems. Depending on the flow regime, the computational fluid dynamics solver hy2Foam or direct-simulation Monte Carlo solver dsmcFoam+ are employed in the study. Because hyStrath is based on OpenFOAM, it allows the use of an unstructured adaptive mesh refinement approach for arbitrary geometries. We obtain excellent results throughout all investigated flow regimes and Mach numbers with an average deviation of 1.5% and 2% from the measured lift and drag coefficients, respectively. The applicability of the framework for accurately modeling both rarefied and continuum Mars reentry problems of complex 3D geometries such as the Viking capsule is demonstrated.

Accurately predicting aerothermal loads on reentering spacecraft will be critical for the design of future Mars mission vehicles and to ensure a safe descent during atmospheric reentry. The steadily increasing interest of (eventually manned) Mars missions such as SpaceX's Starship require ever more complicated designs and reentry approaches, making reliable modeling indispensable.

Compared to the Earth, the significantly thinner atmosphere of Mars results in the spacecraft experiencing rarefied flow until lower altitudes. However, during reentry, the majority of deceleration and aerothermal loads take place in the high-enthalpy hypersonic continuum flow regime, invoking several physical phenomena that occur around the reentering vehicle. These include the formation of a bow shock followed by a zone of high temperature and pressure in front of the heat shield. Here, the fluid is characterized by strong thermal and chemical non-equilibrium, yielding dissociation of molecular species and excitation of electronic and vibrational energy modes. Experimental testing of subscale models in wind tunnels with these conditions requires huge amounts of effort. For this reason, a large proportion of research in the hypersonic area has been devoted to numerical methodology in order to reliably predict the flow physics of reentering spacecraft.

The most widely used approach to modeling rarefied fluid flow ( Kn > 1) is the DSMC method proposed by Bird.1,2 Significant extensions and improvements have since been achieved in terms of modeling chemistry,3 multiple energy modes,4,5 and their relaxation6,7 as well as surface-gas interaction models.8 Commonly employed DSMC codes are SPARTA,9 DAC,10 MONACO,11 DS3V,12 SMILE,13 and dsmcFoam.14 Aerothermal loads during rarefied hypersonic Mars reentry of the MSL mission were extensively investigated by Borner et al.15 Similarily, DSMC simulations concerned with both transitional and rarefied aerodynamic analysis were carried out by Moss et al.16,17 for the Pathfinder and by Blanchard et al.18 for the Viking missions.

The gradual increase in density through loss of altitude results in the fluid flow eventually becoming continuous ( Kn < 0.01). Compared to the scope of research published in the area of rarefied Mars reentry, more research has been devoted to modeling the continuum flow regime. NASA's compressible code LAURA19 was used to study the aerodynamics of the Viking capsule20 as well as its backshell heating.21 Similar analysis using the DLR TAU code22 focused on radiative and convective Viking afterbody heating.23 In a publication by Egorov and Pugach,24 the HSFlow solver was used to model the flow physics of the ExoMars vehicle traveling at Mach 29 and 69 km altitude. Another study employing the LAURA and DPLR25 codes was dedicated to the investigation of brownout and blackout analysis during Mars reentry of the Mars 2020 mission.26 Using a customized version of the commercially available code Ansys Fluent, non-equilibrium flowfield analysis of a manned braking system (MBS) lifting body for manned Mars exploration missions was investigated by Viviani et al.27,28 Detailed studies of continuum surface-gas reactions during Mars reentry were reported from Refs. 29–31 and a technical report summarizing and evaluating different models for modeling Mars reentry using computational fluid dynamics (CFD) including recommendations for computing chemistry, transport properties, and surface–gas interaction were concluded by Noeding.32 

In a recent development by Zeng et al.,33 nonlinear coupled constitutive relations (NCCR) were coupled with the chemical models from Ref. 34 and the two-temperature model from7 for modeling hypersonic Earth reentry in the continuum and transitional flow regime. Compared with existing Navier–Stokes and DSMC results, they were able to obtain improved results for characteristic quantities such as the electron number density or heat flux coefficient based on flight data of Earth reentry experiments such as the RAM C-II test.35 They attribute the better agreement with the improved capturing of extreme (rarefied) non-equilibrium effects when using the NCCR approach. In a follow-up study,36 they investigated jet interaction with high Mach number freestream flow including real gas effects.

In the present study, hypersonic reentry of the Viking 1 lander is simulated for altitudes spanning from 140 km (rarefied, Ma = 25) down to 30 km (continuum, Ma = 13) using the hyStrath library37,38 that contains DSMC and CFD solvers that are built upon existing OpenFOAM39 solvers, namely, rhoCentralFoam and dsmcFoam. They were extended with thermodynamic and chemical non-equilibrium models by Casseau et al.37,38 so that they are able to better capture the physical phenomena of high-enthalpy hypersonic flows. The goal of the present research is to evaluate the dsmcFoam+ and hy2Foam solvers for predicting aerodynamic coefficients of Mars reentry vehicles traveling in the rarefied and continuum hypersonic flow regime at non-zero angle-of-attack. The Viking lander was selected for this purpose due to the good availability of geometric and aerodynamic data. During reentry, various on-board sensors acquired data of the atmospheric properties and spatial accelerations, from which the lift and drag coefficients were derived.40 

The Viking 1 lander capsule successfully landed on the surface of Mars on July 20, 1976 after a 10 min reentry journey. The geometry of the capsule and telemetry of the velocity and freestream temperature is shown in Fig. 1 from the point of reentry at 140 km down to 30 km with the corresponding flow regimes denoted on the right side.

FIG. 1.

Geometry of the Viking lander (left, units in m) and derived velocity, freestream temperature, and (right) flow regimes determined by Blanchard and Walberg40 based on freestream Knudsen number defined in Eq. (1).

FIG. 1.

Geometry of the Viking lander (left, units in m) and derived velocity, freestream temperature, and (right) flow regimes determined by Blanchard and Walberg40 based on freestream Knudsen number defined in Eq. (1).

Close modal

In the present study, the investigated range of reentry cases was divided into steps of 10 km in the rarefied and transitional regimes and 5 km in the continuum regime, yielding a total of 14 performed simulations. The freestream properties and telemetry of the capsule for each altitude were taken from analyses carried out by Blanchard and Walberg40 and Nier et al.41 and set accordingly at the numerical boundary conditions, which are illustrated in Fig. 2. A table detailing freestream properties (species mass fraction, pressure and temperature) and telemetry (angle-of-attack and velocity) for each simulated altitude can be found in  Appendix A of Table I.

FIG. 2.

Viking capsule geometry residing within farfield and outlet boundary conditions. Global coordinate system is located in the center of the heatshield, from where the capsule is rotated by the respective angle-of-attack.

FIG. 2.

Viking capsule geometry residing within farfield and outlet boundary conditions. Global coordinate system is located in the center of the heatshield, from where the capsule is rotated by the respective angle-of-attack.

Close modal
The appropriate method for simulating each respective altitude is selected based on the freestream Knudsen number that determines the flow regime, defined as
Kn = λ L = k B T 2 π d ¯ 2 p L Flow regime : Kn < 0.01 Continuum hy 2 Foam ( CFD ) 0.01 < Kn < 1 Transitional > 1 < Kn Rarefied dsmcFoam + ( DSMC ) .
(1)

In Eq. (1), λ denotes the particle mean free path, kB is the Boltzmann constant, T is the freestream temperature, d ¯ is the mean particle kinetic diameter of the mixture, p is the freestream pressure, and L is the characteristic length, i.e., the capsule diameter.

Both the transitional and rarefied regime ( Kn > 0.01) were simulated using the dsmcFoam+ solver which is an extension of the existing dsmcFoam solver by,14 based on the DSMC method originally proposed by Bird.1 Key features required for modeling hypersonic flow problems were implemented by White et al.42 including (molecular) electro-vibrational energy modes and chemical reactions. One of the main benefits is the design within the OpenFOAM framework, allowing the utilization of fully unstructured meshes for arbitrary geometries. This aspect played a crucial role in our study since the wake of the spacecraft produces a very-low-density rarefied flow zone that would lead to numerical and unphysical instabilities when employing the continuum solver despite Kn 0.1. With an unstructured grid, we could substantially increase the cell volume in the wake of the Viking capsule while decreasing the cell volume in the front zone affected by the bow shock, ensuring both sufficient number of particles per cell and resolution of the mean free path.

DSMC simulations were carried out including 9 species: CO2, N2, Ar, O2, CO, O, N, NO, and C. 16 reactions were considered using the quantum-kinetic (Q-K) chemical reaction model according to Ref. 43 and implemented into dsmcFoam by Scanlon et al.44 A table containing all chemical reactions formulas and types is provided in  Appendix B (see Table II).

Initially, the simulations were run until quasi-steady-state was reached, i.e., N in N out. From there, field averaging was activated until the maximum global cell-wise rate of change was below 0.1%. The time step Δ t and DSMC weighting factor w were set in a manner so that the mean-collision time t mc to time step ratio and mean free path to cell size ratio remained below 1, respectively
t mc Δ t < 1 and λ V i 1 / 3 < 1 .
(2)

The computational meshes were generated with the snappyHexMesh mesher available in the OpenFOAM framework. We implemented a semi-automated procedure to adaptively refine the mesh based on the output files written by the solver since the dsmcFoam+ implementation does not allow for run-time adaptive mesh refinement (AMR). In our approach, the macro field of the Mach and density gradient were evaluated after quasi-steady-state was achieved and logarithmic 3D contour surfaces exported using the post-processing tool ParaView. The 3 D surfaces were used as refinement regions within snappyHexMesh and the DSMC simulation restarted. This algorithm was repeated until the mean-free particle path was adequately resolved in all cells. An overview of the computational meshing approach is shown in Fig. 3, outlining colored boundary conditions and pre/post adaptive mesh refinement section cuts.

FIG. 3.

(a) Baseline mesh used for DSMC simulations and (b) example of adaptively refined mesh for 100 km case.

FIG. 3.

(a) Baseline mesh used for DSMC simulations and (b) example of adaptively refined mesh for 100 km case.

Close modal

Freestream properties were enforced at the farfield and outlet boundaries whereby each crossing particle was deleted. The Variable Soft Sphere (VSS) model introduced by Koura and Matsumoto45 was utilized to model binary collision events, using the species-dependent collision parameters from the study of Ref. 46 (for species-wise VSS properties, see  Appendix C of Table III). Particle-wall interactions were modeled with the diffuse-specular wall boundary condition provided in dsmcFoam+, using a diffuse fraction of 0.9 for rarefied flow cases ( Kn > 1) and 1.0 for transitional flow cases ( 0.1 < Kn < 1) as estimated by Blanchard and Walberg.40 

The number of DSMC parcels varied from 10–500 × 106 for the 140 and 80 km cases, respectively, required to satisfy the minimum number of particles per cell ( N i > 5). The DSMC simulations were run on the HSUper cluster of Helmut Schmidt University on 256–2304 cores, depending on the number of cells and parcels.

Below 70 km altitude, the flow around the Viking 1 capsule became strictly continuous ( Kn < 0.01). For this regime, we employed the hy2Foam solver developed by Casseau et al.38 which in its core relies on the rhoCentralFoam solver from Ref. 47 and rhoReactingFoam solver within the OpenFOAM framework.rhoCentralFoam is a density-based shock capturing solver that discretizes convective fluxes by means of a total variation diminishing (TVD) scheme. The governing equations are the compressible Navier–Stokes equations that describe the conservation of mass, momentum and energy, expressed as
ρ t + · ( ρ u ) = 0 ,
(3)
( ρ u ) t + · ( ρ u u ) = p + · τ ¯ ¯ ,
(4)
E t + · ( ρ u E ) + · ( ρ u ) = · ( u T ) + · J ,
(5)
respectively. In Eqs. (3)–(5), ρ denotes the density defined as ρ = p / ( R s T ) according to the ideal-gas law, u is the velocity vector, t is the time, p is the pressure, τ ¯ ¯ is the viscous stress tensor, E is the total energy, T is the temperature, and J is the heat flux vector. The vector of conserved quantities is given by
U = ( ρ ρ u ρ E ) .
(6)
In order to capture thermal non-equilibrium effects, the hy2Foam solver includes the addition of electro-vibrational energy modes implemented using Park's two-temperature model48 and energy transfers between the vibrational and translational modes based on the Landau–Teller theory49 and work of Ref. 6. This yields additional dimensions in the vector of conserved quantities considering multiple species with subscript s as38 
U = ( ρ ρ s ρ u ρ E v m ρ E ) F inv = ( ρ u i ρ s u i ρ u i u j + p δ i j E v m u i ( E + p ) u i ) .
(7)
In Eq. (7), E v m represents the vibro-electronic energy of the molecule m, F is the (inviscid) flux vector, and δij is the Kronecker delta. From here, the non-equilibrium Navier–Stokes–Fourier (NSF) equation can be formulated for a fluid mixture consisting of s species and m molecules through
U t + F x = W ̇ ,
where W ̇ is the source term vector describing the production rate of species through chemical reactions. A detailed derivation of the presented relations and the implementation in hy2Foam can be found in the notes of Ref. 50 and the publications from Refs. 37 and 38.

Modeling reactions is achieved using finite-rate chemistry with Arrhenius rate coefficients based on Park's recommendation for Mars reentry two-temperature CFD problems.7,51 The complete list of reactions for the CFD cases can be found in Table IV. Species transport properties including viscosity, thermal conductivity, mixing rules, and diffusion were computed according to the works of Refs. 52–55, respectively.

The meshing approach for the CFD cases is outlined in Fig. 4. Similar to the DSMC approach, a baseline mesh was created including boundary layers with N BL = 20, k = 1.15, and y + = 0.1 based on the freestream properties. Then, due to the presence of discontinuities in the continuum regime, the mesh was iteratively refined at the bow shock. The maximum refinement level of the octree was set to 3, yielding a maximum of 20 × 106 cells for the 30 km case in order to satisfy our requirements. Since a fully unstructured meshing approach was employed, choosing appropriate numerical methods was important to obtain robustness, sufficient accuracy, and low numerical diffusion. We used the Crank–Nicolson technique56 to advance the solution in time using a coefficient of 0.9 and the MUSCL scheme proposed by57 as the divergence scheme coupled with the vanAlbada limiter.58 The CFD simulations were performed with an adaptive time-stepping approach based on a fixed local CFL number of 0.1 to ensure numerical robustness. This resulted in a minimum time step of 5 × 10−9 s for the 30 km case and required us to run the hy2Foam simulations also on HSUper on up to 4608 cores.

FIG. 4.

(a) Baseline mesh used for CFD simulations and (b) example of adaptively refined mesh for 50 km case.

FIG. 4.

(a) Baseline mesh used for CFD simulations and (b) example of adaptively refined mesh for 50 km case.

Close modal

The kω SST turbulence model proposed by Menter59 was employed to capture effects of residual viscosity and turbulent flow detachment that might occur around the heat shield and wake of the capsule at lower altitudes. Freestream flow properties are enforced at the farfield boundary (see Fig. 2) using a Dirichlet type fixedValue condition, while at the outlet Neumann type, zeroGradient conditions are used for all flow quantities. The maxwellSlipU wall boundary condition60 was utilized with an accomodation coefficient of α MW = 1.0, while the smoluchowskiJumpT boundary condition61 was employed for modeling temperature jumps occurring in regions of rarefaction close to the wall.

Due to signal strength (i.e., forces on the capsule) and measurement noise, only drag coefficient flight data are available for the rarefied and transitional flow regime and also subject to noise above 120 km. A comparison of predictions of cd from DSMC simulations with flight data is outlined in Fig. 5. The magnitude of mean relative error throughout all altitudes remains small at an average of δ ¯ = 1.8 %. With the exception of the 110 km case, the drag coefficient is generally slightly underpredicted.

FIG. 5.

Drag coefficient predicted by dsmcFoam+ solver compared with flight data.40 The gray area denotes the region where sensor noise introduced noticeable errors on the acceleration data ( h 120 km).

FIG. 5.

Drag coefficient predicted by dsmcFoam+ solver compared with flight data.40 The gray area denotes the region where sensor noise introduced noticeable errors on the acceleration data ( h 120 km).

Close modal

In Fig. 6, contours of different flow quantities are illustrated in the xy-plane cross section for the 80, 100, and 140 km cases. The first column indicates contours at 140 km altitude. It is clearly visible that no bow shock is present at 140 km altitude which can be attributed to very-low-density, rarefied flow around the capsule and partially specular surface–gas interactions. In the area in front of the capsule, multiple regions of very hot gas are formed ( T t r 10 000 K), in which some of the specularly reflected gas particles collide with the incoming freestream. In this flow regime, however, most of the incoming particles in the freestream impact the capsule's heat shield almost undisturbed due to their long mean free path. As a result, the drag coefficient reaches its maximum value at the 140 km case. The wake of the capsule is characterized by an area of extremely low gas density, indicated through a very long mean free path of the particles.

FIG. 6.

Contours of Mach number, translational temperature, pressure, and mean-free path for three different altitudes in the transitional (80 and 100 km) and rarefied flow regime (140 km). Contours generated in ParaView 5.12 using jet colormap with 256 levels, min. and max. values of each plot indicated on colormaps on the right side of row.

FIG. 6.

Contours of Mach number, translational temperature, pressure, and mean-free path for three different altitudes in the transitional (80 and 100 km) and rarefied flow regime (140 km). Contours generated in ParaView 5.12 using jet colormap with 256 levels, min. and max. values of each plot indicated on colormaps on the right side of row.

Close modal

As the capsule loses altitude, the flow becomes transitional and a (thickened) bow shock starts to appear which can be observed at the Mach number contours of the 100 km case. Consequently, the hot temperature region forms directly behind the bow shock and in front of the capsule. Effects of dissociation and increased density reduce the maximum translational temperature to max ( T t r ) = 15 000 K. Looking at the 80 km case, it is evident that the fundamental flow physics do not change as drastical compared to the 100 km case since the flow is still transitional. The most significant changes include the thinning of the bow shock, as the surrounding flow becomes increasingly continuous. In addition, a clearly recognizable recompression shock occurs behind the capsule, which is noticeable by a hot temperature zone in the wake and a shorter and more homogenously distributed mean free path in this region.

When the capsule traversed 70 km altitude, the flow became strictly continuous since Kn < 0.01. This section is concerned with predictions of the hy2Foam solver for altitudes between 70 and 30 km. For this range, flight data with small measurement errors is available as well as numerical data from a previous publication by Edquist20 where the LAURA code was used (48–31 km). Note that after a threshold of a x > 0.05 · g E, the RCS thrusters, which held the capsule at a constant α = 11.1 °, were deactivated to allow the vehicle to freely adjust the angle of attack based on its natural trim angle. The altitude roughly coincided with the beginning of the continuum regime at 70 km when the capsule underwent pitching motion around its center of gravity.

The angle of attack data is shown in Fig. 7(a), and we took the raw data for each altitude since we compare our results to raw acceleration data from which the aerodynamic coefficients were derived. In Figs. 7(b) and 7(c), results of the drag and lift coefficients are compared with flight data and previous numerical results, respectively. In general, predictions of the hy2Foam setup yield slightly lower values of the drag coefficient with a mean relative error of δ ¯ ( c d ) = 1.4 %, while the predicted lift coefficient does not indicate a clear tendency and exhibits a similar mean relative error of δ ¯ ( c l ) = 2.1 %. It is visible that the hy2Foam simulations show an improvement in both accuracy and trend compared to the LAURA results from Ref. 40.

FIG. 7.

(a) Flight data angle of attack. (b) and (c)Drag and lift coefficients predicted by hy2Foam solver compared with flight data40 and LAURA simulations by Edquist et al.62 

FIG. 7.

(a) Flight data angle of attack. (b) and (c)Drag and lift coefficients predicted by hy2Foam solver compared with flight data40 and LAURA simulations by Edquist et al.62 

Close modal

Contours of essential flow quantities for the 70, 50, and 30 km cases are illustrated in Fig. 8. At 70 km altitude, it is apparent that the bow shock has become much thinner and forms a discontinuity around the capsule. The translational temperature increases to almost 9500 K as a step function when the freestream fluid traverses the bow shock and subsequently undergoes thermochemical non-equilibrium. A thin region of hot compressed gas appears behind the bow shock and initiates dissociation of CO2 molecules into CO and O, indicated by the plot in the last row. When the capsule flew below 50 km, the process of maximum deceleration occurred and a large amount of CO2 dissociates due to the combination of high translational-vibrational energy and pressure in front of the capsule. The recompression shock becomes more apparent in the wake and its distance to the capsule reduces. At 30 km altitude, the bow shock moves slightly farther away from the heat shield and becomes wider due to the lower Mach number. The maximum temperature of T t r T v = 3000 K results in a drastic reduction of CO2 dissociation, hence the mass fraction of CO decreases to a maximum of 2%.

FIG. 8.

Contours of Mach number, translational temperature, pressure, and CO mass-fraction for three different altitudes in the continuum regime. Contours generated in ParaView 5.12 using jet colormap with 256 levels, min. and max. values of each plot indicated on colormaps on the right side of row.

FIG. 8.

Contours of Mach number, translational temperature, pressure, and CO mass-fraction for three different altitudes in the continuum regime. Contours generated in ParaView 5.12 using jet colormap with 256 levels, min. and max. values of each plot indicated on colormaps on the right side of row.

Close modal
Finally, a comparison of the maximum pressure coefficient cp is given as
c p = 2 p s ρ u 2 where p s = ρ R T M
(8)
is shown in Fig. 9. Expansion tube data from Ref. 63 suggested that between 40 and 60 km altitude, non-equilibrium effects impose a variation of the pressure coefficient. Our simulations of the 40 km case underline this suggestion since we obtain similar values of cp as well as noticable difference between the translational and vibrational temperature in the stagnation region (see 50 km case in Fig. 8). The relative error of predicted cp decreases at higher altitudes, suggesting the presence of strong (local) turbulent effects occurring at lower altitudes in the stagnation region between the bow shock and the front surface.
FIG. 9.

Comparison of stagnation point pressure coefficient with expansion tube data from Ref. 63 at different altitudes.

FIG. 9.

Comparison of stagnation point pressure coefficient with expansion tube data from Ref. 63 at different altitudes.

Close modal

Detailed views of the pressure coefficient distribution on the front surface of the capsule are presented in Fig. 10. The assumption that turbulent effects play an important role in the stagnation zone at lower altitudes is supported by an asymmetry of the pressure distribution visible at 30 km. Furthermore, the maximum obtained mass fraction CO at this altitude is 1.2%, while the vibrational temperature distribution is very similar to the translational temperature ( T t r max T v max = 2900 K), implying that the thermo-chemical flow characteristics are close to equilibrium and cannot attribute to strong local pressure or temperature fluctuations in the stagnation region.

FIG. 10.

Comparison of instantaneous pressure coefficient distribution on the front surface (axis dimensions in m), stagnation point highlighted by green square.

FIG. 10.

Comparison of instantaneous pressure coefficient distribution on the front surface (axis dimensions in m), stagnation point highlighted by green square.

Close modal

The hyStrath framework was employed to simulate both the hypersonic rarefied and continuum flow regime around the Viking 1 capsule in order to evaluate the framework for predicting the lift and drag coefficients during hypersonic reentry. Using a fully unstructured octree adaptive meshing approach, multi-temperature model, and specialized wall boundary conditions, very good agreement with flight data was achieved. In the high-altitude rarefied and transitional regime, where the dsmcFoam+ solver was employed, the formation of highly non-equilibrium flow and eventual formation of a (thickened) bow shock ( h 100 km) that drastically influence the flow physics and, therefore, the aerodynamic coefficients around the capsule, could be captured. Only the drag coefficient could be evaluated with flight data in the rarefied and transitional regime and showed excellent agreement with flight data analyzed by Blanchard and Walberg.40 The hy2Foam solver used for simulating the hypersonic continuum regime during Viking 1 reentry enjoyed improved validation due to the better availability of data for the continuum hypersonic coefficients. Here, the predictions for both the lift and drag coefficients could be, on the one hand, improved compared to a previous publication by Edquist20 and, on the other hand, the maximum altitude for employing a continuum solver increased from 48 to 70 km. Additionally, stagnation pressure coefficients obtained from a combination of flight data pressure gauges and expansion tube ground tests could be reproduced, indicating that the non-equilibrium flow in the stagnation region could be accurately modeled using the hy2Foam solver. Note that the SST model cannot fully capture the highly transient and turbulent effects in the stagnation zone, which become increasingly significant at lower altitudes and freestream Reynolds numbers. This is a finding that provides a basis for future studies of turbulent effects in the stagnation zone, especially using large eddy simulation (LES).

The present study serves as a basis for further research on reentry vehicles not only bound to Mars reentry. Next steps include evaluation of the framework and meshing approach for predicting aerothermal loads on Mars, Earth, and Titan reentry vehicles as well as simulating currently developed vehicles and future designs.

Computational resources (HPC-cluster HSUper) have been provided by the project hpc.bw, funded by dtec.bw—Digitalization and Technology Research Center of the Bundeswehr. dtec.bw is funded by the European Union—NextGenerationEU.

Support in the improvement of the dsmcFoam+ code and general assistance in the implementation of thermodynamic properties and reactions by Dr. Vincent Casseau is gratefully acknowledged.

The authors have no conflicts to disclose.

Maximilian Maigler: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Investigation (equal); Methodology (equal); Project administration (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Valentina Pessina: Conceptualization (supporting); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (supporting); Writing – review & editing (equal). Jochen Schein: Funding acquisition (lead); Project administration (equal); Supervision (lead).

Raw data were generated at the HSUper large scale facility. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

Atmospheric freestream properties and telemetry of investigated flow cases (Table I).

TABLE I.

Atmospheric freestream properties and telemetry of investigated flow cases.

h, km Flow regimea Kn α, ° u , ms-1 p , Pa T , K CO2, % N2, % Ar, %
30  3 × 10−5  12.1  2700  3.8 × 10  171  96 
40  9 × 10−5  10.8  4013  1.1 × 10  152  96 
50  3 × 10−4  10.5  4450  2.6 × 10  143  96 
60  2 × 10−3  11.0  4550  6.6 × 10−1  166  96 
70  5 × 10−3  12.2  4561  1.6 × 10−1  145  96 
80  3 × 10−2  11.1  4560  3.4 × 10−2  152  96 
90  5 × 10−2  11.1  4554  2.2 × 10−2  168  96 
100  1 × 10−1  11.1  4547  8.1 × 10−3  161  96 
110  6 × 10−1  11.1  4540  2.1 × 10−3  189  96 
120  1 × 100  11.1  4532  1.3 × 10−3  198  96 
130  8 × 100  11.1  4523  9.1 × 10−5  122  96 
140  4 × 101  11.1  4514  2.3 × 10−5  167  94 
h, km Flow regimea Kn α, ° u , ms-1 p , Pa T , K CO2, % N2, % Ar, %
30  3 × 10−5  12.1  2700  3.8 × 10  171  96 
40  9 × 10−5  10.8  4013  1.1 × 10  152  96 
50  3 × 10−4  10.5  4450  2.6 × 10  143  96 
60  2 × 10−3  11.0  4550  6.6 × 10−1  166  96 
70  5 × 10−3  12.2  4561  1.6 × 10−1  145  96 
80  3 × 10−2  11.1  4560  3.4 × 10−2  152  96 
90  5 × 10−2  11.1  4554  2.2 × 10−2  168  96 
100  1 × 10−1  11.1  4547  8.1 × 10−3  161  96 
110  6 × 10−1  11.1  4540  2.1 × 10−3  189  96 
120  1 × 100  11.1  4532  1.3 × 10−3  198  96 
130  8 × 100  11.1  4523  9.1 × 10−5  122  96 
140  4 × 101  11.1  4514  2.3 × 10−5  167  94 
a

C: Continuum (hy2Foam), T: Transitional (dsmcFoam+), R: Rarefied (dsmcFoam+).

Q-K reactions used in dsmcFoam+ simulations (Table II). Collision-averaged VSS model parameters based on Ref. 46 with T ref = 273 K (Table III).

TABLE II.

Q-K reactions used in dsmcFoam+ simulations.

Reaction type Nr. Reaction
Dissociation  CO2 + CO2 → CO + O + CO2 
  CO2 + CO → CO + O + CO 
  CO2 + N2 → CO + O + N2 
  CO2 + Ar → CO + O + Ar 
  N2 + N2 → N + N + N2 
  N2 + CO2 → N + N + CO2 
  N2 + O2 → N + N + O2 
  N2 + Ar → N + N + Ar 
  O2 + O2 → O + O + O2 
  10  O2 + CO2 → O + O + CO2 
  11  O2 + N2 → O + O + N2 
  12  O2 + Ar → O + O + Ar 
Exchange  13  NO + O N + O2 
  14  N2 + O NO + N 
  15  CO + O C + O2 
  16  CO2 + O CO + O2 
Reaction type Nr. Reaction
Dissociation  CO2 + CO2 → CO + O + CO2 
  CO2 + CO → CO + O + CO 
  CO2 + N2 → CO + O + N2 
  CO2 + Ar → CO + O + Ar 
  N2 + N2 → N + N + N2 
  N2 + CO2 → N + N + CO2 
  N2 + O2 → N + N + O2 
  N2 + Ar → N + N + Ar 
  O2 + O2 → O + O + O2 
  10  O2 + CO2 → O + O + CO2 
  11  O2 + N2 → O + O + N2 
  12  O2 + Ar → O + O + Ar 
Exchange  13  NO + O N + O2 
  14  N2 + O NO + N 
  15  CO + O C + O2 
  16  CO2 + O CO + O2 
TABLE III.

Collision-averaged VSS model parameters based on Ref. 46 with T ref = 273 K.

Species d ref (Å) ω (K) α VSS
CO2  4.647  0.693  1.373 
CO  4.101  0.726  1.341 
NO  3.983  0.716  1.425 
N2  3.911  0.693  1.351 
O2  3.773  0.702  1.391 
3.340  0.772  1.471 
3.402  0.753  1.477 
4.042  0.811  1.490 
Ar  3.832  0.700  1.384 
Species d ref (Å) ω (K) α VSS
CO2  4.647  0.693  1.373 
CO  4.101  0.726  1.341 
NO  3.983  0.716  1.425 
N2  3.911  0.693  1.351 
O2  3.773  0.702  1.391 
3.340  0.772  1.471 
3.402  0.753  1.477 
4.042  0.811  1.490 
Ar  3.832  0.700  1.384 

Collision-averaged VSS model parameters based on Ref. 46 with T ref = 273 K (Table IV).

TABLE IV.

Laminar finite-rate Arrhenius reaction rates used in hy2Foam simulations.

Reaction type Nr. Reaction A, m3 mol−1 s−1 β Ta, K
Dissociation  CO2 + M → CO + O + M  6.9 × 1017  −1.50  63 275 
  CO + M → C + O + M  2.3 × 1016  −1.00  129 000 
  O2 + M → O + O + M  2.0 × 1018  −1.50  59 750 
  N2 + M → N + N + M  7.0 × 1018  −1.60  113 200 
  NO + M → N + O + M  8.4 × 109  0.00  19 450 
Exchange  CO + O C + O2  3.9 × 1010  −0.18  69 200 
  CO2 + O CO + O2  2.1 × 1010  0.00  27 800 
  N2 + O NO + N  6.4 × 1014  −1.00  38 370 
  NO + O N + O2  8.4 × 1019  0.00  19 450 
Reaction type Nr. Reaction A, m3 mol−1 s−1 β Ta, K
Dissociation  CO2 + M → CO + O + M  6.9 × 1017  −1.50  63 275 
  CO + M → C + O + M  2.3 × 1016  −1.00  129 000 
  O2 + M → O + O + M  2.0 × 1018  −1.50  59 750 
  N2 + M → N + N + M  7.0 × 1018  −1.60  113 200 
  NO + M → N + O + M  8.4 × 109  0.00  19 450 
Exchange  CO + O C + O2  3.9 × 1010  −0.18  69 200 
  CO2 + O CO + O2  2.1 × 1010  0.00  27 800 
  N2 + O NO + N  6.4 × 1014  −1.00  38 370 
  NO + O N + O2  8.4 × 1019  0.00  19 450 
1.
G. A.
Bird
, “
Approach to translational equilibrium in a rigid sphere gas
,”
Phys. Fluids
6
,
1518
(
1963
).
2.
G. A.
Bird
,
The DSMC Method: Version 1.2
(
CreateSpace
,
2013
).
3.
H.
Deng
,
T.
Ozawa
, and
D. A.
Levin
, “
Analysis of chemistry models for DSMC simulations of the atmosphere of Io
,”
J. Thermophys. Heat Transfer
26
,
36
(
2012
).
4.
E.
Chabut
, “
Internal energy exchange and dissociation probability in dsmc molecular collision models
,”
AIP Conf. Proc.
1084
,
389
(
2008
).
5.
G.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
(
Clarendon Press
,
Oxford
,
1994
).
6.
R. C.
Millikan
and
D. R.
White
, “
Systematics of vibrational relaxation
,”
J. Chem. Phys.
39
,
3209
(
1963
).
7.
C.
Park
,
J. T.
Howe
,
R. L.
Jaffe
, and
G. V.
Candler
, “
Review of chemical-kinetic problems of future NASA missions. II—Mars entries
,”
J. Thermophys. Heat Transfer
8
,
9
(
1994
).
8.
J. F.
Padilla
and
I. D.
Boyd
, “
Assessment of gas-surface interaction models for computation of rarefied hypersonic flow
,”
J. Thermophys. Heat Transfer
23
,
96
(
2009
).
9.
S. J.
Plimpton
,
S. G.
Moore
,
A.
Borner
,
A. K.
Stagg
,
T. P.
Koehler
,
J. R.
Torczynski
, and
M. A.
Gallis
, “
Direct simulation Monte Carlo on petaflop supercomputers and beyond
,”
Phys. Fluids
31
,
086101
(
2019
).
10.
G.
LeBeau
, “
A parallel implementation of the direct simulation Monte Carlo method
,”
Comput. Methods Appl. Mech. Eng.
174
,
319
(
1999
).
11.
S.
Dietrich
and
I.
Boyd
, “
A scalar optimized parallel implementation of the DSMC method
,” in
32nd Aerospace Sciences Meeting and Exhibit
,
1994
.
12.
G. A.
Bird
, “
The DS2V/3V Program Suite for DSMC Calculations
,”
AIP Conf. Proc.
762
,
541
546
(
2005
).
13.
M. S.
Ivanov
,
A. V.
Kashkovsky
,
S. F.
Gimelshein
,
G.
Markelov
,
A. A.
Alexeenko
,
Y. A.
Bondar
,
G.
Zhukova
,
S.
Nikiforov
, and
P.
Vaschenkov
, “
Smile system for 2D/3D DSMC computations
,” in
Proceedings of 25th International Symposium on Rarefied Gas Dynamics
,
St. Petersburg, Russia
, 21 July 2006 (
2006
), pp.
21
28
.
14.
T.
Scanlon
,
E.
Roohi
,
C.
White
,
M.
Darbandi
, and
J.
Reese
, “
An open source, parallel dsmc code for rarefied gas flows in arbitrary geometries
,”
Comput. Fluids
39
,
2078
(
2010
).
15.
A.
Borner
,
J. B. E.
Meurisse
, and
N. N.
Mansour
, “
DSMC simulations of Mars Science Laboratory entry from rarefied to continuum conditions
,” International Symposium on Rarefied Gas Dynamics, ARC-E-DAA-TN59250 (
2018
).
16.
J.
Moss
,
R.
Wilmoth
, and
J.
Price
, “
DSMC simulations of blunt body flows for Mars entries: Mars pathfinder and Mars microprobe capsules
,” J. Spacecr. Rockets
36
(
3
),
330
339
(
1999
).
17.
J. N.
Moss
,
R. C.
Blanchard
,
R. G.
Wilmoth
, and
R. D.
Braun
, “
Mars Pathfinder rarefied aerodynamics: Computations and measurements
,”
J. Spacecr. Rockets
36
,
330
(
1999
).
18.
R. C.
Blanchard
,
R. G.
Wilmoth
, and
J. N.
Moss
, “
Aerodynamic flight measurements and rarefied-flow simulations of Mars entry vehicles
,”
J. Spacecr. Rockets
34
,
687
(
1997
).
19.
K. B.
Thompson
,
C. O.
Johnston
,
B. R.
Hollis
, and
V.
Lessard
, “
Recent Improvements to the LAURA and HARA Codes
,” in
AIAA Aviation 2020 Forum
,
2020
.
20.
K.
Edquist
, “
Computations of Viking lander capsule hypersonic aerodynamics with comparisons to ground and flight data
,” AIAA Paper No. 2006-6137,
2006
.
21.
K.
Edquist
,
A.
Dyakonov
,
M.
Wright
, and
C.
Tang
, “
Aerothermodynamic design of the Mars science laboratory backshell and parachute cone
,” AIAA Paper No. 2009-4078,
2009
.
22.
D.
Schwamborn
,
T.
Gerhold
, and
R.
Heinrich
, “
The DLR TAU-code: Recent applications in research and industry
,” in
ECCOMAS CFD 2006: Proceedings of the European Conference on Computational Fluid Dynamics, Egmond aan Zee, The Netherlands, September 5-8, 2006
(
Delft University of Technology; European Community on Computational Methods
, (
2006
).
23.
D.
Potter
,
S.
Karl
,
M.
Lambert
, and
K.
Hannemann
, “
Computation of radiative and convective contributions to Viking afterbody heating
,” in
44th AIAA Thermophysics Conference
(
AIAA
, 2014) (
2014
).
24.
I. V.
Egorov
and
M.
Pugach
, “
Numerical simulation of flow over space vehicle in Martian atmosphere,”
in
7th European Conference for Aeronautics and Space Sciences
(
2017
).
25.
M. J.
Wright
,
G. V.
Candler
, and
D.
Bose
, “
Data-parallel line relaxation method for the Navier-Stokes equations
,”
AIAA J.
36
,
1603
(
1998
).
26.
D. D.
Morabito
,
E.
Papajak
,
T.
Hedges
,
D.
Saunders
,
P.
Ilott
,
C.
Jin
,
P.
Fieseler
,
M.
Kobayashi
, and
M.
Shihabi
, “
The Mars 2020 entry, descent, and landing communications brownout and blackout at ultra-high frequency
,” IPN Progress Report No. 42-233 (
NASA
,
2023
).
27.
A.
Viviani
,
G.
Pezzella
, and
C.
Golia
, “
Aerothermodynamic analysis of a space vehicle for manned exploration missions to Mars
,” in
27th International Congress of the Aeronautical Sciences
,
2010
.
28.
A.
Viviani
and
G.
Pezzella
, “
Nonequilibrium computational flowfield analysis for the design of Mars manned entry vehicles
,”
Prog. Flight Phys.
5
,
493
(
2013
).
29.
M.
Wright
,
J.
Olejniczak
,
K.
Edquist
,
E.
Venkatapathy
, and
B.
Hollis
, “
Status of aerothermal modeling for current and future Mars exploration missions
,” in
IEEE Aerospace Conference
,
2006
.
30.
K.
Zhong
,
X.
Wang
, and
C.
Yan
, “
Wall temperature effects on hypersonic aerodynamics of the Mars entry capsule,” in
9th International Conference on Mechanical and Aerospace Engineering (ICMAE)
(
IEEE
,
2018
), pp.
141
145
.
31.
K.
Zhong
,
H.
Zhao
,
Z.
Xu
,
F.
Zhang
, and
W.
Zhao
, “
Wall catalytic effects on aerodynamics of the mars entry capsule
,” in
13th International Conference on Mechanical and Aerospace Engineering (ICMAE)
(IEEE,
2022
) pp.
443
447
.
32.
P.
Noeding
, “
Technologies for safe and controlled martian entry
,” Technical Report No. Del. No. D 6.1 (
SACOMAR
,
2011
).
33.
S.
Zeng
,
Z.
Yuan
,
W.
Zhao
, and
W.
Chen
, “
Numerical simulation of hypersonic thermochemical nonequilibrium flows using nonlinear coupled constitutive relations
,”
Chin. J. Aeronaut.
36
,
63
(
2023
).
34.
R. N.
Gupta
,
J. M.
Yos
,
R. A.
Thompson
, and
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
,” NASA Reference Publication 1232 (
1990
).
35.
W.
Jones
and
A.
Cross
, “
Electrostatic-probe measurements of plasma parameters for two reentry flight experiments at 25000 feet per second
,” NASA TN D-6617 (
1972
).
36.
S.
Zeng
,
J.
Yang
,
W.
Zhao
,
Y.
Huang
,
Z.
Jiang
, and
W.
Chen
, “
Computational study of lateral jet interaction in hypersonic thermochemical non-equilibrium flows using nonlinear coupled constitutive relations
,”
Phys. Fluids
35
,
116117
(
2023
).
37.
V.
Casseau
,
R. C.
Palharini
,
T. J.
Scanlon
, and
R. E.
Brown
, “
A two-temperature open-source CFD model for hypersonic reacting flows, part one: Zero-dimensional analysis
,”
Aerospace
3
,
34
(
2016
).
38.
V.
Casseau
,
D. E. R.
Espinoza
,
T. J.
Scanlon
, and
R. E.
Brown
, “
A two-temperature open-source CFD model for hypersonic reacting flows, part two: Multi-dimensional analysis
,”
Aerospace
3
,
45
(
2016
).
39.
H. G.
Weller
,
G.
Tabor
,
H.
Jasak
, and
C.
Fureby
, “
A tensorial approach to computational continuum mechanics using object-oriented techniques
,”
Comput. Phys.
12
,
620
(
1998
).
40.
R. C.
Blanchard
and
G. D.
Walberg
, “
Determination of the hypersonic-continuum/rarefied-flow drag coefficient of the Viking lander capsule 1 aeroshell from flight data
,” NASA Technical Paper No. 1793 (
Langley Research Center
,
1980
).
41.
A. O.
Nier
,
W. B.
Hanson
,
A.
Seiff
,
M. B.
McElroy
,
N. W.
Spencer
,
R. J.
Duckett
,
T. C. D.
Knight
, and
W. S.
Cook
, “
Composition and structure of the martian atmosphere: Preliminary results from Viking 1
,”
Science
193
,
786
(
1976
).
42.
C.
White
,
M.
Borg
,
T.
Scanlon
,
S.
Longshaw
,
B.
John
,
D.
Emerson
, and
J.
Reese
, “
dsmcFoam+: An OpenFOAM based direct simulation Monte Carlo solver
,”
Comput. Phys. Commun.
224
,
22
(
2018
).
43.
G. A.
Bird
, “
The Q-K model for gas-phase chemical reaction rates
,”
Phys. Fluids
23
,
106101
(
2011
).
44.
T. J.
Scanlon
,
C.
White
,
M. K.
Borg
,
R. C.
Palharini
,
E.
Farbar
,
I. D.
Boyd
,
J. M.
Reese
, and
R. E.
Brown
, “
Open-source direct simulation Monte Carlo chemistry modeling for hypersonic flows
,”
AIAA J.
53
,
1670
(
2015
).
45.
K.
Koura
and
H.
Matsumoto
, “
Variable soft sphere molecular model for inverse-power-law or Lennard-Jones potential
,”
Phys. Fluids A
3
,
2459
(
1991
).
46.
M.
Pfeiffer
, “
An optimized collision-averaged variable soft sphere parameter set for air, carbon, and corresponding ionized species
,”
Phys. Fluids
34
,
117110
(
2022
).
47.
A.
Kurganov
and
E.
Tadmor
, “
New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations
,”
J. Comput. Phys.
160
,
241
(
2000
).
48.
C.
Park
, “
Assessment of two-temperature kinetic model for ionizing air
,”
J. Thermophys. Heat Transfer
3
,
233
(
1989
).
49.
L. D.
Landau
and
E.
Teller
, “
On the theory of sound dispersion
,”
Phys. Z. Sowjetunion
10
,
34
43
(
1936
).
50.
G. V.
Candler
and
I.
Nompelis
, “
Computational fluid dynamics for atmospheric entry
,” in Non-Equilibrium Dynamics: From Physical Models to Hypersonic Flights (The von Karman Institute for Fluid Dynamics, Belgium, 2009).
51.
C.
Park
,
Nonequilibrium Hypersonic Aerothermodynamics
(
Wiley
,
1990
).
52.
F. G.
Blottner
,
M.
Johnson
, and
M.
Ellis
, “
Chemically reacting viscous flow program for multi-component gas mixtures
,” Sandia National Lab. Report No. SC-RR-70-754 (
Sandia National Laboratories
,
Albuquerque
,
NM
,
1971
).
53.
A.
Eucken
, “
On the temperature dependence of the thermal conductivity of several gases
,”
Phys. Z
12
,
1101
1107
(
1911
).
54.
B.
Armaly
and
K.
Sutton
, “
Viscosity of multicomponent partially ionized gas mixtures
,” in
15th Thermophysics Conference
,
1980
.
55.
H. J. V.
Tyrrell
, “
The origin and present status of Fick's diffusion law
,”
J. Chem. Educ.
41
,
397
(
1964
).
56.
J.
Crank
and
P.
Nicolson
, “
A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type
,”
Math. Proc. Cambridge Philos. Soc.
43
,
50
(
1947
).
57.
B.
van Leer
, “
Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method
,”
J. Comput. Phys.
32
,
101
(
1979
).
58.
G. D.
van Albada
,
B.
van Leer
, and
W. W.
Roberts
, “
A comparative study of computational methods in cosmic gas dynamics
,” in
Upwind and High-Resolution Schemes
(
Spri
nge
r
Berlin Heidelberg,
Berlin, Heidelberg
,
1997
), pp.
95
103
.
59.
F. R.
Menter
, “
Two-equation Eddy-viscosity turbulence models for engineering applications
,”
AIAA J.
32
,
1598
(
1994
).
60.
J. C.
Maxwell
, “
VII. On stresses in rarified gases arising from inequalities of temperature
,”
Philos. Trans. R. Soc. London
170
,
231
256
(
1879
)
61.
M. V.
Smoluchowski
, “
Über Wärmeleitung in verdünnten Gasen
,”
Ann. Phys. Chem.
64
,
101
(
1898
).
62.
K.
Edquist
,
M.
Wright
, and
G.
Allen
, “
Viking afterbody heating computations and comparisons to flight data
,” AIAA Paper No. 2006-386,
2006
.
63.
C.
Miller
and
L. R.
Center
, “
Shock shapes on blunt bodies in hypersonic-hypervelocity helium, air, and CO2 flows, and calibration results in Langley 6-inch expansion tube
,” Technical Report No. NASA-TN-D-7800 (Washington, D.C.,
1975
).