Breakdown Parameters

The first issue in developing the coupled DSMC-NS method is how to determine both

the appropriate computational domain for the DSMC and NS solvers, and the proper interface boundary between these two solvers. A continuum breakdown parameter, proposed by Wang and Boyd [2003] for hypersonic flows, is employed in the present coupled DSMC-NS method as one of the criteria for selecting proper solvers. The continuum breakdown parameter Knmax is defined as [Wang and Boyd, 2003],

] temperature, respectively. They can be calculated from the following general formula [Wang and Boyd, 2003]

Q Q Kn_Q = λ ∇

(2.17)

where Q is the specific flow property (density, velocity and temperature) and λ is the local mean free path. If the calculated value of the continuum breakdown parameter in a region is larger than a preset threshold value, for example Kn_max^Thr., then it cannot be modeled using the NS equation. Instead, a particle solver like DSMC has to be used for that region.

In addition, another breakdown parameter is used to identify regions that exhibit thermal non-equilibrium among various degrees of freedom. The breakdown parameter of thermal equilibrium is proposed to be the ratio of the difference between translational and rotational temperatures to translational temperature in the present study, if diatomic gas molecules at

moderate temperature are involved in the simulation. Indeed, the definition of this parameter for atomic gas can be easily changed to the ratio of the difference between two translational temperatures to any specific translational temperature. For high-temperature flows, vibrational degrees of freedom may be also used to define this thermal non-equilibrium indicator. In the current study, this thermal non-equilibrium indicator is defined as

Tr

where TTr and TR are translational and rotational temperature, respectively. It is obvious that the nearer the value of PTne is zero, the closer the thermal equilibrium between translational and rotational degrees of freedom is. If the value of the computed thermal non-equilibrium indicator in a region is larger than some preset threshold value, for example

. Thr

PTne in the current study, then this flow region cannot be modeled correctly by the NS equation because it generally assumes thermal equilibrium among various degrees of freedom.

Hence, the DSMC method has to be used for that region instead. Note the parameter PTne

defined in Eq. (3) can be calculated only by DSMC, because HYB3D does not have the multi-temperature modelling capbability. This means it can only serve simultaneously with the Knudsen numbers (Eq. (1)) to control a potential switch back from DSMC to NS for those regions where the thermal non-equilibrium effect disappears after iterating between DSMC and NS methods.

Some comments address the use of the single-temperature Maxwellian distribution function at the interface of DSMC and NS regions in the proposed coupled method as follows.

Use of multiple-temperature NS equation solver could help to model the thermal non-equilibrium in much extended (into DSMC region) NS region; however, two or more energy (or temperature) equations have to be solved in all NS regions. In addition, extended NS region into DSMC region (larger property gradients or more rarefied) may further deteriorate the inherent assumption of Maxwellian distribution function in each degree of freedom in the multiple-temperature NS equation solver. Thus, use of two- or multiple-temperature NS solver may possibly benefit in reducing computational efforts, while it may introduce inaccuracy from physical point of view.

Furthermore, Chapmann-Enskog distribution temperature and density gradients may have to be used at the interface of DMSC and NS regions, where slight thermal non-equilibrium is considered. Its use would, however, greatly increase the computational cost as demonstrated in Garcia and Alder [1998] since more random number calls and computational operations are required. In practice, use of the single-temperature Maxwellian distribution function at the interface of DSMC and NS regions is much easier and with less computational cost. By taking the above intertwining factors into account, the proposed coupled DSMC-NS scheme simply utilizes Maxwellian distribution function at the DSMC-NS interface by properly controlling the magnitude of the breakdown parameters and

the appropriate overlapping regions that extends the DSMC region.

Based on the breakdown parameters, calculated from the preliminary simulation data using the continuum flow solver, and the criteria for the breakdowns of the continuum theory and thermal equilibrium, the domain for suitable DSMC simulation can be determined properly. Detailed procedures of determining the boundary (Boundary-I) between the

DSMC and NS approaches and marking the breakdown domain Ω^A are shown in Algorithm 2.1, which will be explained later. In addition, the mesh resolution across this region (e.g., shock layer) can be increased using mesh refinement approach [Wu et al., 2004] in the PDSC, although it is not employed in the current study for simplicity.

In Algorithm 2.1, distribution of breakdown parameters and the array of right-hand and left-hand cells for each cell interface are first read in. Note that the idea of right-hand and left-hand arrays of each cell interface is schematically shown in Algorithm 2.1. Then, all cells and cell faces are initialized, respectively, to be neither a continuum breakdown region

(Ω^A) nor a Boundary-I face. Note that the notations used in the current study can be found in Fig 2.3 with explanation that will be introduced shortly. All cells and cell faces are then checked, respectively, to decide if they are part of Ω^Aand interface Boundary-I. Using this subroutine, breakdown (Ω^A) and non-breakdown regions (ΩB ∪ΩC ∪ΩD), and Boundary-I can be properly identified for the entire computational domain.