Hydrodynamic and Thermodynamic Nonequilibrium Effects around Shock Waves: Based on a Discrete Boltzmann Method

A shock wave that is characterized by sharp physical gradients always draws the medium out of equilibrium. In this work, both hydrodynamic and thermodynamic nonequilibrium effects around the shock wave are investigated using a discrete Boltzmann model. Via Chapman–Enskog analysis, the local equilibrium and nonequilibrium velocity distribution functions in one-, two-, and three-dimensional velocity space are recovered across the shock wave. Besides, the absolute and relative deviation degrees are defined in order to describe the departure of the fluid system from the equilibrium state. The local and global nonequilibrium effects, nonorganized energy, and nonorganized energy flux are also investigated. Moreover, the impacts of the relaxation frequency, Mach number, thermal conductivity, viscosity, and the specific heat ratio on the nonequilibrium behaviours around shock waves are studied. This work is helpful for a deeper understanding of the fine structures of shock wave and nonequilibrium statistical mechanics.


Introduction
Shock waves are widespread in many fields of physics and engineering [1][2][3][4], and their frequent occurrence is based on the fact that matter is more or less compressible [5]. The shock wave is a kind of disturbance waves that carries abundant energy and propagates at supersonic speed in the medium [6]. At the macroscopic level, it is a high-temperature, high-pressure, high-density surface moving forward at high-speed, through which the pressure, density, temperature, and other physical quantities of the medium undergo rapid changes [7,8]. Shock waves, as an objective phenomenon, are generated by lightning, earthquakes, volcanic eruptions, nuclear or chemical explosions, reentry vehicles, sonic boom of supersonic aircraft and any supersonic flying projectile, etc. [5,8].
In practice, the shock wave is extremely complex, as it covers a wide range of spatio-temporal scales, and incorporates hydrodynamic and thermodynamic nonequilibrium effects [9,10]. In fact, hydrodynamic and thermodynamic nonequilibrium effects often exert significant influences on fluid systems [11]. It is necessary to take these nonequilibrium effects into consideration in order to accurately predict fluid behaviours. To this aim, a possible method is to use molecular dynamics [12][13][14] or direct simulation Monte Carlo [15]. But both of them encounter the issue that the spatio-temporal scales amenable are rather limited due to the high computation cost. On the other side, traditional hydrodynamic methods which require less computation cost usually ignore more abundant, complex and essential thermodynamic nonequilibrium effects caused by interactions at microscopic scales [16,17]. The hydrodynamic models are based on the continuity assumption and usually lack the high physical accuracy for practice flows with sharp physical gradients and/or strong nonequilibrium effects. To address this problem, great efforts have been made and various kinetic models based on the Boltzmann equation have been developed.
Actually, there are two classes of kinetic models for fluid systems. One aims to provide a new algorithm to solving the hydrodynamic governing equations, e.g., Euler equations, incompressible or full Navier-Stokes equations. Such models include the lattice Boltzmann method [18][19][20][21] and the gas kinetic scheme [22,23]. These models could only simulate the fluid behaviours described by the original macroscopic fluid equations, they can not provide any information beyond the original equations. The other type of kinetic models is to capture both hydrodynamic and thermodynamic nonequilibrium behaviours beyond the macroscopic models. The successful physical models include the unified gas kinetic scheme [24][25][26], the discrete unified gas kinetic scheme [27,28], the discrete Boltzmann method (DBM) [29][30][31][32][33][34][35][36][37][38][39], etc. These powerful models are suitable for continuum and rarefied systems with a wide range of Knudsen numbers, and capable to subsonic and supersonic flows with essential nonequilibrium effects.
Compared with the unified gas kinetic scheme and the discrete unified gas kinetic scheme, the DBM requires less discrete velocities and distribution functions, hence it has a higher computational efficiency. From a physical modeling perspective, a DBM is approximately equivalent to a continuous fluid model plus a coarse-grained model of other relevant thermodynamic nonequilibrium effects. At present, the DBM has been developed as a nonequilibrium flow simulation tool for various flow systems, including the fluid instability [29,32,38], multiphase flow [34], shock and detonation [30,31,33,36,37,39]. Roughly speaking, the DBMs can be classed into two categories: the single-relaxation-time DBM [29][30][31][32][33][34][35][36] and the multiple-relaxation-time DBM [37][38][39]. In the single-relaxation-time model, there is only one relaxation time that controls a thermodynamic nonequilibrium system approaching its equilibrium state. In the multiple-relaxation-time DBM, the relaxation parameters for various kinetic modes are independent of each other. Actually, the multiple-relaxation-time DBM reduces to a single-relaxation-time model if all relaxation parameters equal to each other.
In this article, we use the multiple-relaxation-time DBM to study nonequilibrium shock waves, where both specific heat ratio and Prandtl number are adjustable [39]. Thanks to its kinetic nature, the versatile DBM can be adopted to capture and measure both hydrodynamic and thermodynamic nonequilibrium behaviours in an accurate and effective way [39]. The DBM is briefly introduced in Appendix A. The rest of the paper is structured as follows. The study of local nonequilibrium effects is presented in Section 2. Then the investigation of global nonequilibrium effects is shown in Section 3. Finally, Section 4 concludes.

