- Email: [email protected]
Thermal protection of a hypersonic vehicle by modulating stagnation-point heat flux
Journal Pre-proof Thermal protection of a hypersonic vehicle by modulating stagnation-point heat flux
Quan Han, Chengdong Sun, Yi Tao, Zhongwu Li, Yan Zhang, Yunfei Chen
PII:
S1270-9638(19)33092-5
DOI:
https://doi.org/10.1016/j.ast.2019.105673
Reference:
AESCTE 105673
To appear in:
Aerospace Science and Technology
Received date:
29 November 2019
Revised date:
25 December 2019
Accepted date:
26 December 2019
Please cite this article as: Q. Han et al., Thermal protection of a hypersonic vehicle by modulating stagnation-point heat flux, Aerosp. Sci. Technol. (2020), 0, 105673, doi: https://doi.org/10.1016/j.ast.2019.105673.
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.
Thermal Protection of a Hypersonic Vehicle by Modulating Stagnation-Point Heat Flux Quan Han†, Chengdong Sun†, Yi Tao, Zhongwu Li, Yan Zhang*, Yunfei Chen* Jiangsu Key Laboratory for Design and Manufacture of Micro-Nano Biomedical Instruments, School of Mechanical Engineering, Southeast University, Nanjing 211189, China
Abstract: The heat flux at the stagnation-point is usually used to evaluate the heat transfer efficiency between high-temperature gas and a hypersonic flying vehicle. Controlling the stagnation-point heat flux is of critical important for the design of a thermal protection system for a hypersonic flying vehicle. In this work, we demonstrate that the nonequilibrium molecular dynamic (NEMD) method can well describe the shock wave and the stagnation-point heat flux. In particular, our results show that the heat flux is more sensitive to the change of the peak amplitude of the hypersonic flow in the bimodal velocity probability density function (PDF) at the stagnation point. Decreasing the peak amplitude of the hypersonic flow in the PDF can effectively reduce the stagnation-point heat flux. Three types of the outer shapes for the flying objects are employed to compare the heat flux at the stagnation-point. The NEMD simulation results confirmed that the heat flux can be decreased 16.8% for the spherical cavity outer shape compared with the flat shape in the leading edge. Keywords: stagnation-point heat flux; velocity probability density function; hypersonic flow; atom collisions
1.
Introduction A hypersonic flying vehicle has to face the challenge of aerodynamic heating as it flying in the
near atmosphere space. The heating comes from the strong shock wave formed on the leading edges of the vehicle and the viscous boundary layer attached to the whole vehicle body [1]. To avoid the damage of the vehicle, thermal protection systems (TPS) have to be deployed, including the insulated
瀗
澳Q. H. and C.S. contributed equally to the work澳 澳All correspondence should be addressed to Y. Zhang, Email: [email protected]澳or Y. Chen, Email: [email protected]; 澳 1
澽
structure, heat pipe, ablation, aerospike, cavity, jet and some combinations [2-9]. Experimental [10, 11] and numerical studies [12, 13] have proven that these TPS are effective if the flying time is short. However, for hypersonic vehicles with flying time a little bit longer, traditional TPS techniques become powerless due to the extremely high-temperature gas in front of the vehicle. In order to propose more efficient TPS techniques, the intrinsic mechanism of heat transfer between the high-temperature gases and the hypersonic flying vehicle must be figured out first. In a hypersonic flying vehicle, peak heat flux generally occurs at the stagnation point located in the region of nose and wing leading edges, which may induce enormous heat load on the vehicle. Thus, the theory of stagnation-point heat flux attracts lots of attention in the field of hypersonic thermal protection. In general, the flow state is categorized as continuous flow, transitional flow and free molecular flow by the non-dimensional parameter Knudsen number (Kn), defined as the ratio of the mean free path of gas molecules (ߣ) and the characteristic length (L). In 1950s, Lees [14], Fay and Riddell [15] proposed the stagnation-point heat flux formula by analytically solving the boundary equation simplified from the Navier-Stokes (NS) equations in the continuous flow regime (Kn < 0.01). It was estimated that the heat flux q at the stagnation point was inversely proportional to the square root of the nose radius R ( ן ݍඥሺͳȀܴሻ) [16-18]. However, in practical applications, the heat flux does not increase to infinity as predicted from the formula if the nose radius approaches zero. With the continual decrease of the nose radius, the Knudsen number will increase and be larger than 0.01. In this case, the simplified formula fails in predicting the heat flux. Furthermore, in the free molecular flow regime with Kn > 10, Regan and Anandakrishnan [19] introduced an analytical solution for the stagnation-point heat flux by solving the Boltzmann equation without the collision term, and the results indicated that the heat flux was independent of the nose radius. However, in the transitional flow regime as 0.01< Kn <10, it was much hard to obtain an analytical solution or the direct numerical solution for heat flux from the Boltzmann equation [20]. Subsequently, many simplified analytic methods developed from the Boltzmann equation were adopted, such as the linearized Boltzmann equation [21], the moment method [22] and the Bhatnagar-Gross-Krook (BGK) model [23]. In recent two decades, a new gas-kinetic scheme method developed from the BGK model for the entire Knudsen number flow, were constructed to solve the multi-scale problems from the continuous flow to the free molecular flow [24]. Although the theoretical studies of stagnation-point heat flux for continuous flow, 2
transitional flow, and free molecular flow have been applied to engineering applications, the solution for the heat flux is still not accurate enough. Particularly, at high Mach number, the governing equations including the Newton’s law for describing gas viscosity and the Fourier’s law for heat conduction may fail, which lead to the failure of the NS equations in describing the hypersonic gas flow [25]. Also, the linearized Boltzmann equation [21] is only suitable for the low speed flow problems, and the solutions for the moment method including the Burnett equation [26] and thirteen moment equations of Grad [27] are not exact at high Mach number. The widely used BGK model equation [23] is powerless in solving the problems that need to consider momentum and energy exchanges simultaneously, especially in predicting the heat flux and friction force of surface. With the growth of computer processing power, numerical methods are widely used to analyze gas flow at microscale or nanoscale level. In 1976, Bird [28] first implemented the direct simulation Monte Carlo (DSMC) method to delineate the circumstances under the condition in which the continuum model lost its validity. Up to now, the DSMC method has been considered as the most effective method for hypersonic rarefied flow [29-32]. However, the DSMC method heavily relies on the collision models, the molecular chaos hypothesis and dilute gas hypothesis [33], which makes the DSMC model fail in some cases. Besides the DSMC method, the molecular dynamic simulation (MD) is another popular method to study the microstructure of hypersonic shock wave [34], and has been widely used in studying the heat transfer mechanisms at microscale and nanoscale [35]. From the atomic or molecular scale, the MD method uses realistic empirical interatomic potential and numerically solves the Newton's second law of motion to further reduce the approximation. In 1978, Klimenko and Dremin [36] published the first non-equilibrium molecular dynamics (NEMD) simulation of shock waves in dense liquid argon, and the results of NEMD and the NS equations were in good agreement. During the past forty years, the developed NEMD methods had become an accurate tool to study the shock waves in the hypersonic flow [37]. With the help of the NEMD method, the velocity distribution can be effectively extracted and gradually become a powerful tool for studying shock waves. Although the NEMD method has well revealed the microscopic structure of the shock waves, most of the simulated models are based on the shock tube or the high-speed piston [25, 38, 39], while little researches focus on the model of a flying object in the gas environment. Intending to improve the understanding of heat transfer between the gas and a hypersonic flying vehicle at the atomic scale, this work focuses on describing the stagnation-point heat flux with the 3
NEMD method. First, a nanoscale model of a piece of silicon cuboid flying in an argon gas environment is established, and the density, temperature, heat flux and velocity distribution of the shock waves formed over the cuboid can be investigated. Subsequently, by changing the cross-section size of the silicon cuboid, the stagnation-point heat fluxes at different Kn numbers are studied through the atomic scale method and classical aerodynamic analysis methods. Especially, through comparing the bimodal velocity distributions at the stagnation point at different Kn numbers, the influence of the velocity distribution on the heat flux is studied. Next, simulations under different environmental atmospheric pressures and different geometry shapes are conducted to study the factors affecting the velocity distribution and the heat flux. Based on the NEMD model, the outer geometry of a hypersonic flying silicon cuboid is optimized to modulate the stagnation-point heat flux, which provides a technique route to reduce the heat load for a hypersonic vehicle. 2.
Methodology
2.1 Details of the NEMD simulation In this work, a silicon cuboid flying with a hypersonic velocity through an argon gas environment is simulated using the NAMD 2.12 package [40]. Neutral atoms or molecules can be accurately modeled with the 12-6 Lennard-Jones (LJ) potential [41].
σ
σ
VLJ (rij ) = 4ε (( )12 − ( )6 ) rij rij
(1)
where ɂ andɐare the depth of the potential well and distance at which the LJ interaction is zero, and rij is the distance between two particles. For argon atoms, ɂ ൌ ͲǤʹ͵ͺͳkcal/mole andɐ ൌ ͵ǤͶͳͳͺՀ [42]. CHARMM36 force fields [43] are used to describe the interactions between silicon atoms, and the nonbonded interaction for silicon atoms is described by the LJ potential with the parameters ɂ ൌ ͲǤ͵ͳͲͲ kcal/mole and ɐ ൌ ͵ǤͺͲͷՀ [44]. Moreover, the nonbonded interactions between two different particles are calculated using the mixing rules [45]. In our simulations, the cutoff radius is set as 8Հ.
4
Figure 1. Schematic diagram of a silicon cuboid flying in an argon gas environment along the z-direction.
First, all simulations for studying heat flux at different Kn numbers, pressures, as well as the geometry shapes, are performed on the same system. The schematic diagram of the system is displayed in Fig. 1, in which the silicon cuboid flies along the z-direction surrounded by the gaseous argon. Periodic boundary conditions are applied along the three axis directions. For the cases at different Kn numbers, the initial states of the argon gas environment are set as p0 = 10 atm and T0 = 300 K, and the mean free path of the gas is predicted to be 7 nm [46]. For silicon, its lattice constant is fixed at the optimized value of 5.43Հ. Table 1 lists the details of each simulation when the cross-section size of the silicon cuboid varies from 4 to 32 unit cells (UC), and the corresponding global Kn number varies from 3.2 to 0.4 at the tenfold atmospheric pressure. For the cases at different pressures and geometry shapes, the details are described in section 3.2 and section 3.3. Table 1. Size of the silicon cuboid, Kn number, box size, and the number of atoms Case
Pressure
Silicon
Global
Box
Silicon
Argon
No.
(atm)
Size (UC3)
Kn number
Size (Հଷ )
Number
Number
1
10
ͶൈͶൈͶ
3.2
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
512
25969
2
10
ൈൈͶ
2.1
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
1152
25964
3
10
ͺൈͺൈͶ
1.6
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
2048
25958
4
10
ͳʹ ൈ ͳʹ ൈ Ͷ
1.1
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
4608
25940
5
10
ͳ ൈ ͳ ൈ Ͷ
0.8
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
8192
25915
6
10
͵ʹ ൈ ͵ʹ ൈ Ͷ
0.4
ͺͲͲ ൈ ͺͲͲ ൈ ͳͲͲͲ
32768
168197
For each case, the velocity Verlet algorithm with a time step of 1 fs is used for the integration of the motion equations [47]. The system is first minimized for 0.1ns and then equilibrated in an isothermal-isobaric (NPT) ensemble for 2 ns. Then in the canonical (NVT) ensemble, the silicon cuboid is pulled along the positive z-direction with a speed of 2000 m/s by a driving spring, named the 5
Steered Molecular Dynamics method [48]. In order to keep the flying direction along z-direction, the silicon cuboid atoms are harmonically restrained in x and y directions. The argon atoms and silicon atoms are coupled with Langevin thermostats at 300 K with different damping coefficients. The simulation is carried out for at least 7 ns to guarantee the system reaching a steady state, and then the heat flux is collected to certify the convergence of the NEMD calculations. In this work, a series of simulations have been performed to optimize the parameters in constructing the simulation systems (Figure. S1, Supplementary material). Macroscopic observables, including density, temperature, and heat flux are calculated as the spatial and time average in the volume of a prescribed grid cell, with a size defined asȟ ݔൈ ȟ ݕൈ ȟݖ. When plotting the curves of macroscopic parameters along the z-direction, we set the length of ȟ ݖas 7Հ and the size of ȟ ݔൈ ȟ ݕas the cross-section size of the silicon cuboid in the x-y plane (consistent with the grid definition of the stagnation-point region, Section S2, Supplementary material). When drawing the contour of density in the x-z plane, the size of ȟ ݔൈ ȟ ݖis set as ͳ ൈ ͳՀଶ and the length of ȟ ݕis set as the length of the silicon cuboid in the y-direction. In the MD simulations, the local heat flux can be calculated as the time derivative of the sum of the site energy moments over the atoms in each grid cell [49].
G q=
1 d N JG ¦ ri Ei ΔxΔyΔz dt i =1
(2)
N 2 1 JG Ei = m vi − v + ¦Vij 2 j =1
(
)
(3)
where ݒҧ is the averaged velocity of all the atoms in each grid cell, ሬሬሬԦ ݒప is the velocity of each atom i, and Vij is the potential energy of atom i due to atom j. The formula of the local heat flux is finally derived as: G q=
N JJG JG º ª JG 1 1 N JJG ª JG − + E v v f ij < vi − v + v j − v º
(
( (
)
) (
) )
(4)
డೕ where ሬሬሬԦ ݂ ൌ ሬሬሬሬሬԦ is the force on atom i due to atom j. డഢണ
2.2 Model verification To verify that our simulation setup can effectively depict the essential physics in a hypersonic flow, 6
the number density, temperature, and heat flux profile of the gas flow in the steady state are extracted from the NEMD model and displayed in Fig. 2. In this verification case, a silicon cuboid with size of
ͳ ൈ ͳ ൈ Ͷଷ flying in the gaseous argon with a speed v0 of 2000 m/s (about 6.2 Mach number) is simulated, which is listed as case No. 5 in Table. 1. As expected, a shock wave is generated in front of the flying object demonstrated as the normalized number density contour in the inset of Fig. 2(a). In the shock wave, an abrupt change in normalized number density of the gaseous argon can be observed from Fig. 2(a). The normalized number density reaches the peak in position (Ϩ) near the silicon cuboid front surface. With the increase of the distance away from the silicon cuboid front surface, the normalized number density of the gaseous argon starts to decrease. At position (IV), the number density of the gaseous argon decreases to the environmental argon density. This profile of the normalized number density of the gaseous argon matches qualitatively well with the numerical solutions of the augmented Burnett equation [50, 51] and the DSMC results [52]. Besides, the temperature curve of the gaseous argon is provided in Fig. 2(b), and a peak value occurs in the position (ϩ) of the shock wave. From the heat flux computed by Eq. (4), it can be seen that a negative heat flux locates at the left of the position (ϩ) in Fig. 2(b), meaning that heat is transferred from the shock wave to the silicon cuboid, while the positive heat flux indicates that heat is transferred from the shock wave to the gas far away from the silicon cuboid. Both the temperature gradient and heat flux direction indicate the rationality of the simulation model.
7
Figure 2. The NEMD simulation results. (a) Normalized number density n/n0, where n denotes the number density inside each grid and n0 is the initial gaseous number density. (b) Temperature and heat flux calculated in NEMD are given as functions of distance in the shock wave. The normalized number density contour shows a shock wave generated before the hypersonic flying silicon cuboid with a speed of v0 = 2000 m/s. The grey region stands for the silicon cuboid. (c) The probability density function of argon atom velocity in the z-direction in four different positions: near the front surface (Ϩ), shock wave (ϩ) and (Ϫ), and far away from the front surface (ϫ). The black dots stand for the results extracted from the NEMD simulation and the red curves are predicted from the Mott-Smith theory.
Notably, the velocity distribution of the argon atoms in the shock wave does not obey the Maxwell distribution due to the high nonequilibrium flow. The Mott-Smith theory [22] describes a shock wave as a superposition of a hypersonic flow and a subsonic flow, and their proportions vary with the position along the z-direction, which is named as the bimodal velocity distribution function:
f = fs + f h
(5)
G JJG m(v − vα ) 2 m 3/2 fα = nα ( z )( ) exp( − ) 2π k BTα 2k BTα
(6)
where
where f is the velocity distribution function and the subscript s and h denote the subsonic flow and 8
hypersonic flow, respectively. The ߙ represents the subscript s or h, and the݊ఈ ሺݖሻdenotes the number density varied with the position along the z-direction. The gaseous temperature ܶఈ and velocity of gas flowݒ ሬሬሬሬԦ ఈ are assumed to be independent of z, and they are calculated through the Rankine-Hugoniot relations [1]. Fig. 2(c) displays the velocity distribution of argon atoms in the z-direction from the NEMD simulation and the Mott-Smith theory in the four positions (Ϩ), (ϩ), (Ϫ) and (ϫ) indicated in Fig. 2(a). In order to describe the velocity PDF more obviously, we consider the flying silicon cuboid as the frame of reference, in which the velocity for the argon atoms move toward the cuboid is defined to be positive, and the velocity is normalized by the speed v0 of silicon cuboid. In NEMD, the velocity distribution is given as a probability density function (PDF), defined as fv = ni /(n¨v), and where ni denotes the number density in the velocity range from vi/v0 - ¨v/2 to vi/v0 + ¨v/2 and n denotes the number density inside each grid. In the Mott-Smith theory, the velocity PDF is obtained by the velocity distribution function f in the z-direction, which is properly normalized so that ݂ ௩ ݀ݒ௭ over all velocities equals one. First, a nearly undisturbed hypersonic flow locates in the position (ϫ), and its velocity PDF is a single sharp peak centered at the velocity v0 of 2000 m/s. As the flow approaches the silicon cuboid, the argon atoms in position (Ϫ) experience a few collisions with the atoms reflected back by the cuboid obstacle and a subsonic flow appears in the distribution. As expected, more collisions with the reflected atoms occur in position (ϩ), resulting in a higher peak amplitude of subsonic flow. When it comes to the position (Ϩ) near the silicon cuboid front surface, the initial hypersonic flow is almost replaced by the subsonic flow due to the numerous atom collisions. These observations show that when the initial hypersonic flow gradually approaches the object, the flow in shock wave gradually becomes a superposition of the hypersonic flow and subsonic flow, which is called bimodal velocity distribution. However, there is an obvious quantitative difference in the subsonic flow between the NEMD results and the Mott-Smith theory in Fig. 2(c). In the Mott-Smith theory, the temperature and velocity of subsonic flow behind a shock wave are calculated through the Rankine-Hugoniot relations, which are based on the three conservation equations of mass, momentum, energy and the equation of state for ideal gases. Actually, the Rankine-Hugoniot relations fail in accurately depicting the parameters behind a shock wave in this case due to the highly nonequilibrium flow. But it is worth noting that the 9
overall tendency of velocity PDF obtained from NEMD simulation is quite reasonable in Fig. 2(c), further validating the feasibility of the NEMD method in simulating the hypersonic flow.
3.
Results and discussion
3.1 The effect of Kn number on the stagnation-point heat flux In our data analysis process, the grid length is a critical parameter to calculate the macroscopic parameters such as density, temperature and heat flux. To determine the optimal lengthȟݖnear the front surface, we perform a series of statistical analysis with different grid lengths ranging from 3 to 11Հ (Figure. S2, Supplementary material). Length of 7 Հ is chosen for further study, because it can ensure a sufficient number of atoms in a statistical grid cell and is small enough to guarantee the accuracy of the definition that considering position (Ϩ) as the stagnation point in the immediate neighborhood of the silicon cuboid front surface, besides, it produces a convergent heat flux. To investigate the effect of Kn number on the heat flux in the stagnation-point region (i.e., stagnationpoint heat flux), we simulate six cases of different Kn numbers with an atmospheric pressure of 10 atm. In our model, the different Kn number ( ݊ܭൌ ߣȀ )ܮis obtained through changing the cross-section size L of the silicon cuboid as shown in Table. 1.
3.1.1 Details in the stagnation-point region First, we study the influence of the Knudsen number Kn on the physical parameters of gaseous argon in the stagnation-point region, including the number density and velocity distribution. During the flight of the silicon cuboid, more atom collisions between the free-stream flow and the silicon cuboid front surface occur when the cross-section size becomes larger. Thus, more argon atoms are reflected back to the face coming flow, resulting in a larger number density ns at the stagnation point in the steady state. With the Kn number increasing from 0.4 to 3.2, the number density decreases at the stagnation point as listed in Fig. 3. For example, as the Knudsen number Kn increases to 3.2, the number density at the stagnation point decreases to 0.76/nm3. It is believed that the events of the atom collision become less in front of the flying cuboid silicon at the high Kn number, which leads to the decrease of the number density. Similarly, the phenomenon that the events of the atom collision become less also significantly affects the velocity PDF, in which the peak amplitude of the hypersonic flow gradually increases at the larger Knudsen number, and the peak amplitude of subsonic flow 10
decreases correspondingly. It can be predicted that when the Kn number increases infinitely, that is, the cross-section size of silicon cuboid decreases to zero, the velocity PDF will be the same as that of the free-stream flow with a single hypersonic peak.
Figure 3. The velocity PDF of argon atoms in the z-direction at the stagnation point (Ϩ) as the Kn numbers vary from 0.4 to 3.2. The number density ns at the stagnation point are listed in the upper right corner of each figure. The dots represent the data extracted from the NEMD simulations and the curves are the fitted bimodal velocity PDF.
3.1.2 Heat flux at the stagnation point Several models are used to calculate the heat flux at the stagnation-point. If the velocity distribution function f in Eq. (5) is known, various macroscopic parameters including the heat flux at the stagnation point can be accurately obtained by calculating the moment of f. In the method of the “moment of f ”, the flux of translational energy Ei brought to the silicon front surface per area in unit time by all the incident argon atoms is first calculated by the velocity distribution function. Then the flux of translational energy Er carried away by those reflected atoms is computed by the velocity distribution function for the reflected atoms, which is determined by the law of interaction of the argon atoms with the surface. Finally, the heat flux to the silicon cuboid front surface can be expressed as q = Ei – Er [53]. Besides the model of the moment equation, the method of bridging function is also 11
employed to estimate the heat flux from continuous flow to free molecular flow. By using the classical aerodynamic solutions of the stagnation-point heat flux for continuous flow and free molecular flow, the bridging function can be obtained at the entire Kn number by fitting the free parameters a1, a2 and a3 to the NEMD results. The bridging function is written as a function of Kn number [54] :
qb = qFM × Pb + (1 − Pb ) × (qc / (1 − a3 × Kn))
(7)
Pb = sin(ϕ )
(8)
ϕ = π (a1 + a2 × log10 ( Kn))
(9)
where a1, a2, a3 are the fitting parameters. qFM is the analytical solution of the heat flux for free molecular flow [55]. qc is the heat flux calculated from the Fay-Riddell model solving the NS equations in the continuum flow [15], and the effective radius used in the Fay-Riddell model for the silicon cuboid can be computed from ref. [56].
Figure 4. Heat flux coefficients defined as h = q/qFM are given as a function of the Kn number obtained from several theoretical models and the NEMD. The NEMD data is computed from Eq. (4). The model of “moment of f ” denotes that the heat flux is calculated from the fitted velocity distribution f , and the proportions of the heat fluxes contributed by hypersonic flow and subsonic flow at different Kn numbers are presented in the inset. The bridging function is shown as Eq. (7) and the Fay-Riddell model denotes that the stagnation-point heat flux is solved from the NS equations. Free molecular flow model means that the heat flux is calculated from the Boltzmann equation. (Argon gas, T0 = 300K; Silicon cuboid, T0 = 300K, v0 = 2000 m/s, the sizes are shown in Table. 1).
When studying the influence of the Kn number on the stagnation-point heat flux, the heat flux is generally normalized by the value for the free molecular flow, called the heat flux coefficient h. In Fig. 12
4, we compare the coefficient h obtained from the NEMD with the calculations by other methods, in which the green dashed line with a value of one denotes the results for free molecular flow. Due to the limited computing resource, the minimum Kn number of the NEMD models is 0.4. As can be seen in Fig. 4, the heat flux coefficient predicted by all the methods increases with the Kn number, which is consistent with the previous report that the smaller cross-section size with larger Kn number contributes to the larger heat flux [57]. However, the result of the Fay-Riddell model is significantly different from the NEMD results, which indicates the failure of the linear law in computing the heat flux when it comes to the transitional flow. On the contrary, the heat flux coefficient from the bridging function agrees well with the NEMD results in Fig. 4, with the fitting parameters in the expression as a1 = 0.12, a2 = 0.32, a3 = 0.8 in Eq. (7). This consistency shows that the NEMD results can be well depicted by the bridging function. It is worth noting that the results of the “moment of f ” are almost identical to the NEMD results in Fig. 4, indicating that the heat flux can be accurately derived from the velocity distribution function. Actually, this velocity distribution function f is obtained by the fitted velocity PDF fv multiplied by the number density ns in Fig. 3. To further understand the roles of hypersonic and subsonic flow in the stagnation-point heat flux, we study the proportions of heat flux contributed by hypersonic and subsonic flow with the model of “moment of f ”. As shown in the inset of Fig. 4, the heat flux contributed by hypersonic flow increases significantly with the Kn number, and the heat flux contributed by subsonic flow declines rapidly. And in the distribution function f , the peak amplitude of the hypersonic flow still become higher although the number density ns decreases significantly at the larger Kn number (Figure S3, Supplementary material). Therefore, an increase in the stagnationpoint heat flux is observed due to the larger increment of heat flux contributed by the hypersonic flow. In addition, although the peak of the hypersonic flow in Fig. 3 (a) is extremely tiny, the proportion of heat flux contributed by hypersonic flow still reaches to 37% when Kn number equals 0.4. When the Kn number rises to 3.2, the proportion of heat flux contributed by hypersonic flow is up to 91% while the peak of the hypersonic flow in the velocity PDF is still lower than the peak of subsonic flow. These results show that the stagnation-point heat flux is more sensitive to the change of the peak amplitude of the hypersonic flow. Based on this finding, a design criterion is proposed, i.e. an efficient way to reduce the heat flux can be realized by decreasing the peak amplitude of the hypersonic flow in the velocity PDF. 13
3.2 Effect of pressure on velocity distribution and heat flux at the stagnation point In order to investigate the factors affecting the velocity distribution, we carry out the NEMD simulations of silicon cuboid with the size of ͳ ൈ ͳ ൈ Ͷ ଷ under different environmental atmospheric pressures. As can be observed in Fig. 5(a), the normalized number density n/n0 decreases rapidly to 1.0 along the flying direction. Specifically, the normalized number density at the stagnation point decreases rapidly with the increase of the environmental atmospheric pressure, while only a small difference can be seen when z is larger than 20 Հ . Moreover, in the inset of Fig. 5(a), the peak amplitudes (v/v0=1) in the velocity PDF are almost the same under different atmospheric pressures. It is believed that, during the flight of silicon cuboid under different atmospheric pressures, the number of argon atoms reflected back after colliding with the front surface increases proportionally when the number density n0 of the free-stream flow increases with the pressure, which leads to the same influence on the face coming gas flow. Therefore, we infer that the peak amplitude of the hypersonic flow in the velocity PDF is independent of atmospheric pressure but dependent on the cross-section size of the silicon cuboid.
Figure 5. (a) Normalized number density is given as a function of distance in the z-direction. The size of silicon cuboid is ͳ ൈ ͳ ൈ Ͷ ଷ and the initial atmospheric pressures are 1 atm, 5 atm, 10 atm, and 20 atm respectively. The yellow region stands for the stagnation-point region. The number density and velocity PDF at the stagnation point are shown in the inset. The stagnation-point heat flux as a function of different atmospheric pressure is presented in (b). The inset shows the heat flux q from NEMD simulation versus the number density ns at the stagnation point.
As can be seen in Fig. 5(b), the magnitude of the stagnation-point heat fluxes increase with the 14
pressure. Especially, compared with the linearly increasing heat flux for free molecular flow, the NEMD result increases slowly with the environmental atmospheric pressures. However, at low atmospheric pressure, the gap between the heat flux predicted from the free molecular flow and that extracted from the NEMD becomes small because the flow at the low pressure can be approximated as free molecular flow. We also find that the Fay-Riddell model overestimate the heat flux compared with the NEMD results, and this deviation shows the failure of NS equations in capturing the heat flux in the transitional flow regime. Furthermore, the fitted bridging function still agrees well with the NEMD results at these atmospheric pressures. Besides, it is notable that the heat flux is almost linearly proportional to the number density at the stagnation-point as shown in the inset of Fig. 5(b).
3.3 Modulating stagnation-point heat flux by the outer shape of the leading edge In the previous sections, we find that an efficient method for reducing the stagnation-point heat flux is through decreasing the peak amplitude of the hypersonic flow, which is closely related to the atom collisions between the free-stream flow and the silicon cuboid front surface. Since recent researches have shown that the stagnation-point heat flux can be significantly reduced by introducing an axial cavity in the nose region of a hypersonic vehicle [58], we implement three models of the silicon cuboids with different outer shapes on the front head to verify whether the heat flux is proportional to the peak amplitude of the hypersonic flow. The shapes of silicon cuboid are shown in the inset of Fig. 6(a) as S1, a silicon cuboid with a sphere cavity on the front surface, S2, a silicon cuboid with a flat surface, and S3, a silicon cuboid with a sphere front surface, in which the depth for both the sphere cavity and sphere is set as Z* equaled to 2 UC. The statistical grid is shown as the region with red dashed lines in the inset of Fig. 6(a), which is in front of the silicon surface within a length ȟ ݖequaled to 7Հ.
15
Figure 6. (a) Normalized number density of gas atoms is given as a function of distance in the z-direction. The atmospheric pressure is set at 10 atm. The yellow region stands for the stagnation-point region. The flying object is a silicon cuboid with cross-section size of ͳ ൈ ͳ ଶ and having different outer shapes in the head. In the inset of (a), the magnitude of heat flux at the stagnation point from the NEMD simulations are presented. The velocity PDF are presented in (b), and the inset shows the velocity distribution function f of the hypersonic flow.
First, the normalized number density n/n0 for the flying silicon cuboid with three outer shapes are shown in Fig. 6(a). The normalized number density at the shape S1 of a sphere cavity is the largest, while the value at the shape S3 of a sphere is the smallest. In fact, during the flight of the silicon cuboid, the atom collisions between the free-stream flow and the front surface at both cases of S1 and S3 are stronger than that at the case S2 due to the larger surface area. However, more argon atoms slip away the object after colliding with the spherical surface at the case S3, so fewer argon atoms are reflected back to the shock wave layer, resulting in the smallest normalized number density. On the contrary, more argon atoms can be reflected back to form the shock wave layer at the case S1 because of the cavity, so the normalized number density is the largest. Similarly, the peak amplitude of the hypersonic flow at the shape S1 is the lowest because large amount of argon atoms are reflected back, and the peak amplitude of the hypersonic flow in the velocity PDF increases obviously from the shape S1 to S3 as shown in Fig. 6(b). In particular, despite the decreasing of stagnation-point number density from the shape S1 to the shape S3, the peak amplitude of the hypersonic flow still increases in the velocity distribution function f as shown in the inset of Fig. 6(b). Thus, the magnitude of the stagnation-point heat flux increases from the shape S1 to the shape 16
S3 as shown in the inset of Fig. 6(a). And compared to the shape of the flat surface, the stagnationpoint heat flux is reduced by 16.8% due to the decreased peak amplitude of the hypersonic flow at the shape of the sphere cavity, indicating that the cavity shape can be used to reduce the heat flux.
4.
Conclusion In summary, the NEMD simulations are carried out to study the stagnation-point heat flux of the
hypersonic vehicle at nanoscale. Through analyzing the gaseous number density, temperature, heat flux and velocity probability density function along the flight direction, it is validated that the hypersonic flow can be well described by the NEMD simulations. Through the NEMD simulations, it is found that the heat flux increases with the Kn number, which is attributed to the peak amplitude rise of the hypersonic flow in the velocity PDF. Furthermore, when the peak amplitude of the hypersonic flow is the same while the peak amplitude of subsonic flow is different in the velocity PDF, the stagnation-point heat flux increases almost linearly with the stagnation-point number density. Based on the results, a conclusion is derived that the stagnation-point heat flux is more sensitive to the change of the peak amplitude of the hypersonic flow in the velocity PDF due to the much large translational energy stored in the hypersonic flow. Thus, a design criterion is proposed to reduce the stagnation-point heat flux, i.e. decreasing the peak amplitude of the hypersonic flow in the velocity PDF. Three type outer shapes for the leading edge of a hypersonic flying silicon cuboid are employed to verify the design criterion. In addition, we investigate the factors that affect the peak amplitude of the hypersonic flow in the velocity PDF. During the hypersonic flight of the silicon cuboid, argon atoms are reflected back by the leading edge and collide with the face coming flow, resulting in the formation of the shock wave. The results indicate that the peak amplitude of the hypersonic flow in the velocity PDF is independent of the atmospheric pressures but dependent on the geometry size of silicon cuboid. Our study provides a method to decrease the peak amplitude of hypersonic flow to reduce heat flux, which will be helpful to innovate the TPS design at nanoscale in a hypersonic flying vehicle.
Acknowledgement The authors thank the financial support from the National Key R&D Program of China (2018YFB1105400) and the Natural Science Foundation of China (Grants No.51435003, 51575104, 17
51705074).
References [1] A. Viviani, G. Pezzella, Aerodynamic and aerothermodynamic analysis of space mission vehicles, Springer press2015. [2] W. Huang, L. Yan, J. Liu, L. Jin, J.-g. Tan, Drag and heat reduction mechanism in the combinational opposing jet and acoustic cavity concept for hypersonic vehicles, Aerospace Science and Technology, 42 (2015) 407-414. [3] P. Mashaei, M. Shahryari, S. Madani, Analytical study of multiple evaporator heat pipe with nanofluid; A smart material for satellite equipment cooling application, Aerospace Science and Technology, 59 (2016) 112-121. [4] Q. Qin, J. Xu, Numerical evaluation of aerodome and cooling jet for aeroheating reduction, Aerospace Science and Technology, 86 (2019) 520-533. [5] M. Rivier, J. Lachaud, P.M. Congedo, Ablative thermal protection system under uncertainties including pyrolysis gas composition, Aerospace Science and Technology, 84 (2019) 1059-1069. [6] L. Zhu, Y. Li, X. Chen, H. Li, W. Li, C. Li, Hypersonic flow characteristics and relevant structure thermal response induced by the novel combined spike-aerodome and lateral jet strategy, Aerospace Science and Technology, 95 (2019) 105459. [7] P.Q. Elias, N. Severac, J.M. Luyssen, Y.B. Andre, I. Doudet, B. Wattellier, J.P. Tobeli, S. Albert, B. Mahieu, R. Bur, A. Mysyrowicz, A. Houard, Improving supersonic flights with femtosecond laser filamentation, Science Advances, 4 (2018). [8] S. Li, W. Huang, J. Lei, Z. Wang, Drag and heat reduction mechanism of the porous opposing jet for variable blunt hypersonic vehicles, International Journal of Heat and Mass Transfer, 126 (2018) 1087-1098. [9] R.-r. Zhang, W. Huang, L.-q. Li, L. Yan, R. Moradi, Drag and heat flux reduction induced by the pulsed counterflowing jet with different periods on a blunt body in supersonic flows, International Journal of Heat and Mass Transfer, 127 (2018) 503-512. [10] Z.-g. Wang, X.-w. Sun, W. Huang, S.-b. Li, L. Yan, Experimental investigation on drag and heat flux reduction in supersonic/hypersonic flows: a survey, Acta Astronautica, 129 (2016) 95-110. [11] M. Ou, L. Yan, W. Huang, S.-b. Li, L.-q. Li, Detailed parametric investigations on drag and heat flux reduction induced by a combinational spike and opposing jet concept in hypersonic flows, International Journal of Heat and Mass Transfer, 126 (2018) 10-31. [12] X. Sun, W. Huang, M. Ou, R. Zhang, S. Li, A survey on numerical simulations of drag and heat reduction mechanism in supersonic/hypersonic flows, Chinese Journal of Aeronautics, (2019). [13] W. Huang, Z. Chen, L. Yan, B.-b. Yan, Z.-b. Du, Drag and heat flux reduction mechanism induced by the spike and its combinations in supersonic flows: A review, Progress in Aerospace Sciences, (2018). [14] L. Lees, LAMINAR HEAT TRANSFER OVER BLUNT-NOSED BODIES AT HYPERSONIC FLIGHT SPEEDS, Jet Propulsion, 26 (1956) 259-&. [15] J.A. Fay, F.R. Riddell, THEORY OF STAGNATION POINT HEAT TRANSFER IN DISSOCIATED AIR, Journal of the Aeronautical Sciences, 25 (1958) 73-&. [16] N.H. Kemp, Laminar heat transfer around blunt bodies in dissociated air, Journal of the aerospace sciences, 26 (1959) 421-430. [17] R.B. Pope, Stagnation-point convective heat transfer in frozen boundary layers, AiAA journal, 6 (1968) 619-626. [18] K. Sutton, R.A. Graves Jr, A general stagnation-point convective heating equation for arbitrary gas mixtures, 1971. [19] F.J. Regan, S.M. Anandakrishnan, Dynamics of atmospheric re-entry, American Institute of Aeronautics and Astronautics press1993. [20] N. Singh, T.E. Schwartzentruber, Heat flux correlation for high-speed flow in the transitional regime, Journal of Fluid Mechanics, 792 (2016) 981-996. [21] E.P. Gross, E.A. Jackson, Kinetic models and the linearized Boltzmann equation, The physics of fluids, 2 (1959) 43218
441. [22] H.M. Mott-Smith, The solution of the Boltzmann equation for a shock wave, Physical Review, 82 (1951) 885. [23] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Physical review, 94 (1954) 511. [24] K. Xu, A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method, Journal of Computational Physics, 171 (2001) 289-335. [25] E. Salomons, M. Mareschal, Usefulness of the Burnett description of strong shock waves, Physical review letters, 69 (1992) 269. [26] D. Burnett, The distribution of molecular velocities and the mean motion in a nonϋuniform gas, Proceedings of the London Mathematical Society, 2 (1936) 382-435. [27] H. Grad, On the kinetic theory of rarefied gases, Communications on pure and applied mathematics, 2 (1949) 331407. [28] G.A. Bird, Molecular gas dynamics, NASA STI/Recon Technical Report A, 76 (1976). [29] G. Dahlen, M. Macrossan, C. Brundin, J. Harvey, Blunt cones in rarefied hypersonic flow: Experiments and MonteCarlo simulations, (1984). [30] G. Pham-Van-Diep, D. Erwin, E. Muntz, Nonequilibrium molecular motion in a hypersonic shock wave, Science, 245 (1989) 624-626. [31] P. Shufflebotham, T. Bartel, B. Berney, Experimental validation of a direct simulation by Monte Carlo molecular gas flow model, Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures Processing, Measurement, and Phenomena, 13 (1995) 1862-1866. [32] J. All-ccedil, gre, D. Bisch, J. Lengrand, Experimental rarefied density flowfields at hypersonic conditions over 70degree blunted cone, J Spacecraft Rockets, 34 (1997) 714-718. [33] S.M. Deshpande, P.S. Raju, Monte carlo simulation for molecular gas dynamics, Sadhana, 12 (1988) 105-123. [34] B. Holian, Atomistic computer simulations of shock waves, Shock waves, 5 (1995) 149-157. [35] T. Luo, G. Chen, Nanoscale heat transfer–from computation to experiment, Physical Chemistry Chemical Physics, 15 (2013) 3389-3412. [36] V. Klimenko, A. Dremin, Detonatsiya, Chernogolovka, Akademik Nauk Moscow, 1978. [37] P. Valentini, P.A. Tump, C. Zhang, T.E. Schwartzentruber, Molecular dynamics simulations of shock waves in mixtures of noble gases, J Thermophys Heat Tr, 27 (2013) 226-234. [38] F. Uribe, R. Velasco, L. Garcia-Colin, Burnett description of strong shock waves, Physical review letters, 81 (1998) 2044. [39] W.G. Hoover, C.G. Hoover, Shockwaves and Local Hydrodynamics: Failure of the Navier-Stokes Equations,
New
Trends In Statistical Physics: Festschrift in Honor of Leopoldo García-Colín's 80th Birthday, World Scientific2010, pp. 1526. [40] J.C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R.D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, Journal of Computational Chemistry, 26 (2005) 1781-1802. [41] J.E. Jones, On the determination of molecular fields.—II. From the equation of state of a gas, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 106 (1924) 463-477. [42] M.P. Allen, D.J. Tildesley, Computer simulation of liquids, Oxford university press2017. [43] A.D. MacKerell Jr, D. Bashford, M. Bellott, R.L. Dunbrack Jr, J.D. Evanseck, M.J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, All-atom empirical potential for molecular modeling and dynamics studies of proteins, The journal of physical chemistry B, 102 (1998) 3586-3616. [44] J.A. Wendel, W.A.G. Iii, The Hessian biased force field for silicon nitride ceramics: Predictions of thermodynamic and mechanical properties for ¢and £· i3N4, Journal of Chemical Physics, 97 (1992) 5048-5062. [45] Z. Xu, M.J. Buehler, Nanoengineering heat transfer performance at carbon nanotube interfaces, ACS nano, 3 (2009) 19
2767-2775. [46] G. Chen, Nanoscale energy transport and conversion: a parallel treatment of electrons, molecules, phonons, and photons, Oxford University Press2005. [47] D. Frenkel, B. Smit, Understanding molecular simulation: From algorithms to applications, 1996. [48] B. Isralewitz, J. Baudry, J. Gullingsrud, D. Kosztin, K. Schulten, Steered molecular dynamics investigations of protein function, Journal of Molecular Graphics & Modelling, 19 (2001) 13-25. [49] Z. Fan, L.F.C. Pereira, H.Q. Wang, J.C. Zheng, D. Donadio, A. Harju, Force and heat current formulas for many-body potentials in molecular dynamics simulation with applications to thermal conductivity calculations, Physical Review B, 92 (2015) 094301. [50] X.L. Zhong, G.H. Furumoto, AUGMENTED BURNETT-EQUATION SOLUTIONS OVER AXISYMMETRICAL BLUNT BODIES IN HYPERSONIC FLOW, J Spacecraft Rockets, 32 (1995) 588-595. [51] R.K. Agarwal, K.Y. Yun, R. Balakrishnan, Beyond Navier-Stokes: Burnett equations for flows in the continuum-transition regime, Physics of Fluids, 13 (2001) 3061-3085. [52] J.N. Moss, G.A. Bird, Direct simulation of transitional flow for hypersonic reentry conditions (Reprinted from AIAA paper 84-0223, Jan. 1984), J Spacecraft Rockets, 40 (2003) 830-843. [53] M.N. Kogan, Flows at large Knudsen numbers, Rarefied gas dynamics, Plenum press1969. [54] R.G. Wilmoth, R.C. Blanchard, J.N. Moss, Rarefied Transitional Bridging of Blunt Body Aerodynamics,
21st
International Symposium on Rarefied Gas DynamicsMarseille, France, 1998. [55] F.J. Regan, S.M. Anandakrishnan, Dynamics of Atmospheric Re-Entry, Nasa Sti/recon Technical Report A, 93 (2015). [56] J.C. Boison, An experimental investigation of blunt body stagnation point velocity gradient, ARS Journal, 29 (1959) 130-135. [57] Y. Tao, Z. Li, Q. Han, C. Sun, Y. Zhang, Y. Chen, Theory of aerodynamic heating from molecular collision analysis, Physics Letters A, (2019) 126098. [58] B. Sudarshan, S. Saravanan, Heat Flux Characteristics Within and Outside a Forward Facing Cavity in a Hypersonic Flow, Experimental Thermal & Fluid Science, 97 (2018) 59-69.
20
Declaration of interests ܈The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. տThe authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
Quan Han, Chengdong Sun, Yi Tao, Zhongwu Li, Yan Zhang, Yunfei Chen
PII:
S1270-9638(19)33092-5
DOI:
https://doi.org/10.1016/j.ast.2019.105673
Reference:
AESCTE 105673
To appear in:
Aerospace Science and Technology
Received date:
29 November 2019
Revised date:
25 December 2019
Accepted date:
26 December 2019
Please cite this article as: Q. Han et al., Thermal protection of a hypersonic vehicle by modulating stagnation-point heat flux, Aerosp. Sci. Technol. (2020), 0, 105673, doi: https://doi.org/10.1016/j.ast.2019.105673.
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.
Thermal Protection of a Hypersonic Vehicle by Modulating Stagnation-Point Heat Flux Quan Han†, Chengdong Sun†, Yi Tao, Zhongwu Li, Yan Zhang*, Yunfei Chen* Jiangsu Key Laboratory for Design and Manufacture of Micro-Nano Biomedical Instruments, School of Mechanical Engineering, Southeast University, Nanjing 211189, China
Abstract: The heat flux at the stagnation-point is usually used to evaluate the heat transfer efficiency between high-temperature gas and a hypersonic flying vehicle. Controlling the stagnation-point heat flux is of critical important for the design of a thermal protection system for a hypersonic flying vehicle. In this work, we demonstrate that the nonequilibrium molecular dynamic (NEMD) method can well describe the shock wave and the stagnation-point heat flux. In particular, our results show that the heat flux is more sensitive to the change of the peak amplitude of the hypersonic flow in the bimodal velocity probability density function (PDF) at the stagnation point. Decreasing the peak amplitude of the hypersonic flow in the PDF can effectively reduce the stagnation-point heat flux. Three types of the outer shapes for the flying objects are employed to compare the heat flux at the stagnation-point. The NEMD simulation results confirmed that the heat flux can be decreased 16.8% for the spherical cavity outer shape compared with the flat shape in the leading edge. Keywords: stagnation-point heat flux; velocity probability density function; hypersonic flow; atom collisions
1.
Introduction A hypersonic flying vehicle has to face the challenge of aerodynamic heating as it flying in the
near atmosphere space. The heating comes from the strong shock wave formed on the leading edges of the vehicle and the viscous boundary layer attached to the whole vehicle body [1]. To avoid the damage of the vehicle, thermal protection systems (TPS) have to be deployed, including the insulated
瀗
澳Q. H. and C.S. contributed equally to the work澳 澳All correspondence should be addressed to Y. Zhang, Email: [email protected]澳or Y. Chen, Email: [email protected]; 澳 1
澽
structure, heat pipe, ablation, aerospike, cavity, jet and some combinations [2-9]. Experimental [10, 11] and numerical studies [12, 13] have proven that these TPS are effective if the flying time is short. However, for hypersonic vehicles with flying time a little bit longer, traditional TPS techniques become powerless due to the extremely high-temperature gas in front of the vehicle. In order to propose more efficient TPS techniques, the intrinsic mechanism of heat transfer between the high-temperature gases and the hypersonic flying vehicle must be figured out first. In a hypersonic flying vehicle, peak heat flux generally occurs at the stagnation point located in the region of nose and wing leading edges, which may induce enormous heat load on the vehicle. Thus, the theory of stagnation-point heat flux attracts lots of attention in the field of hypersonic thermal protection. In general, the flow state is categorized as continuous flow, transitional flow and free molecular flow by the non-dimensional parameter Knudsen number (Kn), defined as the ratio of the mean free path of gas molecules (ߣ) and the characteristic length (L). In 1950s, Lees [14], Fay and Riddell [15] proposed the stagnation-point heat flux formula by analytically solving the boundary equation simplified from the Navier-Stokes (NS) equations in the continuous flow regime (Kn < 0.01). It was estimated that the heat flux q at the stagnation point was inversely proportional to the square root of the nose radius R ( ן ݍඥሺͳȀܴሻ) [16-18]. However, in practical applications, the heat flux does not increase to infinity as predicted from the formula if the nose radius approaches zero. With the continual decrease of the nose radius, the Knudsen number will increase and be larger than 0.01. In this case, the simplified formula fails in predicting the heat flux. Furthermore, in the free molecular flow regime with Kn > 10, Regan and Anandakrishnan [19] introduced an analytical solution for the stagnation-point heat flux by solving the Boltzmann equation without the collision term, and the results indicated that the heat flux was independent of the nose radius. However, in the transitional flow regime as 0.01< Kn <10, it was much hard to obtain an analytical solution or the direct numerical solution for heat flux from the Boltzmann equation [20]. Subsequently, many simplified analytic methods developed from the Boltzmann equation were adopted, such as the linearized Boltzmann equation [21], the moment method [22] and the Bhatnagar-Gross-Krook (BGK) model [23]. In recent two decades, a new gas-kinetic scheme method developed from the BGK model for the entire Knudsen number flow, were constructed to solve the multi-scale problems from the continuous flow to the free molecular flow [24]. Although the theoretical studies of stagnation-point heat flux for continuous flow, 2
transitional flow, and free molecular flow have been applied to engineering applications, the solution for the heat flux is still not accurate enough. Particularly, at high Mach number, the governing equations including the Newton’s law for describing gas viscosity and the Fourier’s law for heat conduction may fail, which lead to the failure of the NS equations in describing the hypersonic gas flow [25]. Also, the linearized Boltzmann equation [21] is only suitable for the low speed flow problems, and the solutions for the moment method including the Burnett equation [26] and thirteen moment equations of Grad [27] are not exact at high Mach number. The widely used BGK model equation [23] is powerless in solving the problems that need to consider momentum and energy exchanges simultaneously, especially in predicting the heat flux and friction force of surface. With the growth of computer processing power, numerical methods are widely used to analyze gas flow at microscale or nanoscale level. In 1976, Bird [28] first implemented the direct simulation Monte Carlo (DSMC) method to delineate the circumstances under the condition in which the continuum model lost its validity. Up to now, the DSMC method has been considered as the most effective method for hypersonic rarefied flow [29-32]. However, the DSMC method heavily relies on the collision models, the molecular chaos hypothesis and dilute gas hypothesis [33], which makes the DSMC model fail in some cases. Besides the DSMC method, the molecular dynamic simulation (MD) is another popular method to study the microstructure of hypersonic shock wave [34], and has been widely used in studying the heat transfer mechanisms at microscale and nanoscale [35]. From the atomic or molecular scale, the MD method uses realistic empirical interatomic potential and numerically solves the Newton's second law of motion to further reduce the approximation. In 1978, Klimenko and Dremin [36] published the first non-equilibrium molecular dynamics (NEMD) simulation of shock waves in dense liquid argon, and the results of NEMD and the NS equations were in good agreement. During the past forty years, the developed NEMD methods had become an accurate tool to study the shock waves in the hypersonic flow [37]. With the help of the NEMD method, the velocity distribution can be effectively extracted and gradually become a powerful tool for studying shock waves. Although the NEMD method has well revealed the microscopic structure of the shock waves, most of the simulated models are based on the shock tube or the high-speed piston [25, 38, 39], while little researches focus on the model of a flying object in the gas environment. Intending to improve the understanding of heat transfer between the gas and a hypersonic flying vehicle at the atomic scale, this work focuses on describing the stagnation-point heat flux with the 3
NEMD method. First, a nanoscale model of a piece of silicon cuboid flying in an argon gas environment is established, and the density, temperature, heat flux and velocity distribution of the shock waves formed over the cuboid can be investigated. Subsequently, by changing the cross-section size of the silicon cuboid, the stagnation-point heat fluxes at different Kn numbers are studied through the atomic scale method and classical aerodynamic analysis methods. Especially, through comparing the bimodal velocity distributions at the stagnation point at different Kn numbers, the influence of the velocity distribution on the heat flux is studied. Next, simulations under different environmental atmospheric pressures and different geometry shapes are conducted to study the factors affecting the velocity distribution and the heat flux. Based on the NEMD model, the outer geometry of a hypersonic flying silicon cuboid is optimized to modulate the stagnation-point heat flux, which provides a technique route to reduce the heat load for a hypersonic vehicle. 2.
Methodology
2.1 Details of the NEMD simulation In this work, a silicon cuboid flying with a hypersonic velocity through an argon gas environment is simulated using the NAMD 2.12 package [40]. Neutral atoms or molecules can be accurately modeled with the 12-6 Lennard-Jones (LJ) potential [41].
σ
σ
VLJ (rij ) = 4ε (( )12 − ( )6 ) rij rij
(1)
where ɂ andɐare the depth of the potential well and distance at which the LJ interaction is zero, and rij is the distance between two particles. For argon atoms, ɂ ൌ ͲǤʹ͵ͺͳkcal/mole andɐ ൌ ͵ǤͶͳͳͺՀ [42]. CHARMM36 force fields [43] are used to describe the interactions between silicon atoms, and the nonbonded interaction for silicon atoms is described by the LJ potential with the parameters ɂ ൌ ͲǤ͵ͳͲͲ kcal/mole and ɐ ൌ ͵ǤͺͲͷՀ [44]. Moreover, the nonbonded interactions between two different particles are calculated using the mixing rules [45]. In our simulations, the cutoff radius is set as 8Հ.
4
Figure 1. Schematic diagram of a silicon cuboid flying in an argon gas environment along the z-direction.
First, all simulations for studying heat flux at different Kn numbers, pressures, as well as the geometry shapes, are performed on the same system. The schematic diagram of the system is displayed in Fig. 1, in which the silicon cuboid flies along the z-direction surrounded by the gaseous argon. Periodic boundary conditions are applied along the three axis directions. For the cases at different Kn numbers, the initial states of the argon gas environment are set as p0 = 10 atm and T0 = 300 K, and the mean free path of the gas is predicted to be 7 nm [46]. For silicon, its lattice constant is fixed at the optimized value of 5.43Հ. Table 1 lists the details of each simulation when the cross-section size of the silicon cuboid varies from 4 to 32 unit cells (UC), and the corresponding global Kn number varies from 3.2 to 0.4 at the tenfold atmospheric pressure. For the cases at different pressures and geometry shapes, the details are described in section 3.2 and section 3.3. Table 1. Size of the silicon cuboid, Kn number, box size, and the number of atoms Case
Pressure
Silicon
Global
Box
Silicon
Argon
No.
(atm)
Size (UC3)
Kn number
Size (Հଷ )
Number
Number
1
10
ͶൈͶൈͶ
3.2
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
512
25969
2
10
ൈൈͶ
2.1
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
1152
25964
3
10
ͺൈͺൈͶ
1.6
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
2048
25958
4
10
ͳʹ ൈ ͳʹ ൈ Ͷ
1.1
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
4608
25940
5
10
ͳ ൈ ͳ ൈ Ͷ
0.8
͵ͲͲ ൈ ͵ͲͲ ൈ ͳͲͲͲ
8192
25915
6
10
͵ʹ ൈ ͵ʹ ൈ Ͷ
0.4
ͺͲͲ ൈ ͺͲͲ ൈ ͳͲͲͲ
32768
168197
For each case, the velocity Verlet algorithm with a time step of 1 fs is used for the integration of the motion equations [47]. The system is first minimized for 0.1ns and then equilibrated in an isothermal-isobaric (NPT) ensemble for 2 ns. Then in the canonical (NVT) ensemble, the silicon cuboid is pulled along the positive z-direction with a speed of 2000 m/s by a driving spring, named the 5
Steered Molecular Dynamics method [48]. In order to keep the flying direction along z-direction, the silicon cuboid atoms are harmonically restrained in x and y directions. The argon atoms and silicon atoms are coupled with Langevin thermostats at 300 K with different damping coefficients. The simulation is carried out for at least 7 ns to guarantee the system reaching a steady state, and then the heat flux is collected to certify the convergence of the NEMD calculations. In this work, a series of simulations have been performed to optimize the parameters in constructing the simulation systems (Figure. S1, Supplementary material). Macroscopic observables, including density, temperature, and heat flux are calculated as the spatial and time average in the volume of a prescribed grid cell, with a size defined asȟ ݔൈ ȟ ݕൈ ȟݖ. When plotting the curves of macroscopic parameters along the z-direction, we set the length of ȟ ݖas 7Հ and the size of ȟ ݔൈ ȟ ݕas the cross-section size of the silicon cuboid in the x-y plane (consistent with the grid definition of the stagnation-point region, Section S2, Supplementary material). When drawing the contour of density in the x-z plane, the size of ȟ ݔൈ ȟ ݖis set as ͳ ൈ ͳՀଶ and the length of ȟ ݕis set as the length of the silicon cuboid in the y-direction. In the MD simulations, the local heat flux can be calculated as the time derivative of the sum of the site energy moments over the atoms in each grid cell [49].
G q=
1 d N JG ¦ ri Ei ΔxΔyΔz dt i =1
(2)
N 2 1 JG Ei = m vi − v + ¦Vij 2 j =1
(
)
(3)
where ݒҧ is the averaged velocity of all the atoms in each grid cell, ሬሬሬԦ ݒప is the velocity of each atom i, and Vij is the potential energy of atom i due to atom j. The formula of the local heat flux is finally derived as: G q=
N JJG JG º ª JG 1 1 N JJG ª JG − + E v v f ij < vi − v + v j − v º
(
( (
)
) (
) )
(4)
డೕ where ሬሬሬԦ ݂ ൌ ሬሬሬሬሬԦ is the force on atom i due to atom j. డഢണ
2.2 Model verification To verify that our simulation setup can effectively depict the essential physics in a hypersonic flow, 6
the number density, temperature, and heat flux profile of the gas flow in the steady state are extracted from the NEMD model and displayed in Fig. 2. In this verification case, a silicon cuboid with size of
ͳ ൈ ͳ ൈ Ͷଷ flying in the gaseous argon with a speed v0 of 2000 m/s (about 6.2 Mach number) is simulated, which is listed as case No. 5 in Table. 1. As expected, a shock wave is generated in front of the flying object demonstrated as the normalized number density contour in the inset of Fig. 2(a). In the shock wave, an abrupt change in normalized number density of the gaseous argon can be observed from Fig. 2(a). The normalized number density reaches the peak in position (Ϩ) near the silicon cuboid front surface. With the increase of the distance away from the silicon cuboid front surface, the normalized number density of the gaseous argon starts to decrease. At position (IV), the number density of the gaseous argon decreases to the environmental argon density. This profile of the normalized number density of the gaseous argon matches qualitatively well with the numerical solutions of the augmented Burnett equation [50, 51] and the DSMC results [52]. Besides, the temperature curve of the gaseous argon is provided in Fig. 2(b), and a peak value occurs in the position (ϩ) of the shock wave. From the heat flux computed by Eq. (4), it can be seen that a negative heat flux locates at the left of the position (ϩ) in Fig. 2(b), meaning that heat is transferred from the shock wave to the silicon cuboid, while the positive heat flux indicates that heat is transferred from the shock wave to the gas far away from the silicon cuboid. Both the temperature gradient and heat flux direction indicate the rationality of the simulation model.
7
Figure 2. The NEMD simulation results. (a) Normalized number density n/n0, where n denotes the number density inside each grid and n0 is the initial gaseous number density. (b) Temperature and heat flux calculated in NEMD are given as functions of distance in the shock wave. The normalized number density contour shows a shock wave generated before the hypersonic flying silicon cuboid with a speed of v0 = 2000 m/s. The grey region stands for the silicon cuboid. (c) The probability density function of argon atom velocity in the z-direction in four different positions: near the front surface (Ϩ), shock wave (ϩ) and (Ϫ), and far away from the front surface (ϫ). The black dots stand for the results extracted from the NEMD simulation and the red curves are predicted from the Mott-Smith theory.
Notably, the velocity distribution of the argon atoms in the shock wave does not obey the Maxwell distribution due to the high nonequilibrium flow. The Mott-Smith theory [22] describes a shock wave as a superposition of a hypersonic flow and a subsonic flow, and their proportions vary with the position along the z-direction, which is named as the bimodal velocity distribution function:
f = fs + f h
(5)
G JJG m(v − vα ) 2 m 3/2 fα = nα ( z )( ) exp( − ) 2π k BTα 2k BTα
(6)
where
where f is the velocity distribution function and the subscript s and h denote the subsonic flow and 8
hypersonic flow, respectively. The ߙ represents the subscript s or h, and the݊ఈ ሺݖሻdenotes the number density varied with the position along the z-direction. The gaseous temperature ܶఈ and velocity of gas flowݒ ሬሬሬሬԦ ఈ are assumed to be independent of z, and they are calculated through the Rankine-Hugoniot relations [1]. Fig. 2(c) displays the velocity distribution of argon atoms in the z-direction from the NEMD simulation and the Mott-Smith theory in the four positions (Ϩ), (ϩ), (Ϫ) and (ϫ) indicated in Fig. 2(a). In order to describe the velocity PDF more obviously, we consider the flying silicon cuboid as the frame of reference, in which the velocity for the argon atoms move toward the cuboid is defined to be positive, and the velocity is normalized by the speed v0 of silicon cuboid. In NEMD, the velocity distribution is given as a probability density function (PDF), defined as fv = ni /(n¨v), and where ni denotes the number density in the velocity range from vi/v0 - ¨v/2 to vi/v0 + ¨v/2 and n denotes the number density inside each grid. In the Mott-Smith theory, the velocity PDF is obtained by the velocity distribution function f in the z-direction, which is properly normalized so that ݂ ௩ ݀ݒ௭ over all velocities equals one. First, a nearly undisturbed hypersonic flow locates in the position (ϫ), and its velocity PDF is a single sharp peak centered at the velocity v0 of 2000 m/s. As the flow approaches the silicon cuboid, the argon atoms in position (Ϫ) experience a few collisions with the atoms reflected back by the cuboid obstacle and a subsonic flow appears in the distribution. As expected, more collisions with the reflected atoms occur in position (ϩ), resulting in a higher peak amplitude of subsonic flow. When it comes to the position (Ϩ) near the silicon cuboid front surface, the initial hypersonic flow is almost replaced by the subsonic flow due to the numerous atom collisions. These observations show that when the initial hypersonic flow gradually approaches the object, the flow in shock wave gradually becomes a superposition of the hypersonic flow and subsonic flow, which is called bimodal velocity distribution. However, there is an obvious quantitative difference in the subsonic flow between the NEMD results and the Mott-Smith theory in Fig. 2(c). In the Mott-Smith theory, the temperature and velocity of subsonic flow behind a shock wave are calculated through the Rankine-Hugoniot relations, which are based on the three conservation equations of mass, momentum, energy and the equation of state for ideal gases. Actually, the Rankine-Hugoniot relations fail in accurately depicting the parameters behind a shock wave in this case due to the highly nonequilibrium flow. But it is worth noting that the 9
overall tendency of velocity PDF obtained from NEMD simulation is quite reasonable in Fig. 2(c), further validating the feasibility of the NEMD method in simulating the hypersonic flow.
3.
Results and discussion
3.1 The effect of Kn number on the stagnation-point heat flux In our data analysis process, the grid length is a critical parameter to calculate the macroscopic parameters such as density, temperature and heat flux. To determine the optimal lengthȟݖnear the front surface, we perform a series of statistical analysis with different grid lengths ranging from 3 to 11Հ (Figure. S2, Supplementary material). Length of 7 Հ is chosen for further study, because it can ensure a sufficient number of atoms in a statistical grid cell and is small enough to guarantee the accuracy of the definition that considering position (Ϩ) as the stagnation point in the immediate neighborhood of the silicon cuboid front surface, besides, it produces a convergent heat flux. To investigate the effect of Kn number on the heat flux in the stagnation-point region (i.e., stagnationpoint heat flux), we simulate six cases of different Kn numbers with an atmospheric pressure of 10 atm. In our model, the different Kn number ( ݊ܭൌ ߣȀ )ܮis obtained through changing the cross-section size L of the silicon cuboid as shown in Table. 1.
3.1.1 Details in the stagnation-point region First, we study the influence of the Knudsen number Kn on the physical parameters of gaseous argon in the stagnation-point region, including the number density and velocity distribution. During the flight of the silicon cuboid, more atom collisions between the free-stream flow and the silicon cuboid front surface occur when the cross-section size becomes larger. Thus, more argon atoms are reflected back to the face coming flow, resulting in a larger number density ns at the stagnation point in the steady state. With the Kn number increasing from 0.4 to 3.2, the number density decreases at the stagnation point as listed in Fig. 3. For example, as the Knudsen number Kn increases to 3.2, the number density at the stagnation point decreases to 0.76/nm3. It is believed that the events of the atom collision become less in front of the flying cuboid silicon at the high Kn number, which leads to the decrease of the number density. Similarly, the phenomenon that the events of the atom collision become less also significantly affects the velocity PDF, in which the peak amplitude of the hypersonic flow gradually increases at the larger Knudsen number, and the peak amplitude of subsonic flow 10
decreases correspondingly. It can be predicted that when the Kn number increases infinitely, that is, the cross-section size of silicon cuboid decreases to zero, the velocity PDF will be the same as that of the free-stream flow with a single hypersonic peak.
Figure 3. The velocity PDF of argon atoms in the z-direction at the stagnation point (Ϩ) as the Kn numbers vary from 0.4 to 3.2. The number density ns at the stagnation point are listed in the upper right corner of each figure. The dots represent the data extracted from the NEMD simulations and the curves are the fitted bimodal velocity PDF.
3.1.2 Heat flux at the stagnation point Several models are used to calculate the heat flux at the stagnation-point. If the velocity distribution function f in Eq. (5) is known, various macroscopic parameters including the heat flux at the stagnation point can be accurately obtained by calculating the moment of f. In the method of the “moment of f ”, the flux of translational energy Ei brought to the silicon front surface per area in unit time by all the incident argon atoms is first calculated by the velocity distribution function. Then the flux of translational energy Er carried away by those reflected atoms is computed by the velocity distribution function for the reflected atoms, which is determined by the law of interaction of the argon atoms with the surface. Finally, the heat flux to the silicon cuboid front surface can be expressed as q = Ei – Er [53]. Besides the model of the moment equation, the method of bridging function is also 11
employed to estimate the heat flux from continuous flow to free molecular flow. By using the classical aerodynamic solutions of the stagnation-point heat flux for continuous flow and free molecular flow, the bridging function can be obtained at the entire Kn number by fitting the free parameters a1, a2 and a3 to the NEMD results. The bridging function is written as a function of Kn number [54] :
qb = qFM × Pb + (1 − Pb ) × (qc / (1 − a3 × Kn))
(7)
Pb = sin(ϕ )
(8)
ϕ = π (a1 + a2 × log10 ( Kn))
(9)
where a1, a2, a3 are the fitting parameters. qFM is the analytical solution of the heat flux for free molecular flow [55]. qc is the heat flux calculated from the Fay-Riddell model solving the NS equations in the continuum flow [15], and the effective radius used in the Fay-Riddell model for the silicon cuboid can be computed from ref. [56].
Figure 4. Heat flux coefficients defined as h = q/qFM are given as a function of the Kn number obtained from several theoretical models and the NEMD. The NEMD data is computed from Eq. (4). The model of “moment of f ” denotes that the heat flux is calculated from the fitted velocity distribution f , and the proportions of the heat fluxes contributed by hypersonic flow and subsonic flow at different Kn numbers are presented in the inset. The bridging function is shown as Eq. (7) and the Fay-Riddell model denotes that the stagnation-point heat flux is solved from the NS equations. Free molecular flow model means that the heat flux is calculated from the Boltzmann equation. (Argon gas, T0 = 300K; Silicon cuboid, T0 = 300K, v0 = 2000 m/s, the sizes are shown in Table. 1).
When studying the influence of the Kn number on the stagnation-point heat flux, the heat flux is generally normalized by the value for the free molecular flow, called the heat flux coefficient h. In Fig. 12
4, we compare the coefficient h obtained from the NEMD with the calculations by other methods, in which the green dashed line with a value of one denotes the results for free molecular flow. Due to the limited computing resource, the minimum Kn number of the NEMD models is 0.4. As can be seen in Fig. 4, the heat flux coefficient predicted by all the methods increases with the Kn number, which is consistent with the previous report that the smaller cross-section size with larger Kn number contributes to the larger heat flux [57]. However, the result of the Fay-Riddell model is significantly different from the NEMD results, which indicates the failure of the linear law in computing the heat flux when it comes to the transitional flow. On the contrary, the heat flux coefficient from the bridging function agrees well with the NEMD results in Fig. 4, with the fitting parameters in the expression as a1 = 0.12, a2 = 0.32, a3 = 0.8 in Eq. (7). This consistency shows that the NEMD results can be well depicted by the bridging function. It is worth noting that the results of the “moment of f ” are almost identical to the NEMD results in Fig. 4, indicating that the heat flux can be accurately derived from the velocity distribution function. Actually, this velocity distribution function f is obtained by the fitted velocity PDF fv multiplied by the number density ns in Fig. 3. To further understand the roles of hypersonic and subsonic flow in the stagnation-point heat flux, we study the proportions of heat flux contributed by hypersonic and subsonic flow with the model of “moment of f ”. As shown in the inset of Fig. 4, the heat flux contributed by hypersonic flow increases significantly with the Kn number, and the heat flux contributed by subsonic flow declines rapidly. And in the distribution function f , the peak amplitude of the hypersonic flow still become higher although the number density ns decreases significantly at the larger Kn number (Figure S3, Supplementary material). Therefore, an increase in the stagnationpoint heat flux is observed due to the larger increment of heat flux contributed by the hypersonic flow. In addition, although the peak of the hypersonic flow in Fig. 3 (a) is extremely tiny, the proportion of heat flux contributed by hypersonic flow still reaches to 37% when Kn number equals 0.4. When the Kn number rises to 3.2, the proportion of heat flux contributed by hypersonic flow is up to 91% while the peak of the hypersonic flow in the velocity PDF is still lower than the peak of subsonic flow. These results show that the stagnation-point heat flux is more sensitive to the change of the peak amplitude of the hypersonic flow. Based on this finding, a design criterion is proposed, i.e. an efficient way to reduce the heat flux can be realized by decreasing the peak amplitude of the hypersonic flow in the velocity PDF. 13
3.2 Effect of pressure on velocity distribution and heat flux at the stagnation point In order to investigate the factors affecting the velocity distribution, we carry out the NEMD simulations of silicon cuboid with the size of ͳ ൈ ͳ ൈ Ͷ ଷ under different environmental atmospheric pressures. As can be observed in Fig. 5(a), the normalized number density n/n0 decreases rapidly to 1.0 along the flying direction. Specifically, the normalized number density at the stagnation point decreases rapidly with the increase of the environmental atmospheric pressure, while only a small difference can be seen when z is larger than 20 Հ . Moreover, in the inset of Fig. 5(a), the peak amplitudes (v/v0=1) in the velocity PDF are almost the same under different atmospheric pressures. It is believed that, during the flight of silicon cuboid under different atmospheric pressures, the number of argon atoms reflected back after colliding with the front surface increases proportionally when the number density n0 of the free-stream flow increases with the pressure, which leads to the same influence on the face coming gas flow. Therefore, we infer that the peak amplitude of the hypersonic flow in the velocity PDF is independent of atmospheric pressure but dependent on the cross-section size of the silicon cuboid.
Figure 5. (a) Normalized number density is given as a function of distance in the z-direction. The size of silicon cuboid is ͳ ൈ ͳ ൈ Ͷ ଷ and the initial atmospheric pressures are 1 atm, 5 atm, 10 atm, and 20 atm respectively. The yellow region stands for the stagnation-point region. The number density and velocity PDF at the stagnation point are shown in the inset. The stagnation-point heat flux as a function of different atmospheric pressure is presented in (b). The inset shows the heat flux q from NEMD simulation versus the number density ns at the stagnation point.
As can be seen in Fig. 5(b), the magnitude of the stagnation-point heat fluxes increase with the 14
pressure. Especially, compared with the linearly increasing heat flux for free molecular flow, the NEMD result increases slowly with the environmental atmospheric pressures. However, at low atmospheric pressure, the gap between the heat flux predicted from the free molecular flow and that extracted from the NEMD becomes small because the flow at the low pressure can be approximated as free molecular flow. We also find that the Fay-Riddell model overestimate the heat flux compared with the NEMD results, and this deviation shows the failure of NS equations in capturing the heat flux in the transitional flow regime. Furthermore, the fitted bridging function still agrees well with the NEMD results at these atmospheric pressures. Besides, it is notable that the heat flux is almost linearly proportional to the number density at the stagnation-point as shown in the inset of Fig. 5(b).
3.3 Modulating stagnation-point heat flux by the outer shape of the leading edge In the previous sections, we find that an efficient method for reducing the stagnation-point heat flux is through decreasing the peak amplitude of the hypersonic flow, which is closely related to the atom collisions between the free-stream flow and the silicon cuboid front surface. Since recent researches have shown that the stagnation-point heat flux can be significantly reduced by introducing an axial cavity in the nose region of a hypersonic vehicle [58], we implement three models of the silicon cuboids with different outer shapes on the front head to verify whether the heat flux is proportional to the peak amplitude of the hypersonic flow. The shapes of silicon cuboid are shown in the inset of Fig. 6(a) as S1, a silicon cuboid with a sphere cavity on the front surface, S2, a silicon cuboid with a flat surface, and S3, a silicon cuboid with a sphere front surface, in which the depth for both the sphere cavity and sphere is set as Z* equaled to 2 UC. The statistical grid is shown as the region with red dashed lines in the inset of Fig. 6(a), which is in front of the silicon surface within a length ȟ ݖequaled to 7Հ.
15
Figure 6. (a) Normalized number density of gas atoms is given as a function of distance in the z-direction. The atmospheric pressure is set at 10 atm. The yellow region stands for the stagnation-point region. The flying object is a silicon cuboid with cross-section size of ͳ ൈ ͳ ଶ and having different outer shapes in the head. In the inset of (a), the magnitude of heat flux at the stagnation point from the NEMD simulations are presented. The velocity PDF are presented in (b), and the inset shows the velocity distribution function f of the hypersonic flow.
First, the normalized number density n/n0 for the flying silicon cuboid with three outer shapes are shown in Fig. 6(a). The normalized number density at the shape S1 of a sphere cavity is the largest, while the value at the shape S3 of a sphere is the smallest. In fact, during the flight of the silicon cuboid, the atom collisions between the free-stream flow and the front surface at both cases of S1 and S3 are stronger than that at the case S2 due to the larger surface area. However, more argon atoms slip away the object after colliding with the spherical surface at the case S3, so fewer argon atoms are reflected back to the shock wave layer, resulting in the smallest normalized number density. On the contrary, more argon atoms can be reflected back to form the shock wave layer at the case S1 because of the cavity, so the normalized number density is the largest. Similarly, the peak amplitude of the hypersonic flow at the shape S1 is the lowest because large amount of argon atoms are reflected back, and the peak amplitude of the hypersonic flow in the velocity PDF increases obviously from the shape S1 to S3 as shown in Fig. 6(b). In particular, despite the decreasing of stagnation-point number density from the shape S1 to the shape S3, the peak amplitude of the hypersonic flow still increases in the velocity distribution function f as shown in the inset of Fig. 6(b). Thus, the magnitude of the stagnation-point heat flux increases from the shape S1 to the shape 16
S3 as shown in the inset of Fig. 6(a). And compared to the shape of the flat surface, the stagnationpoint heat flux is reduced by 16.8% due to the decreased peak amplitude of the hypersonic flow at the shape of the sphere cavity, indicating that the cavity shape can be used to reduce the heat flux.
4.
Conclusion In summary, the NEMD simulations are carried out to study the stagnation-point heat flux of the
hypersonic vehicle at nanoscale. Through analyzing the gaseous number density, temperature, heat flux and velocity probability density function along the flight direction, it is validated that the hypersonic flow can be well described by the NEMD simulations. Through the NEMD simulations, it is found that the heat flux increases with the Kn number, which is attributed to the peak amplitude rise of the hypersonic flow in the velocity PDF. Furthermore, when the peak amplitude of the hypersonic flow is the same while the peak amplitude of subsonic flow is different in the velocity PDF, the stagnation-point heat flux increases almost linearly with the stagnation-point number density. Based on the results, a conclusion is derived that the stagnation-point heat flux is more sensitive to the change of the peak amplitude of the hypersonic flow in the velocity PDF due to the much large translational energy stored in the hypersonic flow. Thus, a design criterion is proposed to reduce the stagnation-point heat flux, i.e. decreasing the peak amplitude of the hypersonic flow in the velocity PDF. Three type outer shapes for the leading edge of a hypersonic flying silicon cuboid are employed to verify the design criterion. In addition, we investigate the factors that affect the peak amplitude of the hypersonic flow in the velocity PDF. During the hypersonic flight of the silicon cuboid, argon atoms are reflected back by the leading edge and collide with the face coming flow, resulting in the formation of the shock wave. The results indicate that the peak amplitude of the hypersonic flow in the velocity PDF is independent of the atmospheric pressures but dependent on the geometry size of silicon cuboid. Our study provides a method to decrease the peak amplitude of hypersonic flow to reduce heat flux, which will be helpful to innovate the TPS design at nanoscale in a hypersonic flying vehicle.
Acknowledgement The authors thank the financial support from the National Key R&D Program of China (2018YFB1105400) and the Natural Science Foundation of China (Grants No.51435003, 51575104, 17
51705074).
References [1] A. Viviani, G. Pezzella, Aerodynamic and aerothermodynamic analysis of space mission vehicles, Springer press2015. [2] W. Huang, L. Yan, J. Liu, L. Jin, J.-g. Tan, Drag and heat reduction mechanism in the combinational opposing jet and acoustic cavity concept for hypersonic vehicles, Aerospace Science and Technology, 42 (2015) 407-414. [3] P. Mashaei, M. Shahryari, S. Madani, Analytical study of multiple evaporator heat pipe with nanofluid; A smart material for satellite equipment cooling application, Aerospace Science and Technology, 59 (2016) 112-121. [4] Q. Qin, J. Xu, Numerical evaluation of aerodome and cooling jet for aeroheating reduction, Aerospace Science and Technology, 86 (2019) 520-533. [5] M. Rivier, J. Lachaud, P.M. Congedo, Ablative thermal protection system under uncertainties including pyrolysis gas composition, Aerospace Science and Technology, 84 (2019) 1059-1069. [6] L. Zhu, Y. Li, X. Chen, H. Li, W. Li, C. Li, Hypersonic flow characteristics and relevant structure thermal response induced by the novel combined spike-aerodome and lateral jet strategy, Aerospace Science and Technology, 95 (2019) 105459. [7] P.Q. Elias, N. Severac, J.M. Luyssen, Y.B. Andre, I. Doudet, B. Wattellier, J.P. Tobeli, S. Albert, B. Mahieu, R. Bur, A. Mysyrowicz, A. Houard, Improving supersonic flights with femtosecond laser filamentation, Science Advances, 4 (2018). [8] S. Li, W. Huang, J. Lei, Z. Wang, Drag and heat reduction mechanism of the porous opposing jet for variable blunt hypersonic vehicles, International Journal of Heat and Mass Transfer, 126 (2018) 1087-1098. [9] R.-r. Zhang, W. Huang, L.-q. Li, L. Yan, R. Moradi, Drag and heat flux reduction induced by the pulsed counterflowing jet with different periods on a blunt body in supersonic flows, International Journal of Heat and Mass Transfer, 127 (2018) 503-512. [10] Z.-g. Wang, X.-w. Sun, W. Huang, S.-b. Li, L. Yan, Experimental investigation on drag and heat flux reduction in supersonic/hypersonic flows: a survey, Acta Astronautica, 129 (2016) 95-110. [11] M. Ou, L. Yan, W. Huang, S.-b. Li, L.-q. Li, Detailed parametric investigations on drag and heat flux reduction induced by a combinational spike and opposing jet concept in hypersonic flows, International Journal of Heat and Mass Transfer, 126 (2018) 10-31. [12] X. Sun, W. Huang, M. Ou, R. Zhang, S. Li, A survey on numerical simulations of drag and heat reduction mechanism in supersonic/hypersonic flows, Chinese Journal of Aeronautics, (2019). [13] W. Huang, Z. Chen, L. Yan, B.-b. Yan, Z.-b. Du, Drag and heat flux reduction mechanism induced by the spike and its combinations in supersonic flows: A review, Progress in Aerospace Sciences, (2018). [14] L. Lees, LAMINAR HEAT TRANSFER OVER BLUNT-NOSED BODIES AT HYPERSONIC FLIGHT SPEEDS, Jet Propulsion, 26 (1956) 259-&. [15] J.A. Fay, F.R. Riddell, THEORY OF STAGNATION POINT HEAT TRANSFER IN DISSOCIATED AIR, Journal of the Aeronautical Sciences, 25 (1958) 73-&. [16] N.H. Kemp, Laminar heat transfer around blunt bodies in dissociated air, Journal of the aerospace sciences, 26 (1959) 421-430. [17] R.B. Pope, Stagnation-point convective heat transfer in frozen boundary layers, AiAA journal, 6 (1968) 619-626. [18] K. Sutton, R.A. Graves Jr, A general stagnation-point convective heating equation for arbitrary gas mixtures, 1971. [19] F.J. Regan, S.M. Anandakrishnan, Dynamics of atmospheric re-entry, American Institute of Aeronautics and Astronautics press1993. [20] N. Singh, T.E. Schwartzentruber, Heat flux correlation for high-speed flow in the transitional regime, Journal of Fluid Mechanics, 792 (2016) 981-996. [21] E.P. Gross, E.A. Jackson, Kinetic models and the linearized Boltzmann equation, The physics of fluids, 2 (1959) 43218
441. [22] H.M. Mott-Smith, The solution of the Boltzmann equation for a shock wave, Physical Review, 82 (1951) 885. [23] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Physical review, 94 (1954) 511. [24] K. Xu, A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method, Journal of Computational Physics, 171 (2001) 289-335. [25] E. Salomons, M. Mareschal, Usefulness of the Burnett description of strong shock waves, Physical review letters, 69 (1992) 269. [26] D. Burnett, The distribution of molecular velocities and the mean motion in a nonϋuniform gas, Proceedings of the London Mathematical Society, 2 (1936) 382-435. [27] H. Grad, On the kinetic theory of rarefied gases, Communications on pure and applied mathematics, 2 (1949) 331407. [28] G.A. Bird, Molecular gas dynamics, NASA STI/Recon Technical Report A, 76 (1976). [29] G. Dahlen, M. Macrossan, C. Brundin, J. Harvey, Blunt cones in rarefied hypersonic flow: Experiments and MonteCarlo simulations, (1984). [30] G. Pham-Van-Diep, D. Erwin, E. Muntz, Nonequilibrium molecular motion in a hypersonic shock wave, Science, 245 (1989) 624-626. [31] P. Shufflebotham, T. Bartel, B. Berney, Experimental validation of a direct simulation by Monte Carlo molecular gas flow model, Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures Processing, Measurement, and Phenomena, 13 (1995) 1862-1866. [32] J. All-ccedil, gre, D. Bisch, J. Lengrand, Experimental rarefied density flowfields at hypersonic conditions over 70degree blunted cone, J Spacecraft Rockets, 34 (1997) 714-718. [33] S.M. Deshpande, P.S. Raju, Monte carlo simulation for molecular gas dynamics, Sadhana, 12 (1988) 105-123. [34] B. Holian, Atomistic computer simulations of shock waves, Shock waves, 5 (1995) 149-157. [35] T. Luo, G. Chen, Nanoscale heat transfer–from computation to experiment, Physical Chemistry Chemical Physics, 15 (2013) 3389-3412. [36] V. Klimenko, A. Dremin, Detonatsiya, Chernogolovka, Akademik Nauk Moscow, 1978. [37] P. Valentini, P.A. Tump, C. Zhang, T.E. Schwartzentruber, Molecular dynamics simulations of shock waves in mixtures of noble gases, J Thermophys Heat Tr, 27 (2013) 226-234. [38] F. Uribe, R. Velasco, L. Garcia-Colin, Burnett description of strong shock waves, Physical review letters, 81 (1998) 2044. [39] W.G. Hoover, C.G. Hoover, Shockwaves and Local Hydrodynamics: Failure of the Navier-Stokes Equations,
New
Trends In Statistical Physics: Festschrift in Honor of Leopoldo García-Colín's 80th Birthday, World Scientific2010, pp. 1526. [40] J.C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R.D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, Journal of Computational Chemistry, 26 (2005) 1781-1802. [41] J.E. Jones, On the determination of molecular fields.—II. From the equation of state of a gas, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 106 (1924) 463-477. [42] M.P. Allen, D.J. Tildesley, Computer simulation of liquids, Oxford university press2017. [43] A.D. MacKerell Jr, D. Bashford, M. Bellott, R.L. Dunbrack Jr, J.D. Evanseck, M.J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, All-atom empirical potential for molecular modeling and dynamics studies of proteins, The journal of physical chemistry B, 102 (1998) 3586-3616. [44] J.A. Wendel, W.A.G. Iii, The Hessian biased force field for silicon nitride ceramics: Predictions of thermodynamic and mechanical properties for ¢and £· i3N4, Journal of Chemical Physics, 97 (1992) 5048-5062. [45] Z. Xu, M.J. Buehler, Nanoengineering heat transfer performance at carbon nanotube interfaces, ACS nano, 3 (2009) 19
2767-2775. [46] G. Chen, Nanoscale energy transport and conversion: a parallel treatment of electrons, molecules, phonons, and photons, Oxford University Press2005. [47] D. Frenkel, B. Smit, Understanding molecular simulation: From algorithms to applications, 1996. [48] B. Isralewitz, J. Baudry, J. Gullingsrud, D. Kosztin, K. Schulten, Steered molecular dynamics investigations of protein function, Journal of Molecular Graphics & Modelling, 19 (2001) 13-25. [49] Z. Fan, L.F.C. Pereira, H.Q. Wang, J.C. Zheng, D. Donadio, A. Harju, Force and heat current formulas for many-body potentials in molecular dynamics simulation with applications to thermal conductivity calculations, Physical Review B, 92 (2015) 094301. [50] X.L. Zhong, G.H. Furumoto, AUGMENTED BURNETT-EQUATION SOLUTIONS OVER AXISYMMETRICAL BLUNT BODIES IN HYPERSONIC FLOW, J Spacecraft Rockets, 32 (1995) 588-595. [51] R.K. Agarwal, K.Y. Yun, R. Balakrishnan, Beyond Navier-Stokes: Burnett equations for flows in the continuum-transition regime, Physics of Fluids, 13 (2001) 3061-3085. [52] J.N. Moss, G.A. Bird, Direct simulation of transitional flow for hypersonic reentry conditions (Reprinted from AIAA paper 84-0223, Jan. 1984), J Spacecraft Rockets, 40 (2003) 830-843. [53] M.N. Kogan, Flows at large Knudsen numbers, Rarefied gas dynamics, Plenum press1969. [54] R.G. Wilmoth, R.C. Blanchard, J.N. Moss, Rarefied Transitional Bridging of Blunt Body Aerodynamics,
21st
International Symposium on Rarefied Gas DynamicsMarseille, France, 1998. [55] F.J. Regan, S.M. Anandakrishnan, Dynamics of Atmospheric Re-Entry, Nasa Sti/recon Technical Report A, 93 (2015). [56] J.C. Boison, An experimental investigation of blunt body stagnation point velocity gradient, ARS Journal, 29 (1959) 130-135. [57] Y. Tao, Z. Li, Q. Han, C. Sun, Y. Zhang, Y. Chen, Theory of aerodynamic heating from molecular collision analysis, Physics Letters A, (2019) 126098. [58] B. Sudarshan, S. Saravanan, Heat Flux Characteristics Within and Outside a Forward Facing Cavity in a Hypersonic Flow, Experimental Thermal & Fluid Science, 97 (2018) 59-69.
20
Declaration of interests ܈The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. տThe authors declare the following financial interests/personal relationships which may be considered as potential competing interests: