Development of a Hybrid Modular Particle-Continuum (MPC) Method for Nonequilibrium, Hypersonic Flows


Computation of the aerothermodynamics of hypersonic re-entry vehicles along their entire trajectory involves near-continuum conditions at low altitudes and rarefied, or non-equilibrium, conditions at high altitudes. Well-established simulation methods already exist for each of these flow regimes. For example, the continuum Navier-Stokes (NS) equations may be solved using algorithms from Computational Fluid Dynamics (CFD) and the particle-based, direct simulation Monte Carlo method (DSMC) may be used for the non-equilibrium, or rarefied, flows. A hybrid method that blends the CFD and DSMC techniques is an attractive approach for simulation of flows involving a mixture of both continuum and non-equilibrium flow regimes. This is the situation at intermediate altitudes where within a mostly continuum flow, there may be local regions of non-equilibrium flow generated by both the rapid expansion behid a re-entering capsule as well as by strong gradients in shock waves and boundary layers.

 Blunt Body Hypersonic Flow

For dilute gases, such as the earth's atmosphere, the most popular numerical method for simulating non-equilibrium flow is the DSMC method first developed by Graeme Bird. Here, the trajectories of a large number of simulated particles are followed simultaneously through a grid of computational cells. For each iteration, particles are first moved along their trajectories without colliding, after which particles residing within the same cell are randomly selected for a collision process. A major limitation of DSMC is that in order for this collision process to be physically accurate, the cell size must be on the order of the mean free path, λ, while containing at least 20 simulated particles per cell. As a result, 2D and 3D DSMC simulations can require prohibitively high numbers of computational cells and therefore simulated particles, especially in regions where λ is very small. This can be seen in variation in mean free path around a planetary probe (seen below) where the characteristic flow length scale can be very small in the dense fore body region of the probe. It is precisely in these regions, where λ is small, that the continuum approximation is valid and the NS equations can be solved without the same restriction on cell size. This is the motivation behind developing a hybrid particle-continuum numerical scheme for hypersonic non-equilibrium flows.

Basic Hybrid Methodology

A hybrid CFD-DSMC approach, developed at the University of Michigan, is called the Modular Particle Continuum (MPC) method and merges existing CFD and DSMC codes.  Both codes remain virtually unmodified and hybrid data structures and functions are constructed to track macroscopic quantities, passs information between solvers, and perform other required hybrid processes. The accuracy of the MPC method relies on the proper positioning of the DSMC-CFD interfaces. The interface must lie in near equilibrium regions where solution of the NS equations will introduce minimal or no error. Typically, particle and continuum regions are determined by evaluating a continuum breakdown parameter in the flow field. The MPC method uses the gradient-length Knudsen number, KnGL-Q|∇Q/Q|, where Q represents local flow quantities of interest such as density, temperature, or velocity magnitude, λ is the local mean-free-path. It has been shown for hypersonic flows that, in regions of the flow field where KnGL-MAX <0.05, the discrepancy between a NS and DSMC solution is less than 5%. Thus, these regions could be solved using a NS solver with little error.

The MPC method begins with a NS solution of the entire flow field and then uses gradient-length Knudsen number to decompose the domain into CFD and DSMC regions. The method uses state-based coupling to transfer information between particle and continuum regions. After evaluation of the breakdown parameter, the particle region is extended by a few extra cells into the continuum domain to create an overlap region as seen in the figure below. Next, one row of NS and two rows of DSMC boundary cells are initialized. The DSMC and NS domains are then coupled by transferring information across the interfaces. With state-based coupling, this involves updating the boundary conditions of each solver. In this way, information transfer into both the particle and continuum regions is handled through existing boundary procedures already used by both solvers. The DSMC boundary cells are continually filled with particles consistent with the flow properties in the corresponding NS cell using the Chapman-Enskog distribution based on the local macroscopic state and gradients, known from the NS solver. As particles in the DSMC domain interact, their distributions evolve in time towards the distribution predicted by full DSMC.  The MPC method also tracks the macroscopic variation in each DSMC cell. In order to provide these averaged properties with minimal statistical scatter, a mixture of spatial and temporal averaging is used. Specifically, MPC uses the subrelaxation technique proposed by Sun et al. These averaged DSMC properties are then used to update boundary conditions for the NS solver.



