Extension of a Hybrid Particle-Continuum Method

for a Mixture of Chemical Species


Analysis of hypersonic flows requires consideration of multiscale phenomena due to the range of flight regimes encountered. At high altitudes, for example, low densities lead to rarefied flows that are characterized by a large global Knudsen number, Kn, which is defined as the ratio of the freestream mean free path of a gas to a characteristic length scale of the vehicle or phenomenon of interest. However, even at low altitudes where this global parameter may be small, localized regions of thermodynamic nonequilibrium may still exist. Flow gradients that occur over length scales of the order of a mean free path, such as those that arise in strong shock waves, boundary and shear layers, and low density wakes, lead to regions where an equilibrium velocity distribution function (VDF) cannot be maintained due to an insufficient number of collisions among gas particles. These flow phenomena are apparent in the shadowgraph image shown below of a NASA manned capsule concept from 1957. In addition, the high temperatures encountered in these hypersonic flows can lead to excitation of internal energy modes, dissociation, ionization, and chemical reactions among constituent gas species, the finite rates of which may become important due to inherently small residence times of the high-speed gas.

Figure 1: Manned Capsule (NASA)

A computational fluid dynamics (CFD) simulation approach based on the Navier-Stokes equations is computationally efficient and physically accurate for near-continuum flows, but accurate simulation of the aforementioned nonequilibrium effects requires kinetic theory-based techniques, such as the direct simulation Monte Carlo (DSMC) method. Due to restrictions on cell size and time step, although the DSMC method may be physically accurate throughout the entire flow field, it becomes computationally demanding for high-density flows. As such, hybrid particle-continuum methods are an attractive alternative for transitional, hypersonic flows, where a CFD solver is employed where continuum assumptions are applicable, and the DSMC method is utilized in regions delineated by a continuum breakdown parameter. Because of the combined numerical efficiency and physical accuracy of such hybrid particle-continuum methods, the present research efforts aim to extend the capabilities of an existing Modular Particle-Continuum (MPC) method to include multiple chemical species. Such a development not only broadens the applicability of this MPC method to more realistic atmospheric flows, but also serves as a preliminary step towards the method's ability to accurately simulate thermochemical nonequilibrium effects in the future.

Overview of Modular Particle-Continuum Method

The MPC method that has been augmented in this study is a hybrid particle-continuum code that couples the LeMANS CFD code, which was developed at the University of Michigan, to the DSMC code, MONACO. The physical accuracy and numerical efficiency of the MPC code when compared to both full CFD and full DSMC simulation results have been demonstrated for one-dimensional shock waves and for axisymmetric and two-dimensional, steady, hypersonic blunt-body flows. The LeMANS code itself is a finite volume formulation capable of numerically solving the laminar Navier-Stokes equations on axisymmetric, two-dimensional, and three-dimensional grids. Vibrational and rotational nonequilibrium effects can be simulated, along with multispecies gases, weakly ionized flows, and finite rate chemistry. The MONACO code is also capable of simulating rotational and vibrational nonequilibrium effects, in addition to finite rate chemistry. Recent developments in the MPC code have taken advantage of these capabilities by extending the hybrid framework to include rotational and vibrational nonequilibrium effects for a single chemical species. The MPC code is also parallelized using Message Passing Interface (MPI) and METIS graph partitioning.

Beginning with a converged CFD solution from LeMANS, the MPC method calculates the gradient-length local Knudsen number, λ/Q|∇Q|, throughout the domain. This is the continuum breakdown parameter employed to demarcate regions where the continuum assumptions of the Navier-Stokes equations are valid from those regions in which the DSMC method must be employed for physical accuracy. In the above equation, Q represents a macroscopic fluid property, such as mass density, translational temperature, rotational temperature, or speed, and λ is the local mean free path. Below a continuum breakdown threshold value of 0.05, the flow is considered to be in a state of thermodynamic equilibrium, allowing for no more than a 5% difference between results obtained with a continuum method and a particle method.