Local Nonequilibrium Effects
It is noteworthy that the DBM has the capability of measuring the nonequilibrium manifestations in subsonic, sonic and supersonic flows accurately [29][30][31][32][33][34][35][36][37][38][39]. Actually, the DBM is capable to recover the velocity distribution function in two ways. One is to obtain the main features of the velocity distribution function by analysis of the detailed nonequilibrium quantityf neq i that has specific physical meanings. For example,f =f 9 −f eq 9 is twice the nonorganized energy flux in the y direction. In previous works [29,30,40], the DBM has already been utilized to get the characteristics of the velocity distribution function around the detonation wave, shock front, rarefaction wave, and material interface. The other is to recover the main characteristics of the velocity distribution quantitatively through the macroscopic quantities and their spatial and temporal derivatives, which can be achieved by the Chapman-Enskog expansion [35].

Nonequilibrium Manifestations
To have a better understanding of the nonequilibrium effects upon the fine structure of the shock wave, we probe the nonequilibrium variables, effects, and degrees in this subsection. Note that the main characteristics of the velocity distribution function can be obtained from the thermodynamic nonequilibrium manifestations [29,30,40]. Now we consider a shock wave propagating from left to right. At the beginning, the shock front is located at x = 0.01 with the Mach number Ma = 2. The computation is carried out with the mesh number N x × N y = 10,000 × 1. The spatial step is ∆x = ∆y = 10 −4 , the temporal step ∆t = 10 −5 , and relaxation frequency S i = 10 3 . The initial state is set by the Hugoniot relation, ρ, u x , u y , T, P L = (2.667, 1.479, 0, 1.688, 4.5) , ρ, u x , u y , T, P R = (1, 0, 0, 1, 1) , where the subscripts L and R denote 0 < x ≤ 0.01 and 0.01 < x ≤ 1, respectively. In addition, inflow and outflow boundary conditions are employed in the x direction, and the periodic boundary condition is adopted in the y direction.
The profiles of density ρ, pressure P, temperature T, and horizontal velocity u x near the shock front are illustrated in Figure 1a. The solid lines with squares, circles, upper triangles and lower triangles represent the density, pressure, temperature, and horizontal velocity, respectively. Three vertical dashed lines are plotted to guide the locations at x 1 = 0.9545, x 2 = 0.9575 and x 3 = 0.9605, respectively. It's clear that these quantities increase dramatically as the shock wave travels forward through the medium.   are related to the gradients of fluid velocity and temperature (∇u and ∇T), respectively [41]. As shown in Figure 1b, the troughs of the physical gradients ∇u x and ∇T are close to each other. This is the reason why the locations of their minima are near. Figure 1d shows the profiles of nonequilibrium effect (Φ) and deviation degrees (∆ r and ∆ a ). The line with squares denotes the nonequilibrium effect, (1) The solid lines with upper triangles and circles display the absolute and relative deviation degrees defined as, and respectively. Obviously, ∆ a , ∆ r and Φ first grow, then decrease and form peaks near the wavefront. The peak of Φ nearly coincides with those off depend upon the gradients of fluid velocity and temperature. The peak of ∆ r is closer to the peak of Φ, and the peak of ∆ a is located behind them. The hysteresis of this peak arises from the change in density as it passes through the shock wave.

