Next Article in Journal
Correction: Iemma, U. Theoretical and Numerical Modeling of Acoustic Metamaterials for Aeroacoustic Applications. Aerospace 2016, 3, 15
Next Article in Special Issue
Modeling Aerodynamics, Including Dynamic Stall, for Comprehensive Analysis of Helicopter Rotors
Previous Article in Journal
Gaskinetic Modeling on Dilute Gaseous Plume Impingement Flows
Previous Article in Special Issue
Wing Tip Drag Reduction at Nominal Take-Off Mach Number: An Approach to Local Active Flow Control with a Highly Robust Actuator System
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Two-Temperature Open-Source CFD Model for Hypersonic Reacting Flows, Part Two: Multi-Dimensional Analysis †

by
Vincent Casseau
1,2,*,‡,
Daniel E. R. Espinoza
1,2,‡,
Thomas J. Scanlon
1,2,‡ and
Richard E. Brown
2,‡
1
James Weir Fluids Laboratory, University of Strathclyde, Glasgow G1 1XJ, UK
2
Centre for Future Air-Space Transportation Technology, University of Strathclyde, Glasgow G1 1XJ, UK
*
Author to whom correspondence should be addressed.
This paper is an extended version of our paper published in Casseau, V.; Scanlon, T.J.; John, B.; Emerson, D.R.; Brown, R.E. Hypersonic simulations using open-source CFD and DSMC solvers. In DSMC and Related Simulations, Proceedings of the 30th International Symposium on Rarefied Gas Dynamics, Victoria, BC, Canada, 10–15 July 2016; Ketsdever, A.; Struchtrup, H., Eds.; AIP Publishing, Melville, NY, USA; Volume 1786, 050006.[...]
These authors contributed equally to this work.
Aerospace 2016, 3(4), 45; https://doi.org/10.3390/aerospace3040045
Submission received: 30 September 2016 / Revised: 25 November 2016 / Accepted: 8 December 2016 / Published: 14 December 2016
(This article belongs to the Special Issue State-of-the-Art Aerospace Sciences and Technologies in Europe)

Abstract

:
hy2Foam is a newly-coded open-source two-temperature computational fluid dynamics (CFD) solver that has previously been validated for zero-dimensional test cases. It aims at (1) giving open-source access to a state-of-the-art hypersonic CFD solver to students and researchers; and (2) providing a foundation for a future hybrid CFD-DSMC (direct simulation Monte Carlo) code within the OpenFOAM framework. This paper focuses on the multi-dimensional verification of hy2Foam and firstly describes the different models implemented. In conjunction with employing the coupled vibration-dissociation-vibration (CVDV) chemistry–vibration model, novel use is made of the quantum-kinetic (QK) rates in a CFD solver. hy2Foam has been shown to produce results in good agreement with previously published data for a Mach 11 nitrogen flow over a blunted cone and with the dsmcFoam code for a Mach 20 cylinder flow for a binary reacting mixture. This latter case scenario provides a useful basis for other codes to compare against.

Graphical Abstract

1. Introduction

The Knudsen number, defined as the ratio of the mean free path of the gas particles to the characteristic length of the problem, is commonly employed to gauge the degree of rarefaction of a gas. During the entry of a planetary atmosphere at hypervelocities, a craft will traverse the full range of the Knudsen number, from the free-molecular regime down to the continuum regime. Practically, this translates into the need for different numerical techniques to resolve the flow-field past such hypersonic bodies.
The direct simulation Monte Carlo (DSMC) method developed by Bird [1] is a particle-based methodology that is particularly well-suited for computing high Knudsen number flows, typically above 0.05, while the conventional computational fluid dynamics (CFD) approach that solves the Navier–Stokes–Fourier (NSF) equations is generally adopted for the lower range, below 0.005. In between lies an intermediate zone where DSMC is computationally prohibitive and where conventional CFD fails due to the presence of non-continuum regions within the flow-field. To extend the range of applicability of the NSF equations towards the continuum–transition regime, Park formulated the two-temperature CFD model [2], thus distinguishing the trans-rotational energy pool from the vibro-electronic energy pools and modelling the energy exchange processes via vibrational rate equations.
Open-source solvers dedicated to the study of the hypersonic regime are very scarce. Among such codes are hy2Foam and dsmcFoam. hy2Foam [3] is a new two-temperature CFD solver developed at the University of Strathclyde [4,5] within the framework of the open-source CFD platform OpenFOAM [6]. It aims at (1) giving open-source access to a state-of-the-art hypersonic CFD solver to students and researchers; and (2) providing a foundation for a future hybrid CFD-DSMC code [7].
For rarefied gas flow environments, OpenFOAM possesses a DSMC solver called dsmcFoam that has been developed and validated at the same University [8,9]. The two main features of dsmcFoam are that vibrational-translational post-collision energy redistribution is executed using a quantum version of the Larsen–Borgnakke procedure [8,10] and that chemical reactions are described by Quantum-Kinetic (QK) theory, as initially proposed by Bird [11].
The primary objective of this paper concerns the verification and validation of hy2Foam for multi-dimensional case scenarios. The secondary goal is the assessment, for a two-dimensional gas flow, of one of the conclusions formulated in [5], namely that the combination of a particular chemistry–vibration model with chemical rates derived from QK theory provides a satisfactory consistency between the CFD and DSMC codes.

2. Methodology

This section essentially describes the state-of-the-art numerical procedures that have been implemented in hy2Foam to solve hypervelocity flow-fields in near-thermal equilibrium.