CFD-DSMC Mesh Structure

The DSMC-CFD interface locations for Mach 20 flow over a sting-mounted planetary probe with a global Knudsen number of 0.01 are shown below along with the variation in flow length scales.  Due to the overlap regions, the interface locations have moved from the initial position that was calculated using a full Navier-Stokes solution to the final interface location based on the physically correct, DSMC solution.  The MPC method is able to simulate steady-state, hypersonic flow that is in thermal (rotational and vibrational) nonequilibrium.  It is parallelized for distributed memory systems using dynamic domain decomposition with scaled efficiency remaining above 80% for typical near-continuum 2D and axi-symmetric flows.

Continuum-Rarefied Interface 

Recent Results

The following figures compare translational and rotational temperature contours predicted by full DSMC, full CFD, and the MPC method for Mach 20 flow over a sting-mounted, planetary probe with a global Knudsen number of 0.01. In general the MPC method has improved agreement with full DSMC results from the initial CFD solution across the entire flow field. Even in regions where the CFD module is used, the MPC method has obtained improved boundary conditions from DSMC, which has allowed the CFD module to shift its result to excellent agreement with full DSMC.


In addition to improving agreement with DSMC for macroscopic quantities, the MPC method can accurately reproduce full DSMC results down to the velocity and energy probability density functions.  The figures below compare velocity and rotational energy probability density functions predicted by full DSMC, full CFD, and the MPC method at point A shown in the figure above. Due to the high degree of collisional nonequilibrium within the shock, the CFD velocity distribution function, which is generated from gradients and the first order Chapman-Enskog expansion, does not contain sufficient information to correctly generate the velocity distribution function predicted by full DSMC. In contrast, the MPC method is able to remain in very good agreement with DSMC throughout the entire velocity space. Despite the macroscopic rotational temperature predicted by CFD to be within 5% of the full DSMC result, the rotational energy distribution function predicted by full CFD is in poor agreement with the DSMC result throughout the entire rotational energy space. Similarly to the velocity distribution function, the MPC method remains in excellent agreement with full DSMC results for the rotational energy probability density function.

The MPC method is also in near exact agreement with full DSMC for surface properties.  The following figure shows surface heat transfer to the planetary probe predicted by full DSMC, full CFD, and the MPC method compared with available experimental measurements.  Along the fore body where the flow is highly collisional, all three methods are in good agreement with each other and the experimental measurements. Despite this highly collisional flow, CFD still slightly over predicts DSMC and MPC results. This is due to the inability to correctly model the Knudsen layer within the CFD solver. As the flow expands around the corner, full CFD over predicts both DSMC and experimental measurements by over an order of magnitude. In contrast, the MPC method remains in good agreement with both the experimental measurements and the full DSMC results. Similarly, the MPC method remains in excellent agreement with full DSMC and experimental measurements along the sting mount, while full CFD over predicts full DSMC along the entire sting mount.

Heat Transfer SR3_C2

In addition to reproducing full DSMC results, an objective of the MPC method is to reproduce DSMC at a reduced compuational cost compared to full DSMC.  All flows considered have experienced a speedup compared to full DSMC.  Speedup factors have ranged from 1.4 times faster for mostly rarefied flow to up to 30 times faster than DSMC for 2D and axi-symmetric, near-continuum flows, while still maintaining excellent agreement with full DSMC results.  It is expected that current work of extending the method to simulate full 3 dimensional flow will experience even larger differentials for compuational cost between the MPC method and full DSMC.

Recent Publications


Financial support for this work is provided by the Constellation University Institute Project (CUIP), under NASA grant NCC3-989.