Recovery of Velocity Distribution Function
Next, we derive the formula of the velocity distribution function through the macroscopic quantities and their spatial and temporal derivatives for the multiple-relaxation-time DBM. Compared with the derivation in Ref. [35], the specific heat ratio is under consideration here. Then the equilibrium and nonequilibrium velocity distribution functions are recovered and analyzed at three different locations around the shock front that propagates forward.
Let us consider the case with the Prandtl number Pr = 1, which means the relaxation frequency S i equals to each other in the multiple-relaxation-time discrete Boltzmann model. Additionally, the relaxation time is defined as τ i = 1/S i . Via the Chapman-Enskog analysis, the first order approximation of the velocity distribution function f can be obtained, as below with the equilibrium distribution function stands for the space dimension, I indicates extra degrees of freedom due to vibration and/or rotation, and η corresponds to vibrational and/or rotational energies. Appendix B provides the detailed derivation process.   intuitively show the 3D distribution functions in velocity space (v x , v y ) at different points. Figure 2(a2-c2) plot the two-dimensional (2D) contours of distribution functions. Figure 2(a3-c3) are the 1D distribution functions, where the solid and dashed lines represent the velocity distribution function f (v x ) = f dv y dη and its equilibrium counterpart f eq (v x )= f eq dv y dη, respectively. Figure 2(a4-c4) display the other 1D profiles of the distribution functions, where the solid lines denote f v y = f dv x dη and the dashed lines indicate f eq (v y )= f eq dv x dη. From Figure 2, we can obtain the following points.
(I) There are clear peaks in the Figure 2(a1-c1). From left to right, the peaks get "smaller" and more "pointy". Because the shock wave passes through the medium, the physical quantities increase rapidly. The velocity distribution functions become progressively sharper from post-wave to pre-wave.
(II) The contours presented in Figure 2(a2-c2) are symmetrical in the v y direction and asymmetric in the v x direction at various locations. From left to right, the contours gradually look like eggs due to the nonequilibrium effects. The reason is that the system is periodic in the y direction and there is flow or flux just in the x direction.
(III) As Figure 2(a3-c3) show, f eq (v x ) is symmetric and the symmetry axis is u x , while f (v x ) is asymmetric due to the nonequilibrium effects. The velocity distribution function f (v x ) is "lower" and "fatter" than the equilibrium velocity distribution function f eq (v x ), which indicates the nonorganized energy f neq 5 > 0 (see Figure 1). As v x deviates far from u x , the distribution function Figure 1).
(IV) Both f (v y ) and f eq (v y ) in Figure 2(a4-c4) are symmetrical and the axes of symmetry are u y = 0, because there is no flow or flux in the y direction. Additionally, the velocity distribution function f (v y ) "higher" and "thinner" than the equilibrium velocity distribution function f eq (v y ).
, and f eq (v y ) decrease gradually from left to right. In fact, the area corresponds to density ρ, which increases as the shock wave travels from left to right. The widths of f (v x ), f (v y ), f eq (v x ), and f eq (v y ) are related to the temperature T, which increases from the first to the third columns. As a result, the curves become "higher" and "leaner".
It should be mentioned that the numerical results presented in Figure 1 are consistent with the analytic solutions in Figure 2. Up to this point, from Figures 1 and 2, we can obtain an intuitive understanding of the deviation between the nonequilibrium and equilibrium velocity distribution functions around the shock wave. More simulations are conducted in the following section in order to have a better understanding of the physical mechanisms of the nonequilibrium shock waves.

Global Nonequilibrium Effects
We define the following nonequilibrium quantities to further study the nonequilibrium effects around the shock wave, which denote the global nonequilibrium effects, the global nonorganized energy in the x direction, and the global nonorganized energy flux in the x direction, respectively. Next, the DBM is employed to study the effects of the relaxation frequency, Mach number, thermal conductivity, viscosity, and specific heat ratio on those global nonequilibrium variables around shock waves, respectively.