2.1. Non-Equilibrium Navier–Stokes–Fourier Equations

The computation of transient hypersonic reacting flows in the continuum regime traditionally employs the non-equilibrium Navier–Stokes–Fourier (NSF) equations. These are shown below in flux–divergence form using a Cartesian coordinate system for a mixture composed of N s species, including N m molecules [12]:
U t + F i , i n v F i , v i s x i = W ˙ .
The vector of conserved quantities, U , is defined as
U = ρ , ρ s , ρ u , ρ v , ρ w , E v e , m , E T , s N s , m N m ,
where u, v, and w are the components of the velocity vector. ρ is the mass density of the fluid and ρ s is the partial density of species s. The flux vectors are split into inviscid and viscous contributions and are written as follows:
F i , i n v = ρ u i ρ s u i J ρ u i u + δ i 1 p ρ u i v + δ i 2 p ρ u i w + δ i 3 p E v e , m u i J ( E + p ) u i r e h r J r , i ,
F i , v i s = 0 J s , i τ i 1 τ i 2 τ i 3 q v e , i , m e v e , m J m , i τ i j u j q t r , i q v e , i r e h r J r , i ,
where the index i refers to one of the three dimensions of space and δ is the Kronecker delta. E and E v e , m represent, respectively, the total energy and the total vibro-electronic energy for molecule m. In the remainder of the article, the index t r denotes the trans-rotational energy mode, v e the vibro-electronic energy mode and e the electron energy mode. h s is the enthalpy per unit mass of species s, while the pressure, p, is recovered from the partial pressures using Dalton’s law
p = s e p s + p e = s e ρ s R s T t r + ρ e R e T v e , ref ,
where R s is the specific gas constant. Finally, the electron temperature is set to the vibro-electronic temperature of a reference particle, T v e , ref .
In Equation (4), τ i j represent the components of the viscous stress tensor, which can be written as follows:
τ i j = μ u i x j + u j x i + ( λ + μ b ) u k x k δ i j ,
where μ is the shear viscosity, λ is the second viscosity coefficient, and μ b is the bulk viscosity. This work assumes that Stokes’ hypothesis is valid, namely that μ b = 0 and that the shear and second viscosity are not independent quantities but given by the relation λ = 2 3 μ . In hy2Foam, the species shear viscosities can be considered as temperature-independent, or either set to follow Blottner’s formula [13] or a power law. The spatial components of the heat conduction vector are assumed to follow Fourier’s law
q t r , i = s q t r , i , s = s κ t r , s T t r x i for s N s \ { e }
and
q v e , i = s q v e , i , s = s κ v e , s T v e , s x i for s N s .
The thermal diffusivities κ t r , s and κ v e , s can be considered as temperature-independent, or set to follow Eucken’s formula [12,14]. Mixture quantities are recovered from species quantities using a mixing rule. The Wilke [15], Gupta [16], and Armaly and Sutton [17] mixing rules have been implemented, and the latter one is being used in this work, as recommended in the review of Palmer and Wright [18].
The diffusion fluxes, J s , i , are governed by Fick’s law with a correction term to ensure mass conservation [19]
J s , i = I s , i Y s r e I r , i , for s , r N s \ { e } , = M e r e C r × J r , i M r otherwise ,
with
I s , i = ρ D s Y s x i ,
where Y is the mass-fraction and with the effective diffusion coefficient of species s in the mixture defined as [20]
D s = 1 X s r s X r D s , r 1 for s , r N s \ { e } ,
in which X is the molar-fraction, C r represents the charge of species r and M is the molecular weight. The binary diffusion coefficients, D s , r , are modelled using binary collision integrals Δ s , r as [21]
D s , r = k B T t r p Δ s , r .
Finally, the source term vector, W ˙ , is written as
W ˙ = 0 , ω ˙ s , 0 , 0 , 0 , ω ˙ v , m , 0 T ,
where ω ˙ s is the net mass production of species s and ω ˙ v , m , for m N m , is defined as follows:
ω ˙ v , m = Q m , V T + Q m , V V + Q m , C V + Q m , e V , if   the   reference   molecule   for   e   is   not m , = Q m , V T + Q m , V V + Q m , C V + Q m , e V + Q h e + Q e i + U · p e otherwise ,
where U represents the velocity vector. The meaning of the different vibrational source terms that appear in Equation (14) is given in the subsequent paragraphs.

2.2. Energy Transfers