Once initial interfaces between the particle and continuum domains are determined, the particle domain is extended and populated with DSMC simulator particles consistent with a Chapman-Enskog VDF described by the macroscopic properties as calculated by the CFD solver. Particle domains are extended into the continuum region by a specified number of "buffer" cells in order to allow for adjustment of the hybrid interfaces and provide more accurate boundary conditions to the DSMC module, although the ability of the MPC method to handle ill-posed boundary conditions has been demonstrated. State-based coupling procedures are used to transfer information between the CFD and DSMC code modules such that existing boundary condition routines are essentially unmodified. The hybrid interface geometry and terminology are illustrated in the figure below.

Figure 2: Illustration of hybrid interface

After a specified number of DSMC time steps, when the number of newly created DSMC cells at each evaluation of the continuum breakdown parameter has decreased below a reasonable level, information is transferred to the CFD domain. In order to mitigate the inherent statistical scatter of the DSMC method when updating the CFD boundary conditions, the MPC code employs a combination of spatial and temporal averaging. With these boundary conditions from the DSMC module, the CFD module then proceeds to perform a specified number of iterations, after which information is passed back to the DSMC module in the form of representative simulator particles. This reciprocal information exchange between the two domains continues, with recurring calculation of the continuum breakdown parameter, for a sufficient number of iterations until the hybrid interfaces remain stationary. At this point in the hybrid procedure, the interfaces are locked and standard DSMC sampling and CFD time-stepping commence.

Recent Results

Recently, results have been published for a Mach 10, Kn 0.01 flow of 50% N2/50% N by mole over a two-dimensional cylinder. The freestream temperature is 217.45 K, the freestream velocity is 3594 m/s, and the cylinder surface is specified to maintain a constant temperature of 1000 K. As it is the objective of the MPC method to obtain the same level of physical accuracy as the DSMC method at a fraction of the computational cost, the MPC method is compared to full CFD and full DSMC simulations in terms of flow field and surface properties, as well as computational performance metrics. Shown in the figure below are the initial hybrid interfaces obtained by applying the aforementioned continuum breakdown parameter to the full CFD solution. In addition, the final interfaces are also shown. As expected, by employing the DSMC method to simulate the strong bow shock, boundary layer, and wake region of this flow field, and transferring more accurate boundary condition information back to those regions in which CFD is used, the nonequilibrium regions surrounding the shock and wake are seen to expand in size as the MPC method proceeds. This results in a thicker shock, which the Navier-Stokes equations alone are unable to predict, and corrected wake region.

Figure 3: Hybrid interfaces

These features are evident in the next figure, where the contours of velocity magnitude given by the MPC method are in agreement with the predictions of the DSMC method, and show that the onset of the shock is upstream of that predicted by the full CFD solution. It is also shown that the MPC method is able to accurately reproduce the wake region size obtained with full DSMC. Because full DSMC and the DSMC module of the MPC method, which is being used in the boundary layer and wake region, allow for velocity slip at the cylinder surface, whereas the no-slip condition is enforced in the full CFD simulation, the wake predicted by CFD is larger than that predicted by the DSMC and MPC methods. These differences are also seen in the contours of translational temperature predicted by the MPC method, full DSMC, and full CFD.

(a) Contours of velocity magnitude

(b) Contours of translational temperature

Figure 4: 2D cylinder flow field properties

Of particular interest in these simulations involving a multispecies gas is the concentration of constituent chemical species throughout the flow field. The figure below shows a comparison of the N2 mole fraction along the stagnation streamline given by full DSMC, full CFD, and the MPC method. The full CFD result gives a nearly constant mole fraction consistent with the freestream composition along the entire streamline, except for a slight decrease in the N2 mole fraction in the shock and a slight increase at the cylinder surface. Because these regions have been labeled as particle domains by the MPC method due to continuum breakdown parameter values in excess of 0.05, it is expected that such phenomena cannot be properly simulated with CFD alone. As such, both the full DSMC results and those from the MPC method indicate a marked decrease in the N2 mole fraction in the shock, a higher concentration of N2 relative to N between the shock and boundary layer, and a pronounced increase in the N2 mole fraction in the boundary layer. It is important to notice the very close agreement between the MPC method and full DSMC in the fore body flow, specifically, the region between the shock and boundary layer in which the CFD module of the MPC method is employed. Even though this is labeled as part of the continuum domain, improved boundary conditions provided by the DSMC module of the MPC method being used in surrounding regions have enabled the CFD module to yield results that are in much better agreement with full DSMC than could have been achieved with CFD alone.