Impact of the Relaxation Frequency
In the physical sciences, relaxation usually means the return of a perturbed system into equilibrium. In our simulations, the relaxation frequency S i governs the speed of relaxation process from a nonequilibrium state to its equilibrium state. It plays a key role in the nonequilibrium manifestations. Here, let us investigate the nonequilibrium effects around the shock wave under different relaxation frequencies. The relaxation frequencies are chosen as S i = 1 × 10 3 , 2 × 10 3 , 4 × 10 3 , 8 × 10 3 , 1.6 × 10 4 , and 3.2 × 10 4 , respectively.
First of all, to give comparisons between Riemann solutions and DBM results, the profiles of density ρ, horizontal velocity u x , pressure P, and temperature T near the shock front are shown in Figure 3. The black solid line represents the result of the Riemann solutions, and the other lines stand for the DBM results under various relaxation frequencies. It is evident that the Riemann solutions show a straight vertical line near the shock wave, while the DBM results display a smooth curve around the shock front. Because the DBM contains viscosity and thermal conductivity as well as other detailed thermodynamic nonequilibrium effects, while the Riemann solutions are based on the Euler equations that ignore all thermodynamic nonequilibrium effects. As the relaxation frequency increases, the nonequilibrium effects decrease, and the DBM results are close to the Riemann solutions.    Figure 4b-d illustrate the global nonequilibrium quantities ln Ψ, ln Ψ 2,xx , and ln Ψ 3,1,x versus ln S i . Circles, triangles, and squares denote the simulation results, and the lines represent the fitting equations: ln Ψ = 2.892 − 2 × ln S i , ln Ψ 2,xx = −0.423 − 2 × ln S i and ln Ψ 3,1,x = 1.008 − 2 × ln S i , respectively. It is evident that the fitting results agree well with the simulation solutions, and ln Ψ, ln Ψ 2,xx and ln Ψ 3,1,x decrease as ln S i increases linearly. Physically, with the increasing relaxation frequency, the relaxation time fromf i tof eq i decreases, and the nonequilibrium effects around the shock wave reduce [31].

Impact of the Mach Number
Mach number Ma is a dimensionless quantity in fluid dynamics that represents the ratio of flow velocity past a boundary to the local speed of sound [42]. Now, let us demonstrate the relationship between the nonequilibrium effects and Mach number. We simulate a shock wave with various Mach numbers from Ma = 1.5 to 5 at time t = 0.4, with Pr = 1. The Hugoniot relation sets the initial states in all cases. Figure 5 exhibits the nonequilibrium effects versus the Mach numbers.

Impact of the Specific Heat Ratio
In thermal physics and thermodynamics, the specific heat ratio is the ratio of the heat capacity at constant pressure to heat capacity at constant volume. In this subsection, the simulation of shock wave is carried out for different specific heat ratios. The specific heat ratio is γ=(D + I + 2)/(D + I).
In our simulations, the specific heat ratio γ is adjusted by varying the extra degrees of freedom I. The extra degrees of freedom are I = 7, 6, 5, 4, 3, 2, 1, and 0, and the corresponding specific heat ratios are γ = 1.22, 1.25, 1.29, 1.33, 1.4, 1.5, 1.67, and 2, respectively. is related to the gradient of temperature and thermal conductivity. Consequently, with the increases of extra degrees of freedom, the specific heat ratio decreases, the temperature and its gradient decrease, and the thermal conductivity increases. The decreasing temperature gradient and increasing thermal conductivity have opposite impacts on the nonorganized energy flux, and temperature gradient plays a dominant role. So the increasing specific heat ratio result in an increase of Ψ 3,1,x [41].