The decomposition of the total energy, E, as shown in Equation (15), allows the isolation of different energy pools, each of them being described by a specific temperature and exchanging energy with the other pools:
E = 1 2 ρ i u i 2 + s e E t , s + s e E r , s + s e E v , s + s e E e l , s + E e + s e ρ s h s ° .
In order of appearance in Equation (15) are the kinetic, translational (index t), rotational (r), vibrational (v), electronic ( e l ), electron (e), and chemical energies. u i are the spatial components of the velocity vector and h s ° is the standard enthalpy of formation of species s. The relationship between the total energy of a specific mode and the energy per unit mass of species s is given in Equation (16) for the electronic energy mode
E e l = s e E e l , s = s e ρ s e e l , s
and the different energies per unit mass of species are detailed below for each mode:
e t = 3 2 R s T t r ,
e r = R s T t r ,
e v , s = R s θ v , s e x p θ v , s T v e , s 1 ,
e e = R e T v e , ref ,
e e l , s = R s i = 1 g i , s θ e l , i , s e x p θ e l , i , s / T v e , s i = 0 g i , s e x p θ e l , i , s / T v e , s .
The simple harmonic oscillator model [11] is utilized for the vibrational mode, with θ v , s being the characteristic vibrational temperature for species s. θ e l , i , s and g i , s are the characteristic electronic temperature and the degeneracy degree of a given electronic state i for species s, respectively. The two-temperature model formulation adopted in hy2Foam then consists in rewriting Equation (15) as follows:
E = 1 2 ρ i u i 2 + E t r + s e E v e , s + E e + s e ρ s h s °
with
E t r = E t + E r
and
E v e , s = E v , s + E e l , s .
The vibrational source terms introduced into the NSF equations in the preceding paragraph are now described in sequence. The energy exchange between the trans-rotational and the vibro–electronic energy modes (V–T), designated by Q m , V T , is dictated by the Landau–Teller equation [22] and may be written as
Q m , V T = ρ m e v e , m ( T v e , m ) t = ρ m e v e , m ( T t r ) e v e , m ( T v e , m ) τ m , V T m N m ,
where τ m , V T is the molar-averaged V–T relaxation time. This former quantity is evaluated as the summation of the Millikan and White semi-empirical correlation [23] and Park’s correction factor [2], and later denoted as MWP. It becomes
τ m , V T = s e X s s e X s / τ m s , V T m N m ,
where X is the species molar-fraction and τ m s , V T is the interspecies relaxation time expressed as
τ m s , V T = 1 p e x p A m , s T t r 1 / 3 B m , s 18.42 + 1 c ¯ m σ v , m n m , s with   p   in   atm .
For a given colliding pair ( m , s ), the tabulated values of A m , s and B m , s used in this paper can be found in [4]. The collision-limited relaxation time introduced by Park [2] is a function of the average molecular speed, c ¯ m , the limited collision cross-section, σ v , m , and the number density of the colliding pair, n m , s .
The second V–T model tested in this work is the SSH theory named after Schwartz, Slawsky, and Herzfeld [24]. The coefficients of this latter model are taken from the work of Thivet [25] and a blended model is created with the MWP formulation for molecule-atom collisions.
Vibrational–vibrational (V–V) energy transfer is denoted by Q m , V V and modelled according to the formulation given by Knab et al. [26,27], as described in [4]. In Equation (14), Q h e , Q m , e V , and Q e i stand for the energy exchange between free-electrons and heavy-particles, free-electrons and the vibrational mode, and the vibrational energy removal due to electron impact ionisation, respectively. They are assumed to be zero in this work. The term in U · p e represents an approximation to the work done on the electrons by the electric field induced by the electron pressure gradient [28]. Finally, the vibro-electronic energy added or removed by chemical reactions to species m, and represented by Q m , C V , is modelled either with the Park TTv model [2] or the coupled vibration-dissociation-vibration (CVDV) model in hy2Foam [4].

2.3. Departure from the Continuum Regime

The local gradient-length Knudsen number introduced by Boyd [29] is often used as a breakdown parameter in the literature (e.g., in [30,31,32,33]). For a local macroscopic flow quantity ϕ, it is written as
K n G L L ϕ = λ ϕ | ϕ | ,
where λ is the local mean free path of the gas molecules and ϕ can either be the gas density, the temperature, or the magnitude of the velocity vector | U | (or m a x | U | , a in the denominator, where a is the speed of sound for low-speed regions [34]). The degree of local continuum breakdown is then evaluated as the maximum for all of the aforementioned flow quantities
K n G L L = m a x K n G L L ρ , K n G L L T t r , K n G L L T v e , K n G L L | U | ,
and the regions where the continuum assumption does not hold are identified when K n G L L exceeds a given breakdown value, K n B r , typically being taken equal to 0.05. This parameter is key in hybrid CFD-DSMC codes.

3. Results and Discussion

3.1. Mach 11.3 Blunted Cone

In this section, a non-reacting nitrogen flow past a blunted cone is examined at Mach 11.3. The case is composed of a 6.35 mm-radius nose followed by a flat plate forming a 25° angle with the free-stream flow direction and whose streamwise extension is 5 cm. The 2D axisymmetric mesh that has been employed is shown in Figure 1. The structured grid is aligned with the bow shock and consists of 600 by 200 cells. The first spacing at the wall surface is set to 2 μ m by default and to 10 μ m for comparison purposes.
The initial conditions of this case scenario are given in Table 1. The free-stream velocity corresponding to a Mach number of 11.3 is 2764.5 m/s. The free-stream temperature and pressure are 144.4 K and 21.9 Pa, respectively. The wall is assumed to be isothermally heated to a temperature of 297.2 K. This simulation uses the MWP formulation for V–T energy exchange and the Blottner and Eucken formulas to compute transport properties. The variable hard sphere model is chosen for the calculation of the mean-free-path and the Knudsen number is computed using the streamwise extent of the cone as the characteristic length. The non-equilibrium boundary conditions employed at the cone surface are the first-order Smoluchowski temperature jump [35] and Maxwell velocity slip [36], with the accommodation coefficients being taken as equal to unity. The set-up for this case can be found in the Supplementary Materials.
The hy2Foam solver uses the numerics present in the official OpenFOAM application rhoCentralFoam [37], namely the KNP central-upwind schemes of Kurganov, Noelle and Petrova [38]. Hence, hy2Foam is first order accurate in time and second order accurate in space. Improving the numerics of the solver to better capture the near-wall region notably is an ongoing area of research. Convergence was achieved after 2.8 h of computations on the ARCHIE-WeSt (Academic and Research Computer Hosting Industry Enterprise in the West of Scotland) High Performance Computer [39]. The run used 24 Intel Xeon X5650 2.66 GHz cores (Santa Clara, CA, USA) with 48 GB RAM and 4xQDR Infiniband Interconnect computer-networking communications.
This particular configuration has already been studied by Wang and Boyd [40], using the MONACO DSMC code and a Navier–Stokes CFD solver from the University of Michigan. The results from these simulations are reported in the subsequent graphs with the denomination DSMC: MONACO and CFD: Michigan. Moreover, experimental data is also shown in Figure 2d,f and correspond to the run 31 of the CUBRC (CUBRC Inc., Buffalo, NY, USA) experiments [41].
The emphasis is placed on stagnation line data and then on surface properties such as the pressure coefficient, C p , friction coefficient, C f , and Stanton number, S t , respectively, given by
C p = p p 0.5 ρ U 2 ,
C f = τ 0.5 ρ U 2 ,
S t = q 0.5 ρ U 3 ,
those essential aerothermodynamic quantities being shown in Figure 2d–f.
In Figure 2a–c, the pressure, temperature, and velocity stagnation line solutions given by the CFD code of Wang and Boyd are very similar to the ones produced by hy2Foam. It can be seen in Figure 2a that the vibrational mode is barely excited for this case scenario. The single-temperature model version of hy2Foam gave the same results. The small discrepancies concerning the shock stand-off distance can easily be explained by the difference in grid point density along the symmetry axis (indeed, the spatial extension of the domain normal to the body adopted here is about 40% larger than the one in [40]) and by the use of a different viscosity model (the power law was used in [40]). As it is the case in most simulations of hypersonic flow-fields, the bow shock thickness given by the NSF equations is clearly under-predicted, as shown by the different DSMC profiles.
Once again in Figure 2d–f, an excellent agreement is found between hy2Foam and the Michigan NSF solver as the profiles of the surface quantities are shown to be superimposed with a first spacing of 10 μ m. If the mesh is further refined to 2 μ m, a small decrease in the Stanton number is observed. In conclusion, it is thought that the mesh used in [40] had a first spacing close to 10 μ m.
The variable hard sphere mean-free-path computed by hy2Foam in the wall vicinity is of the order of 2 μ m. Good DSMC practice dictates that the mean-free-path to cell-size ratio should exceed one so one could argue that the DSMC mesh employed might not be fine enough near the wall to accurately capture the surface aerothermodynamic coefficients along the body, e.g., the peak amplitude of the skin-friction.
Finally, this case is a good illustration that a hypersonic simulation can now be carried out using open-source packages all the way from pre- to post-processing using Gmsh (version 2.11.0) [42] as a mesher, hy2Foam as a solver, Paraview [43] (version 4.1.0, Sandia Corporation, Albuquerque, NM, USA and Kitware Inc., Clifton Park, NY, USA) as a visualization utility, and Gnuplot [44] (version 5.0) as a grapher.

3.2. Mach 20 Cylinder

