A S ingle Cell with 5 S pecies with 34 Reactions

In this section, simulations with dissociation, recombination and exchange reactions (34 chemical reactions without ionizations) of a single cell are presented. The flow conditions are as follows; the initial total number density is 0.01 of standard atmospheric density and the initial mole fractions of nitrogen (N₂) and oxygen molecule

(O₂) are 0.8 and 0.2, respectively. The total number of particles is 100,000. The cell size of each simulation depends on the number density, which the ratio of cell size to mean free path is about 0.2. The order of time-step is 2.E-8. The transient time-step before sampling and the total time-step are 2,900,000 and 3,000,000, respectively. The real time is 0.06 seconds. The temperature range of these simulations is 2,000 to 16,000 K.

The original and new fitting input files (chem.dat) are simulated for comparison. For the new fitting rate constants of O₂ +N →NO+O reaction is

/kT) 10 exp(-4.970 T

10

3.91× ^{-19 0.65279} × ^-20

f =

k . Tables 6.24 are the mole percentage of

each species with original and new fitting rate constant. The mole percent is the number of moles of each species as a percentage of the total number of moles, which is also the percentage of number of particles.

Figure 6.25 shows the detailed composition as a function of temperature for each species. Figures 6.25(a) and 6.25(b) are original and new fitting input data, respectively.

The symbols in these figures are the steady results, which are obtained from a solution of the rate equations using the same equilibrium constants. Some conclusions are listed as follows;

1. No chemical reaction occurs for the present data at 2,000 K, which is reasonable from the results of the previous single direction chemical reaction in Section 6.2. The dissociation probability of oxygen molecules (O₂) is too low for lower temperature. However, some dissociation of oxygen molecules occurs in the rate equation solution at 2,000 K.

2. The simulated results of nitrogen molecule (N₂) are overestimated when the temperature is higher than 10,000 K. It may be because most of nitrogen molecules are dissociated under higher temperature and leads to too few particles to obtain correct results.

3. The oxygen molecules (O₂) are scattered because the particles are too rare when the temperature is higher. It needs more simulated particles and time-steps to predict the results.

4. The trends of nitric oxide (NO), nitrogen atom (N) and oxygen atom (O) are in very good agreement with the rate equation analysis.

5. If all the recombinations are neglected in the simulation, all the nitrogen molecules (N₂) and oxygen molecules (O₂) will be dissociated as nitrogen atoms (crossed circle) and oxygen atoms (crossed square) when temperature is

16,000K.

6.3 Concluding Remarks

This chapter has presented a chemical reaction module for the DSM C code. The total collision energy (TCE) model is used to calculate the reaction probabilities of dissociation and exchange. An extended TCE model, which is the three-body collision model, is applied to obtain the recombination probability. The function is successful verified by several benchmarks and most of those results agree with theoretical data. In general, the probabilities of dissociation increase with increasing temperatures and the recombination probabilities decrease with increasing temperatures from the single 2-D cell test cases. Those cases show the chemical reactions appear and need to be considered when the temperature is very high. There are some suggestions to the users of using chemical reaction module, which are listed in the following;

1. It needs a large total number of time-steps and a long transient time before sampling.

2. In general, the recombination reaction has lower probability to react. Thus, the recombination simulations need more time-steps to obtain enough sampling and accurate results.

Chapter 7 Applications and Examples

The proposed DSM C method combining mesh adaptation, variable time-step scheme, dynamic domain decomposition, conservative weighting scheme and chemical reaction has been verified successfully by the test problems in Chapters 3-6, respectively. Thus, to demonstrate the powerful capability of the current simulation tool, we will apply it to compute three rather challenging flows, which are an underexpanded twin-jet flow, a 3-D Apollo re-entry vehicle at 100 km altitude and a 3-D sphere shape re-entry vehicle at an altitude of 90 km. The results will then compare with previous simulated or experimental studies available in the literature. Note that either only preliminary results or problem definition will be described in the following.

7.1 Near-Continuum Underexpanded Twin-jet Interaction Flow and Simulation Conditions

For a truly three-dimensional flow, two parallel, near-continuum, under-expanded jets issuing from sonic orifices into a near-vacuum environment is selected as the first application since the experimental data is available [87]. Sketch of the twin-jet interaction is shown in Fig. 7.1(a). It is expected that a secondary jet form between the two primary parallel jets under certain flow conditions. Corresponding flow conditions represent a challenging problem since it involves flow regimes from near-continuum at the inlet to near free-molecular flow at the outlet, where the DSM C method may be the only available tool for analyzing this problem. Related flow conditions are listed as follows: the test gas is nitrogen gas; stagnation pressure P_o=870 Pa; stagnation temperature T_o=285 K; background pressure P_b=3.7 Pa; resulting pressure ratio P_o/P_b=235. Background pressure effect is either neglected or included in the simulation, in which previous simulation by Dagum and Zhu [88] has ignored and vacuum boundary is used instead for simplicity. The corresponding Knudsen number Kn_th (=λ_throat/D; D is the throat diameter, which equals 3 mm) is 0.00385.

Computational domain is reduced to 1/4 of the original physical domain by taking advantage of the geometrical symmetry of the flow. Height (H), width (W) and length (L) of the simulation domain are taken long enough as 10D, 10 D and 20 D, respectively, as shown in Fig. 7.1(b), where the positive z-direction is out of the paper. Not that the

distance between the centers of two primary jet centerline is 3 times the diameter of the orifice. Internal energy is relaxed through the Borgnakke-Larsen model [72], using a fixed rotational collision number of 5. The variable time-step method and an h-refinement mesh are used to reduce the computational time further and to increase the accuracy of the solution, respectively. An h-refinement three-dimensional mesh with mesh quality control is used in this simulation (32 processors on IBM -SP2 system at NCHC) to increase the accuracy of the solution. The resulting simulated particles are about 10 millions at steady state and the number of cells is approximately 0.66 million after two levels of mesh refinement (Fig. 7.2). Figure 7.2(a) shows the surface mesh along x-y and y-z planes, while Fig. 7.2(b) illustrates the exploded view of the surface mesh in the same planes near the sonic orifice, where the mesh is refined.

Corresponding local cell Knudsen number distribution is shown in Fig. 7.3, where most of the flow field satisfies the general cell-size requirement except the locations very near the sonic orifice (0.5<Kn_cc<1). Besides, 20,000 time-steps are used to sample the particles for obtaining macroscopic properties.

Dynamic Domain Decomposition

Figure 7.4 illustrates the initial and final domain decomposition on the surfaces of x-y, y-z and z-x planes for the parallel twin-jet interaction using vacuum outflow boundary. The initial domain decomposition (Fig. 7.4(a)) is obtained by assigning equal weight to each cell, which is the same as the three-dimensional sphere flow. Final domain decomposition changes dramatically as compared with the initial domain decomposition. In this case, the computational time is saved up to 100% using dynamic domain decomposition, as compared with that using static domain decomposition.

Similar trend is found for the case using pressure outflow boundary.

Properties Contour

Figure 7.5 is the particle distribution by using constant time-step and variable time-step schemes. As the figure illustrates, the particle distribution by using variable time-step scheme is less than those by constant variable time-step scheme. And the particle distribution is more uniform. Figure 7.6 presents the normalized density contours (with respect to the density at the sonic orifice) on x-y and y-z planes using vacuum outflow boundary. Note that y-z plane is the symmetric plane between the two primary jets. It is clear that a secondary jet is formed and centered along y-axis between the two primary jets. In addition, the density expends similar to a single jet on the side, where there is no jet interaction. Similar evidence for the formation of a secondary jet

between the two jets exists also for temperature (translational and rotational) contours in Fig. 7.7, which is label in degree of Kelvin.

Centerline Properties Distributions

Figure 7.8 illustrates the density and rotational temperature distribution along the symmetric centerline (y-axis), respectively, using different outflow boundary conditions (vacuum and pressure). Experimental data measured using laser induced fluorescence by Soga [87] and the DSM C simulation data by Dagum and Zhu [88] using vacuum outflow boundary condition are also included for comparison. It is clear that current simulation data using vacuum outflow boundary condition agree favorably with experimental data [87] and previous DSM C simulation data [88] using the same outflow boundary condition. It is, however, that the centerline density increases again for y/d>10 when pressure outflow boundary condition is used. Similar trend is also found for rotational temperature distribution in Fig. 7.8(b). Note that the background pressure in the vacuum chamber is maintained at 3.7 Pa, where no measurement location is clearly specified in the paper [87]. The only possible way to fully duplicate the experimental conditions is to conduct chamber-scale simulation, which is not possible due to the unclear experimental conditions provided in the paper [87].

7.2 Hypersonic Flow around a Three-Dimensional Apollo Re-entry Vehicle Flow and Simulation Conditions

The second application is a three dimensional Apollo case, which is at 100 km altitude and simulated with 34 chemical reactions [89]. Half of the whole Apollo geometry is simulated to save the computational time. Related flow conditions are listed as follows: the altitude of this case is 100 km and the angle of attack is 25^o; the initial number densities of inflow gases are determined by the M SIS-E-90 Atmosphere model (http://nssdc.gsfc.nasa.gov/space/model/models/msis.html), which are nitrogen, oxygen and oxygen atoms with 8.467E18, 2.025E18 and 3.995E17, respectively; the velocity and temperature are 9,592.64 m/s and 190.6 K, respectively; the wall temperature is fixed and is predicted by the Stephan-Boltzmann Law, which is 1,378 K; the corresponding Knudsen number Kn(based on the radius of aft compartment) is 0.033.

The cell number and particle number of this simulation are about 750,000 and 15,000,000, respectively. The mesh is shown as Fig. 7.9(a). Figure 7.9(b) is the particle locations on the slide near the symmetric plane. Each dot means a simulated particle.

There have more particles gathering in front of the object because the bow shock. And

fewer particles existed behind the object.

Properties Contour

Figure 7.10 is number density contour of each species, which is a slide of X-Z plane that can easily understand the flow field of Apollo case. Figure 7.10(a)~7.10(e) represent the number densities of N2, O2, NO, N and O, respectively. In these figures, a clear bow shock is formed in front of the object and the weak region is at under part behind the object. The chemical reactions are occurred at those regions to produce No and N because the temperature is heated by the bow shock. The temperature and velocity contours are shown as Fig. 7.11. Figure 7.11(a)~7.11(c) are the contours of translation, rotation and vibration, respectively. The maximum translational temperature is located on the bow shock near the upper tip of the object. The value reaches to 33,000 K which leads to chemical reactions occurred. These three temperature contours show strong degree of thermal non-equilibrium at the shock and wake regions. Figure 7.11(d)~7.11(f) illustrate the velocity contours in three coordinate directions near X-Z plane. In Fig. 7.11(d), the U-velocity along X-direction is decreased by the bow shock and then increased by rarefaction. The V-velocity along Y-direction in Fig. 7.11(e) is almost zero because it is a symmetric case. The W-velocity along Z-direction in Fig.

7.11(f) is interesting. Particles flowing through the lower part of the Apollo are faster than those through the upper part. And the separation exits behind the object.

7.3 A Three-Dimensional Re-entry S phere Flow Flow and Simulation Conditions

The third application is a three dimensional sphere case, which is at 90 km altitude and simulated with 19 chemical reactions [89]. According to the verifications of a single cell, the order of recombination probability is much less than the dissociation and exchange reactions. Thus, recombination is negligible in this simulation. Only one in sixteenth physical domain is simulated because its axis-symmetric geometry. Related flow conditions are listed as follows [69] the altitude of this case is 90 km and the free-stream density is 3.43E-6 kg/m³; the initial mole fraction of oxygen molecule, nitrogen molecule and oxygen atoms are 0.209, 0.788 and 0.004, respectively; the free-stream velocity and temperature are 7,500 m/s and 188 K, respectively; the wall temperature is 350 K with full diffuse surface; the corresponding Knudsen number Kn and Reynolds number Re(based on the diameter of the sphere, which is 1.6 meter) are 0.01 and 3,243, respectively. The two-level adaptive cell number and particle number of

this simulation are about 670,000 and 3,000,000, respectively. The adaptive mesh refinement can easily capture the bow shock in the front of the sphere because it is a hypersonic flow, which is shown as Fig. 7.12.

Properties Contour

Figure 7.13(a)~(c) illustrate flow field contours of this simulation, which including normalized density, normalized overall temperature and M ach number on the X-Z plane, respectively. The normalized properties are based on the free-stream condition. As we can see, there is a very strong bow shock stands in the front of the sphere and the wake region exists behind the object. The normalized density and overall temperature are increased along the X-axis first because the bow shock effect, which the maximum value are about 120 and 60, respectively. Figure 7.14 represents the mole fractions for five species, which are N2, O2, NO, N and O, along the stagnation line. Nitrogen and oxygen molecules stay constant before the bow shock and then decreasing near the shock region because the temperature goes very high, which can process the dissociation reaction. Thus, the mole fractions of nitric oxide, nitrogen atom and oxygen atom are increased. Figure 7.15 shows the surface coefficients, including pressure, skin-friction and heat transfer coefficients, along the symmetric line of the surface of the sphere. This figure is plotted as a function of the circumferential angle (θ) measured clockwise from the stagnation point and the simulation data of Dogra et al. [69] is also plotted into this figure for comparison. The pressure (Fig. 7.15(a)) and heat transfer (Fig.

7.15(c)) coefficient decreased with increasing circumferential angle and the maximum value are at stagnation point θ=0^o. The simulated result agrees very well with the reference [69]. The skin-friction coefficient increased then decreasing along the sphere surface and the maximum value is near θ=40^o. As we can see the results of skin-friction and heat transfer coefficients are pretty scattering because the particle sampling is not enough. M ore particle number and longer time-steps can overcome this problem and obtain more accurate simulated results. Although the results are not exact the same, the results are acceptable by comparing with reference data.

7.4 Concluding Remarks

This chapter presents several three-dimensional applications, which are a near-continuum parallel twin-jet, the Apollo re-entry vehicle flow and a hypersonic re-entry sphere flow. All these cases are interesting and important but it is very

complicated and the computation is time-consuming. By using the present PDSC, it can efficiently simulate and obtain accurate results by using reasonable computing time, which demonstrates its capability.

Chapter 8 Conclusions

8.1 Summary

The direct simulation of M onte Carlo (DSM C) method is used to simulate gas flows under rarefied gas environment in these four decades. M esh resolution, huge computational cost and flow involving specific problems are three main drawbacks by using the DSM C method. To overcome these disadvantages, a general-purpose parallel DSM C code (PDSC) is presented and verified by comparing with experiment and references. The current PDSC has following features;

1. It is a three-dimensional DSM C code with unstructured tetrahedral mesh, which is easier and flexible to deal with the flow with complex and irregular geometry. A cell-by-cell particle tracking technique can easily and correctly track particle movements.

2. A general unstructured adaptive mesh refinement with variable time-step scheme is developed to save computational time and to obtain accurate results.

The particle number of a 3-D hypersonic sphere flow by using constant time-step scheme is 1.7 million, which can be reduce to 0.34 million particles by applying variable time-step scheme and the transient time is decreased to only 25%. The computational efficiency can be speeded up to 10 times and more accurate result can be obtained if simulation uses suitable mesh resolution and time-step.

3. Parallelization of a 3-D DSM C method is developed to save the computing time. The parallel D SM C code with dynamic domain decomposition is also implemented to speed up the computational efficiency. The dynamic domain decomposition function can alleviate the unbalancing loading between each processor. This method can save 30-100% for the driven cavity flow by using static domain decomposition method.

4. Conservative weighting scheme (CWS) for flows with trace species is incorporated in PDSC to obtain reasonable number of simulated particles. A quasi 2-D hypersonic cylinder flow is used to verify this function. The total particle number by using constant weighting scheme is up to 10 million, which can be reduced to 0.24 million if the conservative weighting scheme is

activated.

5. Chemical reaction is developed for hypersonic reactive flows. It is tested by several cases, which includes comparisons the reaction probability (Pr) and rate constant (k_f) of each single reaction, degree of dissociation of pure species and the mole fraction of air (5 species) for a 2-D single cell. This function can correctly simulate dissociation, exchange and recombination reactions and the results agree with theoretical data.

Finally, several applications, which are a 3-D parallel underexpanded twin-jets, a hypersonic Apollo re-entry vehicle and a hypersonic re-entry sphere flow, are simulated to demonstrate potential and ability of the PDSC.

8.2 Recommended Future S tudies

There is still has something to make the PDSC complete. The diagram of recommended future studies is shown in Fig. 2.4 and described in the following:

1. To develop hybrid mesh code in PDSC. This can be a hybrid code with structured/unstructured and tetrahedral/hexahedral mesh system, which can save the cell number and reduce the time of tracking particle;

2. To couple with other numerical solves which can extend its capability of simulate complex flows. Combining with computational fluid dynamics (CFD) solver can solve the flow has continuum flow region. Combining with particle-in-cell (PIC) method can simulate flows with plasma. Incorporating with particle flux method (PFM ) can simulate inviscid flow.

References

1. Lee, Y. K. and Lee, J. W., “Direct simulation of compression characteristics for a simple drag pump model”, Vacuum, 47, pp. 807-809, 1996.

2. Lee, Y. K. and Lee, J. W., “Direct simulation of pumping characteristics for a model diffusion pump”, Vacuum, 47, pp. 297-306, 1996.

3. Chang, Y. W., et al., “Pumping performance analyses of a turbo booster pump by direct simulation M onte Carlo method”, Vacuum, 60, pp. 291-298, 2001.

4. Plimpton, S. and Bartel, T., “Parallel particle simulation of low-density fluid flows,” U.S. Department of Energy Report No. DE94-007858, 1993.

5. Yang, X. and Yang, J. M ., “M icromachined membrane particle filters”, Sensors and Actuators A, 73(1-2), pp. 184-191, 1999.

6. Piekos, E. S. and Breuer, K. S., “Numerical modeling of micromechanical devices using the direct simulation M onte Carlo method”, Transaction for ASM E Journal of Fluids Engineers, 118(3), pp. 464-469, 1996.

7. Nance, R. P., et al., “Role of boundary conditions in Monte Carlo simulation of microelectromechanical systems”, AIAA Journal, 12(3), pp. 447-449, 1998.

8. Tseng, K.-C., “Analysis of M icro-scale Gas Flows with Pressure Boundaries Using The Direct Simulation M onte Carlo M ethod”, M echanical Engineering, National Chiao-Tung University, Taiwan, M aster Thesis, 2000.

9. Schaff, S. and Chambre, P., Fundamentals of Gas Dynamics, Princeton University Press, Princeton, NJ, 1958, Chapter H.

10. Bird, G. A., M olecular Gas Dynamics, Clarendon Press, Oxford, UK, 1976.

10. Bird, G. A., M olecular Gas Dynamics, Clarendon Press, Oxford, UK, 1976.