Conclusions
Shock wave is relevant to many fields of physics and engineering. When the shock wave passes through the medium, it draws the system out of equilibrium. In this article, we use the multiple-relaxation-time DBM in order to study hydrodynamic and thermodynamic nonequilibrium effects around the shock wave. Via the Chapman-Enskog analysis, we derive the formula of the velocity distribution function through the macroscopic quantities and their spatial and temporal derivatives with the flexible specific heat ratio. Moreover, we probe the nonequilibrium variables and effects in order to have a better understanding of the nonequilibrium effects upon the fine structure of the shock wave. In addition, to give an intuitive description, we define the absolute deviation degree ∆ a and relative deviation degree ∆ r in order to describe the departure of velocity distribution function from its equilibrium counterpart. It is found that the deviation degrees and nonequilibrium effects first grow, then decrease, and form peaks at the shock wave. Additionally, there are small or remarkable distances among those peaks, because the nonequilibrium behaviours are associated with the physical gradients whose peaks are near or far from each other around the shock wave.
We illustrate 1D, 2D, and 3D velocity distribution functions at three different locations around the shock front that propagates forward in the x direction in order to quantitatively study the characteristics of the distribution functions. The 3D velocity distribution functions become progressively sharper from post-wave to pre-wave. The 2D velocity distribution functions near the shock wave are symmetrical in the v y direction and asymmetrical in the v x direction. The contours look like eggs gradually from left to right due to the nonequilibrium effects. As for the 1D velocity distribution functions, the areas of f (v x ), f (v y ), f eq (v x ), and f eq (v y ) correspond to density ρ, and they gradually increase as the shock wave travels from left to right. Meanwhile, the widths of f (v x ), f (v y ), f eq (v x ), and f eq (v y ) that are related to the temperature T increases step by step when the shock wave propagates forward. As a result, the 1D velocity distribution function curves become "higher" and "leaner" from post-wave to pre-wave.
Besides, comparisons are made between the Riemann solutions and DBM results. Physically, the DBM contains viscosity and thermal conductivity, as well as other detailed thermodynamic nonequilibrium effects, while the Riemann solutions are based on the Euler equations, which ignore all of the thermodynamic nonequilibrium effects. As the relaxation frequency increases, the nonequilibrium effects decrease, and the DBM results become close to the Riemann solutions. Moreover, we study the impacts of relaxation frequency, Mach number, thermal conductivity, viscosity, and specific heat ratio upon the nonequilibrium effects. The nonequilibrium quantities under consideration include the global nonequilibrium effects Ψ, the global nonorganized energy in the x direction Ψ 2,xx and the global nonorganized energy flux in the x direction Ψ 3,1,x . From the numerical simulations and corresponding fitting equations, the following points can be obtained: (I) With the increasing relaxation frequency, the amplitude of the nonequilibrium effect remains constant, while the width of the nonequilibrium region increases proportionally. Additionally, ln Ψ, ln Ψ 2,xx and ln Ψ 3,1,x decrease as ln S i increases linearly. (II) As the Mach number increases, all of the nonequilibrium effects increase, the global nonequilibrium effects Ψ increases in a linear relationship approximately, while the global nonorganized energy Ψ 2,xx and global nonorganized energy flux in the x direction Ψ 3,1,x grow exponentially. (III) With the decrease of the thermal conductivity, the global nonequilibrium effects Ψ gradually decreases and tends to be a constant, and Ψ 3,1,x grows linearly with the increasing κ. (IV) With the decrease of viscosity, the global nonequilibrium effects Ψ decreases gradually and then approaches a constant, and there is a linear growth relationship between Ψ 2,xx and µ. (V) As the specific heat ratio increases, the global nonorganized energy Ψ 2,xx decreases in a linear ship, while the nonorganized energy flux Ψ 3,1,x increases in an exponential form. The physical mechanisms are analyzed and explained for the above phenomena. This work is helpful for obtaining a deeper understanding of the fine structures of shock waves and the nonequilibrium statistical mechanics.

Appendix A
In this work, the governing equation of the DBM is the multiple-relaxation-time discrete Boltzmann equation as follows [39], where t is the time, r α is the physical space, v i denotes the discrete velocity with i = 1, 2, · · · , N, and N = 16 is the total number of discrete velocities, f i and f eq i represent the discrete distribution function and its equilibrium counterpart in velocity space, respectively.f i andf eq i represent discrete distribution function and its equilibrium counterpart in kinetic moment space, respectively. The parameter S lk is the element of the matrix S = diag [S 1 , · · · , S k , · · · , S N ] and controls the relaxation speed off i approachingf eq i . Besides, we choose the two-dimensional-sixteen-velocity (D2V16) model, with η i = η α , η b , η c , and η d for 1 ≤ i ≤ 4, 5 ≤ i ≤ 8, 9 ≤ i ≤ 12, and 13 ≤ i ≤ 16, respectively. Figure A1 delineates the sketch of D2V16 . In addition, the equilibrium distribution function f eq can be expressed by ρ, u α and T, the derivative of f eq can turn into the derivative of ρ, u α and T.