This section focuses on the hypersonic Mach 20 flow of nitrogen over a two-dimensional circular cylinder of radius R = 1 m. A symmetry plane exists about the y = 0 plane allowing the modelling solely of the upper half of the domain to be considered. The streamwise extent of the computational domain spans from −1.8 m to 5 m and the initial conditions are listed in Table 2.
A free-stream velocity of 6047 m·s 1 is applied with a free-stream pressure of 0.89 Pa and temperature T = 220 K. Such a temperature is high enough to result in a vibrationally-excited and chemically active flow-field. The cylinder wall is held at a uniform temperature of 1000 K. The overall Knudsen number of 0.0022 lies in the lower range of the continuum-transition regime; however, the gas locally may lie in the transition regime.
The mesh used in this investigation is shown in Figure 3 and was constructed using Ansys ICEM (Integrated Computer Engineering and Manufacturing) CFD (Canonsburg, PA, USA) [45]. This mesh consists of 155,000 cells with the first cell spacing at the cylinder wall set to 2 microns. Both Maxwellian velocity slip and Smoluchowski temperature jump boundary conditions were applied at the walls with the accommodation coefficient equal to 1.
For the shear viscosity and the thermal conductivity, the Blottner and Eucken formulae were applied, respectively, while the mixing rule employed was that of Armaly and Sutton. Both MWP and SSH formulations are successively used for V–T energy transfer for comparison purposes. The different set-ups are summarised in Table 3 together with a run identification number.
Two chemical reactions were also considered, these being the irreversible molecule–molecule and molecule–atom dissociation of nitrogen
N 2 + N 2 2 N + N 2 , N 2 + N 2 N + N .
The Arrhenius rate constants are shown in Table 4 in which the units of A and T a are given in m 3 ·mol 1 ·s 1 and Kelvin, respectively. They are derived from the QK theory [8] and Park’s rates for use in a two-temperature CFD solver [2].
Finally, two configurations associating a chemistry–vibration model with a set of chemical rates have been studied. The first one is called CVDV-QK and has been tested for zero-dimensional heat bath scenarios only [4]. The second configuration, subsequently named Park, is the one conventionally used in hypersonics and combines the Park TTv model with Park’s rates. The set-up for this case can be found in the Supplementary Materials.
The DSMC configuration is provided in [5], and good DSMC practice was satisfied for both the mesh and time-step. In the DSMC mesh, a total of 5.5 million cells were employed, which were filled with over 80 million equivalent DSMC particles at steady-state. Non-reacting case scenarios were also executed in [5], and the analysis will not be repeated here, although some of the results will be reported in the following graphs and tables.
The departure from local thermodynamic equilibrium is shown in Figure 4a for the CFD run number 2 using the local gradient-length Knudsen number. Dark grey and black colours indicate K n G L L values beyond 0.05 and cover, as expected, the bow shock and the near-wake areas. Similar results were found for the other runs.
The Mach and trans-rotational temperature fields given by hy2Foam and dsmcFoam in Figure 4b,c compare well in the compression region, for x < 0 . In the wake of the cylinder, however, T t r is generally smaller using DSMC. Unlike the non-reacting simulation that demonstrated a good agreement for the vibrational temperature field in the compression area [5], the discrepancies in Figure 4d are much larger this time in the whole domain when compared with DSMC. This can be explained by the application of the QK theory to capture the chemistry–vibration coupling and the use of the quantum Larsen–Borgnakke method in dsmcFoam, which promotes a quicker energy harmonisation in expansion regions (where T v > T t r ) as reported in [4] using a zero-dimensional analysis. The Park combination used in run 1 is the instance that provides less accurate results with local vibrational temperatures above 10,000 K.
Figure 5a–c compares the stagnation line profiles of Mach number, temperature, and number density given by hy2Foam and dsmcFoam. It is normally evident that the bow shock is more diffuse when using the DSMC method compared with CFD. However, it is evident that the shock stand-off distances are almost identical using both solvers and are approximately equal to 0.25 m, which is about 5 cm closer to the body than for the non-reacting case [5]. The peak in trans-rotational temperature is correctly determined using the CVDV-QK combination for both runs 1 and 2 and is slightly over-predicted using the Park combination. As shown previously in Figure 4d, the trend in vibrational temperature shows a steeper increase across the shock wave. In Figure 5c, the evolution of the species number densities for all runs is in satisfactory agreement with dsmcFoam outside of the K n G L L > 0.01 band. The early production of atomic nitrogen within the shock is not captured in the CFD solver due to the slight difference in shock thickness prediction.
Surface properties of pressure coefficient, skin friction and heat transfer are shown in Figure 5d–f for both non-reacting (NR) and reacting simulations. There is a reasonable agreement between the CFD and DSMC solvers for C p and C f .
The drag coefficient for each simulation is provided in Table 5, and this coefficient estimated by hy2Foam represents less than 2% error when compared with dsmcFoam, and the reacting environment does not significantly affect its magnitude. In addition, the values of hy2Foam and dsmcFoam are reasonably close to the one predicted by the Newtonian theory that is 4 / 3 . The integrated heat flux, C H , is however showing larger discrepancies between the two codes. There are several factors that could explain this observation: (1) a K n G L L number greater than 0.1 all around the cylinder; (2) a different treatment of the vibration-translational energy transfer between CFD and DSMC codes; and (3) the use of the KNP numerical schemes in hy2Foam that are known for being too dissipative in the near-wall region. It is also shown that the CVDV-QK association is producing a 27% larger integrated heat flux as compared with dsmcFoam, while the Park combination overpredicts C H by 39%.

4. Conclusions

The newly-coded open-source two-temperature CFD solver hy2Foam has been extended to simulate hypersonic multi-dimensional case scenarios. hy2Foam has shown to perform as well as a Navier–Stokes CFD solver from the University of Michigan for a Mach 11 blunted cone, and to be in satisfactory agreement with the dsmcFoam code for a Mach 20 flow of a binary reacting mixture around a circular cylinder. For this latter case, the aerothermodynamic loads were better estimated using the novel CVDV-QK model combination rather than the conventional Park combination, when compared to the dsmcFoam code. This result reaffirms the predictions found using zero-dimensional cases in [4]. Finally, it is considered that the cylinder case scenario presented in this paper provides a useful basis for other codes to compare against.
Future work for the extension of hy2Foam will include the development of an 11-species plasma model for application in weakly-ionized flow environments. The discrepancies observed between the CFD and DSMC codes for K n G L L numbers outwith the continuum regime will motivate the development of an open-source hypersonic hybrid hydrodynamic-molecular gas flow solver using both the dsmcFoam and hy2Foam codes [7].

Supplementary Materials

The following are available online at www.mdpi.com/2226-4310/3/4/45/s1.

Acknowledgments

Vincent Casseau wishes to thank the James Weir Fluids Laboratory (University of Strathclyde, Glasgow, UK) and the centre for Future Air-Space Transportation Technology (cFASTT) (University of Strathclyde, Glasgow, UK) who supported this work. Vincent Casseau would also like to thank Benzi John and David R. Emerson from the Science and Technology Facilities Council (STFC) Daresbury Laboratory in Warrington, UK for setting-up and running the DSMC computations and acknowledge Jimmy-John Hoste from cFASTT for his help in testing the hy2Foam solver. The CFD results were obtained using the Engineering and Physical Sciences Research Council (EPSRC) funded ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk/). EPSRC grant no. EP/K000586/1.

Author Contributions

