Abstract
Numerical studies of microscale shock-vortex interaction were conducted by particle-based direct simulation of Monte Carlo (DSMC). The enstrophy is found to be increased in the strong microscale shock-vortex interaction, which is not observed in the previous DSMC studies within the limited cases. Investigations also show that the increase of the enstrophy results in an increase in dissipation rate during the strong interaction. The incoming Mach number, vortex size, and vortex Mach number turn out to play a critical role in the strength of interaction, which in turn govern the change in the dissipation rate and the increase or decrease in enstrophy during the microscale shock-vortex interaction. It is also observed that the incoming Mach number is the most dominant parameter, followed by vortex size and vortex Mach number, during the microscale shock-vortex interaction.
1. Introduction
Shock-vortex interactions have been extensively investigated in the past decades since they can provide fundamental knowledge on flow physics of high speed gas flows. Most of the previous studies on this problem have been conducted for the purpose of understanding the noise production mechanism as well as deformation of the shock wave and vortex. The sound wave structures are found successively by Ellzey et al. [1], Inoue and Hattori [2], and Zhang et al. [3]. In addition, Chatterjee and Vijayaraj [4] captured multiple acoustic waves, quadrupolar in nature and with successive layout of phase. Later, Chang et al. [5] studied various cases of shock-vortex interactions. Then, Chaudhuri et al. studied the complex shock-obstacle interactions through the immersed boundary method, which involved the typical characteristics of shock-vortex interaction [6]. In these studies, the induced expansion wave and shock wave are found to be the most remarkable flow elements in the shock-vortex interaction.
Although numerous studies have been conducted on the macroscale shock-vortex interaction, the shock-vortex interaction in microscale is not yet well understood. While a few cases had been conducted on the compression effect or turbulent in microchannel flow using DSMC method [7–9], unfortunately, these studies did not consider the interaction of vortex and shock wave. And, up to now, only a few microscale cases involving limited Mach and Knudsen numbers have been conducted by Koffi et al. using DSMC [10, 11]. The study showed that, within the range of the parameters considered, the attenuation overwhelmed the enstrophy generation, which is in stark contrast to the enstrophy production in the macroscale shock-vortex interaction. However, critical questions of the microscale shock-vortex interaction, including whether enstrophy attenuation persists at higher incoming Mach numbers or a larger vortex, remain unsolved.
Since the DSMC has been proven to be successful in rarefied gas simulation [12–14], in the present work, it was selected to study the microscale shock-vortex interaction in which the vortex radial length is within the same order of the shock thickness. And, by extending previous flow cases treated by the DSMC [10, 11], interactions of shock waves up to Mach 3.5 with supersonic or subsonic vortex are considered in much greater detail.
2. Physical System and Simulation Method
2.1. Physical System
The physical system of the two-dimensional microscale shock-vortex interaction, with a moving vortex and a stationary shock, is shown in Figure 1. The microscale discrete vortex is formed by prescribing its initial flow to be that of a composite vortex [5, 11, 15]. The physical size of the composite vortex in area is about one-fourth that of a Taylor-type vortex with the same core size. And so, it is chosen for its computational advantages. The rotation center of the vortex is initially stationary and a velocity distribution between a core radius () and an outer radius () is prescribed for the vortex. In this physical system, the maximum tangential velocity is found at the core radius. Inside the core (), the velocity goes linearly to zero at . Outside (), the tangential velocity is set to zero. The tangential velocity distribution is set as follows:where
The temperature and pressure in the quiescent field surrounding the vortex are prescribed. Inside the vortex, pressure, density, and energy are determined by balancing the pressure gradients with the centripetal force, which is equivalent to solve the following system:
Rectangular grids are used for defining the computational domain. The upper and lower boundaries of the flow field are prescribed as periodic boundaries. And, the left and the right boundary are set to be the pressure boundary condition. The simple monatomic gas considered in the simulations is argon, and its initial quiescent state surrounding the vortex has a pressure of 1013 Pa and a temperature of 273 K. In this initial condition, the mean free path before the shock is m. In the present study, the composite vortex outer radius is double the size of the core radius, which means that .
2.2. DSMC Method
When the molecule’s mean free path is approximated to the characteristic length of the flow field, the hypothesis of continuum breaks down and the NSF equations are questionable. In this kind of problems, the gas flows are determined by Knudsen number (Kn), which is defined byHere is regarded as the characteristic length of the gas flows, and can be density, temperature, pressure, or velocity [16]. Based on the previous studies, the hypothesis of continuum breaks down when or even less [17].
Unlike traditional computational fluid dynamics (CFD) which is based on Navier-Stokes-Fourier (NSF) equations, the DSMC method, introduced by Bird, is a kind of numerical method to simulate flow field by solving the molecules’ moving and collision. And, it has become to be common numerical method in rarefied gas flows.
In DSMC method, the flow progress is divided into two independent steps: motions and collisions. And, time is also discretized so that motions and collisions processes can be executed independently. Also, the DSMC method reduces the number of actual gas particles to an exercisable level by employing simulation molecules to represent real gas particles. Every simulation molecule has position and velocity vector. After a time, the position of each simulation molecule is changed as the product of velocity vector and time. Collisions step simulated the collision process by considering the velocity vector and the position of simulation molecules. Then, macroscopic flow variables can be “sampled” through those simulation molecules. In DSMC method, only binary collision is considered. And, based on the conversation of mass, momentum, kinetic energy, and potential energy, several collision models have been proposed to adapt to different flow conditions. More details of DSMC method can be found in [18] and the program flowchart of a typical DSMC simulation is given in Figure 2.
The sample process can compute not only the macroscopic flow variables such as density, velocity, and temperature, but also the pressure tensor and viscous stress tensor. The local fluid velocity components and pressure are calculated through statistical method as follows:
The pressure tensor and viscous stress tensor can also be computed through statistical method:
Computational resource is a challenge for DSMC method. In the present work, a rectangular cell structure up to is used to simulate the flow field. However, more cells bring large statistical error which cause dramatically deviation. In order to reduce the statistical error and keep acceptable precision, the number the cell is cut to be . The amount of simulation particles in flow field is . Sampling process was performed in every ten steps. In the aspect of collision, we use the variable hard sphere (VHS) model and RSF sample method to consider the binary collision as well as the number of collisions, respectively.
3. Validation of the Methods
In order to validate the current DSMC models, we select two benchmark cases of microscale shock-vortex interaction investigated by Koffi et al. [10, 11]. Argon gas in its quiescent state surrounding the vortex at an initial temperature of 273 K and an initial pressure of 1013 Pa is considered. The core radius varies from to .
Unlike the interaction of macroscale shock-vortex interaction as shown in Figure 3, a substantial attenuation of enstrophy with time is observed as shown in Figure 4. It shows that the results of the validation in case 1 and case 2 in Table 1, which are the same as case 1 and case 4 in the present study as shown in Table 2, are qualitatively the same as Koffi’s DSMC results. And it is the most important contribution in the previous DSMC study by Koffi et al. and also it has been regenerated again in the present study. These validation studies confirm that the present DSMC method can be used for the study of microscale shock-vortex interaction.
4. Net Vorticity in Microscale Shock-Vortex Interaction
A substantial attenuation of enstrophy with time is also observed in the validation cases. These results are in stark contrast with those of the macroscale interaction that shows an increase in enstrophy when the vortex crosses a shockwave [19]. Quantitative studies of these features can be performed by investigating the area-weighted vortex dynamics. Seven representative cases are shown as follows.
4.1. Enstrophy and Dissipation Rate
The mechanisms leading to the generation or attenuation of vorticity in the interaction can be investigated by considering the time evolution of the enstrophy [11, 20], which is defined asAnd, the viscous effect is investigated by introducing the area-weighted dissipation rate of the kinetic energy:Here represents the dissipation rate per unit volume and is defined aswhere .
The time evolution of the enstrophy and the dissipation rate are plotted in Figures 5 and 6 for different shock Mach numbers of cases 1–3. Cases 2 and 3, with higher incoming Mach numbers, show quite different trend compared with the result of case 1 with lower incoming Mach numbers. An increase in the enstrophy is observed in cases 2 and 3 during the interaction process (100–600 ns). But, the enstrophy of case 1 appears to remain constant throughout the entire interaction process.
The difference is also evident in the time evolution of the dissipation rate, as illustrated in Figure 6. The dissipation rates of cases 2 and 3 with higher incoming Mach numbers are much greater than that of case 1. In addition, there are three peak values in the dissipation rate for cases 2 and 3 during the interaction process. But, the dissipation rates for case 1 appear to get less increment compared with those of case 2 and case 3 during the entire interaction process.
Figures 7 and 8 show the time evolution of the enstrophy and dissipation rate for cases 1, 4, and 5 with the same incoming Mach number and vortex Mach number but different vortex radius. The results of these three cases show different trend which exposes that vortex radius also plays significant role in the interaction. Case 4 and case 5 with a smaller core radius result in a decrease of enstrophy as shown in Figure 8. However, case 1 with a larger core radius remains constant during the interaction process.
Figure 8 shows that the time evolution dissipation rate and the increase are observed in all the three cases. The only difference is that case 1 with a bigger core radius causes a distinct increasing level during interaction process.
The effect of Mach number of vortex on shock-vortex interaction is shown in Figures 9 and 10. It shows that the vortex Mach number has a slight effect on the tendency of the enstrophy. Vortices, with different vortex Mach numbers, lead to a similar tendency in Figure 9, which is remaining constant through interaction process.
On the contrary, different vortex Mach numbers cause different trends in time evolution of dissipation rate. The dissipation rate of case 1 with greater vortex Mach number shows a strong increase. But, case 6 and case 7 with smaller vortex Mach number show a little increasing level and decreasing trend in interaction process.
Regarding time evolution of enstrophy, the primary cause of generation (increase) or attenuation (decrease) may be that of the viscous stress and its change due to high Knudsen numbers which dominate the flow structure during the interaction [21]. Strong shock wave, with high incoming Mach number, leads to a significant decrease of . But, at the same time, it has a little influence on . These cause the vorticity and enstrophy generation.
For dissipation rate, (9) can be written asBy coupling (7), (8), and (10), the alternative form of dissipation rate can be expressed as
In right hand side of (11), the first term describes the dissipation caused by vorticity and the second term describes the additional contribution of rotation to the dissipation rate. Then the third term represents the dissipation caused by the effect of compressibility [11]. This expression can be further used to analyze the trend of dissipation rate and its relationship with enstrophy.
According to (11), both the first term and the second term dominate the tendency of dissipation rate when the compressibility effect is similar. Cases 1, 4, and 5 have the same shock Mach number and vortex Mach number, which means they have equal compressibility effect. However, vortex with a small radius such as case 4 leads to a larger and which cause a greater amount of additional contribution of rotation (the second term of (11)). That is why the trend of dissipation rate still has a tiny generation while the enstrophy is decreasing during all the interaction process for case 4. Cases 1, 6, and 7 have the same shock Mach number and vortex radii, but different vortex number causes different compression effect. Lower vortex number leads to smaller compression effect; furthermore, it can also cause small and since cases 1, 6, and 7 have the same vortex radii. Both low compression effect and small additional contribution of rotation cause the attenuation in case 6 and then case 6’s tendency of enstrophy almost remains constant during interaction process. Cases 1, 2, and 3 have different shock Mach number; then strong compression effect makes the rise of dissipation rate and the higher shock Mach number is the greater dissipation rate generation, such as case 3. Besides, after the first peak value, additional two peaks of dissipation rate curve are observed for cases 2 and 3. Figures 11–18 present the pressure contour of case 3 during 380–440 ns. For convenience, high and low pressure zones are numbered as shown in Figure 18. As it is shown in Figure 17, two new high pressure zones of number 4 and number 5 are shaped into flow field. It is a new phenomenon in shock-vortex interaction with high shock Mach number. The reason of the formation of number 5 high pressure zone may be due to the mass accumulation caused by the strong mass transport before shock encountering the low velocity gas flow affected by vortex after the shock wave. However, number 4 high pressure zone has another form reason. Because the vortex is clockwise, two kinds of mass transport exist and point to number 4 high pressure zone. The first one is that the mass transport leads to vortex which has a great negative direction of the -axis velocity and small negative direction of the -axis velocity. The second one is the strong mass transport before the shock wave which has a big negative direction of the -axis velocity. Two kinds of mass transport come across in number 4 high pressure zone and make a notable mass accumulation. Because shock wave cannot have an effect on -axis velocity, number 4 high pressure zone develops towards the bottom of flow field.
In flow field, the intensity of high and low pressure zones dominates dissipation rate, and high or low pressure zones with less strength lead to a weak dissipation rate. However, the zones with greater strength cause a strong dissipation rate, which makes the trend of dissipation rate increase. The reason of the second peak in Figure 6 is the combined action of the dissipation of number 1 and number 2 low pressure areas and the formation as well as development of number 3, number 4, and number 5 high pressure areas. Although number 4 and number 5 high pressure areas have a strong intensity which can be seen in Figure 17, small area limits the increasing level of dissipation rate, which cause a tiny attenuation after time of 430 ns. The continual development of number 4 and number 5 high pressure region brings out the third peak of dissipation rate.
5. Conclusion
The present DSMC study successfully characterizes extensively microscale shock-vortex interaction cases, of a planar shock up to Mach 3.5 with a transverse composite microvortex on the shock thickness scale. New physical features are found in the microscale shock-vortex interaction. Strong interactions with high shock Mach number or large vortex radius cause significant increase in the enstrophy during the interaction, which results in an increase in the dissipation rate. It is not observed in the previous DSMC studies because of the limited studied cases. Weak interactions with low incoming or vortex Mach numbers cause a slight increase in the enstrophy during the interaction, which results in the decrease in dissipation rate throughout the entire process for weak interactions. It differs from previous macroscale shock-vortex interaction studies.
It is also observed that the incoming Mach number is the most dominant, followed by vortex size and vortex Mach number, during interaction. On the other hand, shock Mach number is the most important parameter during the whole interaction process.
Nomenclature
| Sound speed | |
| : | Time step |
| : | Vortex Mach number |
| : | Shock Mach number |
| : | The number of molecules in each cell |
| : | Pressure |
| : | Pressure tensor |
| : | Vortex radius |
| : | Vortex core radius |
| : | Vortex outer radius |
| : | Maximum tangential velocity |
| : | Molecule velocity |
| : | Molecule thermal velocity |
| : | Vortex tangential velocity |
| : | Density |
| : | Mean free path |
| : | Stress tensor |
| : | Kronecker tensor |
| : | Viscosity coefficient |
| : | The second viscosity coefficient. |
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The present research is supported by the Fundamental Research Funds for the Central Universities.