Figure 5: Molecular nitrogen mole fraction along stagnation streamline

Perhaps the most important comparison of this analysis is that of the aerothermodynamic surface properties along the two-dimensional cylinder, as these determine the characteristics of the thermal protection system, control surfaces, and structure of a hypersonic vehicle. The agreement between full DSMC and the MPC method in terms of the heat transfer coefficient is evident on the left of the following figure. Unlike these two sets of results, the temperatures near the cylinder in the CFD simulation are forced to match the isothermal temperature of the wall, thus causing larger temperature gradients and a higher heat transfer along the entire surface of the cylinder. CFD also over predicts the surface pressure coefficient relative to both full DSMC and the MPC method after an angular position of approximately 90 degrees, as shown on the right of this figure.

Figure 6: Heat transfer coefficient (left axis) and Surface pressure coefficient (right axis)

The larger wake region of the CFD solution predicted by the contours of velocity magnitude is reaffirmed by the shear stress profile shown in the figure below. Here the change in sign of the shear stress coefficient, as indicated by a sharp decrease in the profile, occurs at an angular position of 159 degrees, as opposed to the full DSMC and MPC method profiles, in which the recirculation zone begins at an angular position of 163 degrees. In addition, the CFD profile of shear stress coefficient is generally higher than those of the full DSMC and MPC methods along the majority of the cylinder surface. Again, this is due to the no-slip condition enforced at the wall in CFD, but not in the other two simulation techniques.

Figure 7: Shear stress coefficient

Although the MPC method has been shown to reproduce full DSMC results in terms of both flow field and surface properties, it is nearly as important that the MPC method obtain this level of accuracy at a reduced computational cost. To assess the computational requirements, simulations employing full CFD, full DSMC, and the MPC method were run on the Pleiades supercomputer maintained by the NASA Advanced Supercomputing Division. The number of Nehalem processors requested for each simulation was tailored to each technique so that an optimal balance between parallel efficiency and wall time was maintained. The cost of a DSMC simulation scales nearly linearly with the number of simulator particles. In the present study, the DSMC simulation requires 41.2 million particles, whereas the particle domains of the MPC simulation require 20.1 million particles. The same grid and numerical weights are used in both of these cases. The ideal speedup, as calculated by the ratio of particles required by DSMC to the number of particles required by MPC, is then 2.05. Thus, the ideal speedup assumes that the computational requirements of the CFD module of the MPC method are negligible. The actual speedup, which is the ratio of computational time required to complete the DSMC simulation (approximately 1100 CPU hours) to the time required by the MPC simulation (approximately 560 CPU hours), is calculated to be 1.95. It is important to remember that the same number of time steps is used to sample the particle information in each simulation. However, due to the fact that the MPC simulation begins with a converged CFD solution, many more time steps are usually required in the full DSMC simulation prior to sampling than are required by the MPC method. In the present MPC simulation, approximately 20,000 time steps are required for the hybrid interfaces to adjust to their final locations, whereas 110,000 time steps are required for the number of simulator particles in the full DSMC simulation to reach a steady-state value prior to sampling. In addition, the MPC method uses only 85% of the memory required by the full DSMC simulation for this set of flow conditions.


This research was made possible through a National Defense Science and Engineering Graduate (NDSEG) Fellowship provided by the Department of Defense. In addition, we gratefully acknowledge the use of computational resources through the NASA Advanced Supercomputing Division.

Recent Publications