Vincent Casseau implemented the hy2Foam solver; Vincent Casseau and Daniel E. R. Espinoza performed the computations; Vincent Casseau, Daniel E. R. Espinoza, Thomas J. Scanlon and Richard E. Brown analyzed the data; and Vincent Casseau and Thomas J. Scanlon wrote the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bird, G.A. Molecular Gas Dynamics and the Direct Simulation of Gas Flows; Clarendon: Oxford, UK, 1994. [Google Scholar]
  2. Park, C. Nonequilibrium Hypersonic Aerothermodynamics; Wiley International: New York, NY, USA, 1990. [Google Scholar]
  3. Github Repository of the hy2Foam Solver. Version 0.01, Commit 98032f1. Available online: https://github.com/vincentcasseau/hyStrath/ (accessed on 5 December 2016).
  4. Casseau, V.; Palharini, R.C.; Scanlon, T.J.; Brown, R.E. A Two-Temperature Open-Source CFD Model for Hypersonic Reacting Flows, Part One: Zero-dimensional Analysis. Aerospace 2016, 3, 34. [Google Scholar] [CrossRef]
  5. Casseau, V.; Scanlon, T.J.; John, B.; Emerson, D.R.; Brown, R.E. Hypersonic Simulations using Open-Source CFD and DSMC Solvers. In Proceedings of the AIP Conference, Victoria, BC, Canada, 10–15 July 2016; Volume 1786.
  6. OpenFOAM Official Website. Available online: http://www.openfoam.org/ (accessed on 2 August 2016).
  7. Espinoza, D.E.R.; Casseau, V.; Scanlon, T.J.; Brown, R.E. An Open Source Hybrid CFD-DSMC Solver for High-Speed Flows. In Proceedings of the AIP Conference, Victoria, BC, Canada, 10–15 July 2016; Volume 1786.
  8. Scanlon, T.J.; White, C.; Borg, M.K.; Palharini, R.C.; Farbar, E.; Boyd, I.D.; Reese, J.M.; Brown, R.E. Open source Direct Simulation Monte Carlo (DSMC) Chemistry Modelling for Hypersonic Flows. AIAA J. 2015, 53, 1670–1680. [Google Scholar] [CrossRef]
  9. Palharini, R.C. Atmospheric Reentry Modelling Using an Open-Source DSMC Code. Ph.D. Thesis, University of Strathclyde, Glasgow, UK, 2014. [Google Scholar]
  10. Borgnakke, C.; Larsen, P.S. Statistical collision model for simulating polyatomic gas with restricted energy exchange. In Rarefied Gas Dynamics; DFVLR Press: Porz-Wahn, Germany, 1974; Volume 1, p. A7. [Google Scholar]
  11. Bird, G.A. The DSMC Method, 2nd ed.; CreateSpace Independent Publishing Platform: Seattle, WA, USA, 2013. [Google Scholar]
  12. Candler, G.V.; Nompelis, I. Computational Fluid Dynamics for Atmospheric Entry. In Non-Equilibrium Dynamics: From Physical Models to Hypersonic Flights; The von Karman Institute for Fluid Dynamics: Rhode-Saint-Genèse, Belgium, 2009. [Google Scholar]
  13. Blottner, F.G.; Johnson, M.; Ellis, M. Chemically Reacting Viscous Flow Program for Multicomponent Gas Mixtures; Report No. SC-RR-70-754; Sandia Laboratories: Albuquerque, NM, USA, 1971. [Google Scholar]
  14. Vincenti, W.G.; Kruger, C.H. Introduction to Physical Gas Dynamics; Krieger Publishing Company: Malabar, FL, USA, 1965. [Google Scholar]
  15. Wilke, C.R. A Viscosity Equation for Gas Mixtures. J. Chem. Phys. 1950, 18, 517–519. [Google Scholar] [CrossRef]
  16. Gupta, R.; Yos, J.M.; Thompson, R.A.; Lee, K.P. A Review of Reaction Rates and Thermodynamic and Transport Properties for an 11-Species Air Model for Chemical and Thermal Nonequilbrium Calculations to 30,000 K; NASA RP-1232; NASA: Washington, DC, USA, 1990.
  17. Armaly, B.F.; Sutton, L. Viscosity of Multicomponent Partially Ionized Gas Mixtures. In Proceedings of the 15th Thermophysics Conference, Snowmass, CO, USA, 14–16 July 1980. AIAA Paper 80-1495.
  18. Palmer, G.E.; Wright, M.J. Comparison of Methods to Compute High-Temperature Gas Viscosity. J. Thermophys. Heat Transf. 2003, 17, 232–239. [Google Scholar] [CrossRef]
  19. Sutton, K.; Gnoffo, P.A. Multi-Component Diffusion with Application to Computational Aerothermodynamics. In Proceedings of the 7th AIAA/ASME Joint Thermophysics and Heat Transfer Conference, Albuquerque, NM, USA, 15–18 June 1998. AIAA Paper 98-2575.
  20. Hirschfelder, J.O.; Curtiss, C.; Bird, R.B. Molecular Theory of Gases and Liquids; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 1954. [Google Scholar]
  21. Alkandry, H.; Boyd, I.D.; Martin, A. Comparaison of Models for Mixture Transport Properties for Numerical Simulations of Ablative Heat-Shields. In Proceedings of the 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition Grapevine (Dallas/Ft. Worth Region), Grapevine, TX, USA, 7–10 January 2013. AIAA Paper 2013-0303.
  22. Landau, L.; Teller, E. On the theory of sound dispersion. Phys. Z. Sowjetunion 1936, 10, 34. [Google Scholar]
  23. Millikan, R.C.; White, D.R. Systematics of Vibrational Relaxation. J. Chem. Phys. 1963, 39, 3209–3213. [Google Scholar] [CrossRef]
  24. Schwartz, R.N.; Slawsky, Z.I.; Herzfeld, K.F. Calculation of Vibrational Relaxation Times in Gases. J. Chem. Phys. 1952, 20, 1591–1599. [Google Scholar] [CrossRef]
  25. Thivet, F.; Perrin, M.Y.; Candel, S. A unified nonequilibrium model for hypersonic flows. Phys. Fluids 1991, 3, 2799–2812. [Google Scholar] [CrossRef]
  26. Knab, O.; Frühauf, H.H.; Jonas, S. Multiple Temperature Descriptions of Reaction Rate Constants with Regard to Consistent Chemical-Vibrational Coupling. In Proceedings of the 27th Thermophysics Conference, Nashville, TN, USA, 6–8 July 1992. AIAA Paper 92-2947.
  27. Knab, O.; Frühauf, H.H.; Messerschmid, E.W. Theory and Validation of the Physically Consistent Coupled Vibration-Chemistry-Vibration Model. J. Thermophys. Heat Transf. 1995, 9, 219–226. [Google Scholar] [CrossRef]
  28. Gnoffo, P.A.; Gupta, R.; Shinn, J.L. Conservation Equations and Physical Models for Hypersonic Air Flows in Thermal and Chemical Nonequilibrium; Technical Report NASA-TP-2867; NASA Langley: Hampton, VA, USA, 1989.
  29. Boyd, I.D.; Chen, G.; Candler, G.V. Predicting failure of the continuum fluid equations in transitional hypersonic flows. Phys. Fluids 1995, 7, 210–219. [Google Scholar] [CrossRef]
  30. Schwartzentruber, T.E.; Scalabrin, L.C.; Boyd, I.D. A modular particle-continuum numerical method for hypersonic non-equilibrium gas flows. J. Comput. Phys. 2007, 225, 1159–1174. [Google Scholar] [CrossRef]
  31. Abbate, G.; Thijsse, B.J.; Kleijn, C.R. Computational Science—ICCS 2007. Part I, Chapter Coupled Navier–Stokes/DSMC Method for Transient and Steady-State Gas Flows. In Proceedings of the 7th International Conference, Beijing, China, 27–30 May 2007; Springer: Berlin/Heidelberg, Germany, 2007; Volume 225, pp. 842–849. [Google Scholar]
  32. Darbandi, M.; Roohi, E. Applying a hybrid DSMC/Navier–Stokes frame to explore the effect of splitter catalyst plates in micro/nanopropulsion systems. Sens. Actuators A 2013, 189, 409–419. [Google Scholar] [CrossRef]
  33. Roveda, R.; Goldstein, D.B.; Varghese, P.L. Hybrid Euler/particle approach for continuum/rarefied flows. J. Spacecr. Rocket. 1998, 35, 258–265. [Google Scholar] [CrossRef]
  34. Schwartzentruber, T.E.; Scalabrin, L.C.; Boyd, I.D. Hybrid Particle-Continuum Simulations of Nonequilibrium Hypersonic Blunt-Body Flowfields. J. Thermodyn. Heat Transf. 2008, 22, 29–37. [Google Scholar] [CrossRef]
  35. Von Smoluchowski, M. Über wärmeleitung in verdünnten Gasen. Ann. Phys. Chem. 1898, 64, 101–130. [Google Scholar] [CrossRef]
  36. Maxwell, J.C. On stresses in Rarefied Gases Arising from Inequalities of Temperature. Phil. Trans. R. Soc. Lond. 1879, 170, 231–256. [Google Scholar] [CrossRef]
  37. Greenshields, C.J.; Weller, H.G.; Gasparini, L.; Reese, J.M. Implementation of semi-discrete, non-staggered central schemes in a collocated, polyhedral, finite volume framework, for high-speed viscous flows. Int. J. Numer. Methods Fluids 2010, 63, 1–21. [Google Scholar]
  38. Kurganov, A.; Noelle, S.; Petrova, G. Semi-discrete central-upwind schemes for hypersonic conservation laws and Hamilton-Jacobi equations. SIAM J. Sci. Comput. 2001, 23, 707–740. [Google Scholar] [CrossRef]
  39. ARCHIE-WeSt High Performance Computer. Available online: https://www.archie-west.ac.uk/ (accessed on 20 November 2016).
  40. Wang, W.L.; Boyd, I. Hybrid DSMC-CFD Simulations of Hypersonic Flow over Sharp and Blunted Bodies. In Proceedings of the 36th AIAA Thermophysics Conference, Orlando, FL, USA, 23–26 June 2003. AIAA Paper 2003-3644.
  41. Holden, M.S. Experimental Database from CUBRC Studies in Hypersonic Laminar and Turbulent Interacting Flows including Flowfield Chemistry; RTO Code Validation of DSMC and Navier–Stokes Code Validation Studies CUBRC Report; Calspan-University at Buffalo Research Center: Buffalo, NY, USA, 2000. [Google Scholar]
  42. Geuzaine, C.; Remacle, J.F. Gmsh: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 2009, 79, 1309–1331. [Google Scholar] [CrossRef]
  43. Ahrens, J.; Geveci, B.; Law, C. ParaView: An End-User Tool for Large Data Visualization; Elsevier: Amsterdam, The Netherlands, 2005. [Google Scholar]
  44. Williams, T.; Kelley, C. Gnuplot 5.0: An Interactive Plotting Program. Available online: http://gnuplot.sourceforge.net/ (accessed on 20 November 2016).
  45. ANSYS, Inc. ANSYS ICEM CFD Academic Research, Release 15.0. Available online: http://www.ansys.com (accessed on 20 November 2016).
Figure 1. Mesh for the blunted cone (each 5th line is represented in each direction).
Figure 1. Mesh for the blunted cone (each 5th line is represented in each direction).
Aerospace 03 00045 g001
Figure 2. (ac) Stagnation line profiles, and (df) surface quantities along the blunted cone; (a) normalised temperature; (b) normalised mass density; (c) normalised velocity; (d) pressure coefficient; (e) friction coefficient; (f) Stanton number.
Figure 2. (ac) Stagnation line profiles, and (df) surface quantities along the blunted cone; (a) normalised temperature; (b) normalised mass density; (c) normalised velocity; (d) pressure coefficient; (e) friction coefficient; (f) Stanton number.
Aerospace 03 00045 g002
Figure 3. Mesh for the cylinder (each 5th line is represented in each direction).
Figure 3. Mesh for the cylinder (each 5th line is represented in each direction).
Aerospace 03 00045 g003
Figure 4. Computational fluid dynamics-direct simulation Monte Carlo (CFD-DSMC) flow-field comparisons for run number 2. In (bd): the dsmcFoam solution is represented in the upper half and the hy2Foam solution in the lower half. (a) Local gradient-length Knudsen number; (b) Mach number; (c) trans-rotational temperature; and (d) vibrational temperature.
Figure 4. Computational fluid dynamics-direct simulation Monte Carlo (CFD-DSMC) flow-field comparisons for run number 2. In (bd): the dsmcFoam solution is represented in the upper half and the hy2Foam solution in the lower half. (a) Local gradient-length Knudsen number; (b) Mach number; (c) trans-rotational temperature; and (d) vibrational temperature.
Aerospace 03 00045 g004
Figure 5. (ac): Stagnation line profiles. CFD run 1: black lines, run 2: red lines, run 3: blue lines. (df): surface quantities around the cylinder. (a) Mach number; (b) temperature; (c) number density; (d) pressure coefficient; (e) skin-friction coefficient; (f) surface heat flux.
Figure 5. (ac): Stagnation line profiles. CFD run 1: black lines, run 2: red lines, run 3: blue lines. (df): surface quantities around the cylinder. (a) Mach number; (b) temperature; (c) number density; (d) pressure coefficient; (e) skin-friction coefficient; (f) surface heat flux.
Aerospace 03 00045 g005
Table 1. Initial conditions for the Mach 11.3 blunted cone.
Table 1. Initial conditions for the Mach 11.3 blunted cone.
QuantityValueUnit
Free-stream velocity, U 2764.5m/s
Free-stream pressure, p 21.9139Pa
Free-stream density, ρ 5.113 × 10 4 kg/m 3
Free-stream temperature, T 144.4K
Free-stream mean-free-path, λ 1.01 × 10 4 m
Overall Knudsen number, K n o v 0.002-
Wall temperature, T w 297.2K
Table 2. Initial conditions for the Mach 20 cylinder.
Table 2. Initial conditions for the Mach 20 cylinder.
QuantityValueUnit
Free-stream velocity, U 6047m/s
Free-stream pressure, p 0.89Pa
Free-stream density, ρ 1.363 × 10 5 kg/m 3
Free-stream temperature, T 220K
Free-stream mean-free-path, λ 4.45 × 10 3 m
Overall Knudsen number, K n o v 0.0022-
Wall temperature, T w 1000K
Table 3. Computational fluid dynamics (CFD) simulations performed.
Table 3. Computational fluid dynamics (CFD) simulations performed.
Run NumberV–T TransferElectronic ModeCV ModelRates
1MWPnoCVDVQK
2SSHnoCVDVQK
3MWPnoPark TTvPark
Table 4. Parameters for the evaluation of the forward rate constant.
Table 4. Parameters for the evaluation of the forward rate constant.
Reaction RateReactionArrhenius Law Constants
Colliding PartnerAβ T a
Quantum-Kinetics (QK)N2 2.47 × 10 18 −0.62113,500
N 6.02 × 10 18 −0.68113,500 
ParkN2 7.0 × 10 21 −1.6113,200
N 3.0 × 10 22 −1.6113,200
Table 5. Aerothermodynamic coefficients.
Table 5. Aerothermodynamic coefficients.
CFD Run Number C D C H ( kW )
CFDDSMCCFDDSMC
NR1.31.286106115
11.3021.28481.063.3
21.30280.5
31.30488.1
C D is the drag coefficient and C H represents the integrated heat flux.

Share and Cite

MDPI and ACS Style

Casseau, V.; Espinoza, D.E.R.; Scanlon, T.J.; Brown, R.E. A Two-Temperature Open-Source CFD Model for Hypersonic Reacting Flows, Part Two: Multi-Dimensional Analysis. Aerospace 2016, 3, 45. https://doi.org/10.3390/aerospace3040045

AMA Style

Casseau V, Espinoza DER, Scanlon TJ, Brown RE. A Two-Temperature Open-Source CFD Model for Hypersonic Reacting Flows, Part Two: Multi-Dimensional Analysis. Aerospace. 2016; 3(4):45. https://doi.org/10.3390/aerospace3040045

Chicago/Turabian Style

Casseau, Vincent, Daniel E. R. Espinoza, Thomas J. Scanlon, and Richard E. Brown. 2016. "A Two-Temperature Open-Source CFD Model for Hypersonic Reacting Flows, Part Two: Multi-Dimensional Analysis" Aerospace 3, no. 4: 45. https://doi.org/10.3390/aerospace3040045

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop