Solution to Propagation of Shock Wave Based on Continuum Hypothesis

Based on the continuum hypothesis, problem of compressing wave was studied analytically. By exploring the temperature distribution, the propagation velocity, and the thickness of the transition region, the developing process of compressing wave turning into shock wave was revealed, and the following conclusions were reached: (1) the governing equations of compressing wave can be turned into two autonomous equations, and the phase diagram can be used as an effective tool for the analysis of the compressing wave; (2) the solution of compressing wave was composed of two parts: the inner solution and the outer solution; when Pr = 3/4, the analytical inner solution can be obtained; (3) as the disturbance velocity increases, the thickness of compressing wave will decrease, and the propagation velocity of the compressing wave will increase, and the compressing wave will become shock wave when the disturbance velocity is large enough; (4) velocity and temperature across compressing wave change monotonically and continuously, but the entropy generation changing tendency is closely related to Pr; therefore, the inner solution reveals the mechanism of irreversibility happening in compressing waves.


Introductions
Compressing wave is an inherent phenomenon for compressible fluid flow.As all actual fluids have compressibility [1], compressing wave exists universally.According to the difference of intensity, compressing wave can be divided into sonic waves and shock waves; the latter ones, as irreversible processes, have great effects on a lot of equipment such as supersonic aircrafts and turbine blades and have been studied extensively.Currently, researches on the two kinds of compressing wave are carried out by different models; the classical solving method for weak shock wave is to assume it as an isentropic process and then get the solutions such as propagation velocity.But, for shock wave, the methods are different and more complicated.One simplified method is to assume it as a discontinuous plane and then get the relationship between parameters before and after shock wave, for example, the classical Prandtl formula and Rankine-Hugoniot relation formula [2], which are widely used in the research of oblique shock wave, reflection of shock wave, and interaction between shock waves and between shock wave and expansion wave.For example, Ng et al. [3] used this method when researching collusion of normal shock wave and detonation wave.It is convenient to use the discontinuous model to get the parameters outside shock wave.But if the variation process of the parameters within the shock wave was aimed to be revealed, the structure of the wave needs to be solved.As for this, two methods were proposed: one utilized the molecular dynamics to research the microstructure of the shock wave, and the other was based on the continuum hypothesis.Usually the thickness of shock wave is tiny, the magnitude of which is molecule mean free path.Thus, molecular dynamics model has great advantages when solving shock wave problems.It has been employed by many scholars.Malkov et al. [4] solved Boltzmann function of shock wave when Mach number is 4 and 8 based on hard sphere model.Bisi et al. [5] solved the structure of the solution for binary mixed shock wave based on molecular dynamics.Although the thickness of shock wave is extremely small, the second method is still a convenient and effective method, which typically uses the NS equations as description of fluid flow in the shock wave.Chen et al. [6][7][8]  of shock wave in rarefied gas and discussed the difference between that and molecular dynamics theory.Anand and Yadav [9] researched the structure for shock wave solution of magnetic fluid in nonideal gas; Elizarova et al. [10] and Elizarova et al. [11] analyzed the structure of shock wave solution in nitrogen and helium and other gases based on NS equation; Taniguchi et al. [12,13] solved the structure of shock wave solution of thin polyatomic gas based on NS Fourier equation and found that there were three different forms of shock wave solution under different Mach numbers.The significant advantage of the method based on the continuum hypothesis is that the analytical method can be applied, and some important results have been achieved.Morduchow and Libby [14] and Becker [15] got the accurate analysis solution of shock wave when Pr = 3/4, respectively.After that, a number of scholars had made deep researches on this foundation; Von Mises [16] simplified NS equations into two autonomous equations and studied them with phase diagrams; Hamad [17] had a further research on the influence of Pr on the entropy generation process in shock wave.It has been found that there is an extremum of entropy generation when Pr = 3/4, which is that along the shock wave; the entropy generation increases first and then decreases to the value obtained by Rankine-Hugoniot formula and Prandtl formula.In summary, solving compressing problem based on continuum hypothesis is still an important method, and it can also be used to explore the relationship and the transition process between sonic wave and shock wave.Now that the analytical solution of shock wave for Pr = 3/4 has been obtained and Prs of most actual gases are just the magnitude of 3/4, this article aims to reveal the principles for the process of the sonic wave turning to shock wave as the disturbance increases by means of analytical solution, and the parameters such as the entropy generation and velocity will be scrutinized.

Physical Model and Simplification of Governing Equations
2.1.Physical Model.The physical model established in this paper is shown in Figure 1.In a semi-infinite pipe, a piston moves to the right with a velocity  and compresses the gas in the right side; because of inertia, the gas cannot be disturbed immediately but is compressed gradually from left to right.In reality, the compressed region expands with a certain velocity , which is different from but closely related to the piston velocity ; the bigger  is, the faster the compressed region expands.The propagation of this disturbance is just the compressing wave.
Between the compressed region and the unaffected region, there exists a transition region, where parameters such as velocity, temperature, and pressure change continuously from the undisturbed state such as pressure p 1 and density  1 to disturbed state, p 2 and  2 .The existence of this transition region is the result of the momentum and energy exchange between the unaffected region and the disturbed region, and its thickness  is also related to the disturbance velocity V.When the disturbance velocity is small, the thickness of the region is large, and the compressing wave is kind of sonic wave; as the disturbance velocity increased, the compression became severe, and the compressing wave will become shock wave, and the thickness of the transition region will be much thinner; from the viewpoint of macro scope, the transition region seems to be a plane, across which the parameters change abruptly.
In order to construct the governing equations of the compressing wave, the conservation laws should be considered.First, as the piston moves and the compressing wave propagates, the gas is not leaked, so the mass conservation law must be satisfied; for the one-dimension flow, it can be written as below, which is called continuity equation in fluid dynamics.
Second, as the gas is compressed, the force is exerted between fluid elements, so the conservation law of momentum must be satisfied, which can be written as below, which is the one dimension form of the classical NS equations in fluid dynamics.In the equation, the compressibility of the fluid and the diffusion of the momentum were both considered.
Third, as the temperature of the compressed gas increases, the heat transfer between fluid elements exists, so the conservation law of energy must be considered, and its one dimension form can be written as below, in which the viscous dissipation item was considered and included in the last item.

Transformation of Governing Equation.
The governing equations ( 1)-(3) contain partial differential to time t, so the primary model illustrated in Figure 1 is an unsteady form.For the sake of convenience, the model can be turned into steady one by transformation of governing equation.If supposing the solution takes the form then a new coordinate  can be defined as In the equations above, velocity c is constant, and the partial differential terms appearing in ( 1)-( 3) can be transformed into ordinary differential ones: Then a new transformed velocity w can be defined as Equations ( 1)-( 3) can be transformed into ordinary differential equations: Thus, the unsteady problem discussed in Figure 1 can be transformed to a steady problem, in which the coordinate and the velocity are both relative ones, based on ( 5) and ( 7), respectively.Figure 2 shows the transformed steady model, which is more widely used for shock wave discussion.
In this model, compressible fluid flowed in a velocity of c initially, and after the compressing wave, it was compressed, and the velocity decreased to c − V, and the pressure and temperature increased accordingly.

Boundary Conditions.
To solve ( 8)- (10), two boundary conditions are needed for w and T, and one is needed for  and p.As for the velocity boundary condition, its values in the unaffected region and the affected region are both known, so it can be set as follows: As for the pressure boundary condition, its initial value in the unaffected region is known, so it can be set as follows: And so it is with the density boundary condition: As for the temperature, its initial value in the unaffected region is known, but its value in the compressed region is unknown.In order to construct two boundary conditions for the temperature, it is supposed that its gradient far from the transition region is 0, so the temperature boundary condition can be set as follows:

The Analytical Solution of the Compressing Wave
3.1.The Initial Integral of Governing Equations.The nonlinear equations composed of ( 8)∼ (10) can be simplified by one step of integral to reduce the order.First, the continuity equation ( 8) became an algebraic equation after integral as below, in which the integral constant m is the mass flow rate per unit area.
Assuming that the dynamic viscosity  is constant, the momentum equation ( 9) can be integrated as below, in which N is the integral constant which can be called the total pressure, as it has the dimension of pressure, and it equals the sum of pressure and momentum of the fluid when the gradient of relative velocity w is zero.Parameter  e is defined as 4/3 times the dynamic viscosity .
For the energy equation (10), assuming that the thermal conductivity is constant and gas is ideal gas, it can be integrated as below, in which the integral constant M can be called the total energy, as it has the dimension of energy, and it equals the sum of enthalpy and kinetic energy when the gradient of temperature is zero.

Autonomous Equations about T and w and Phase Diagram.
In order to discuss expediently, the ideal gas state equation is introduced, and p in momentum equation can be eliminated as follows: So the closed equations about T and w are established as ( 17)∼(18), and they are autonomous equations, which do not include independent variable , so the phase diagram is an effective method for analysis, which is a figure using the dependent variable as the coordinates and the direction field is drawn simultaneously.The direction field is a group of vectors at each point, and its direction and length are calculated by the autonomous equations.So, from the phase diagram, the tendency of the solution for each initial condition can be easily observed.For the problem discussed in this article, the coordinates are w and T, and in the phase diagram, each point corresponds to a combination of values of T and w, and the components of the direction vector at each point are calculated by ( 17)-(18), so the direction of the vector can be calculated as follows: The phase diagram for this problem discussed in this article is shown as Figure 3, and Figure 3(a) is for the case when Ma = 2.65, and Figure 3(b) is for the case when Ma = 1.22.The definition of Ma is as the equation below, the physical meaning of which is the ratio of relative velocity w to sonic velocity a: In Figure 3, there are two critical points, P and Q, which are two singular points.P is divergent, while Q is convergent, and directions of the direction vectors between these two points are all from P to Q.In fact, the two points correspond to the two steady states of the compressing wave, the unaffected state and the compressed state, and there must be a solution curve for the gas state changing from steady state P to steady state Q; the direction vectors on the solution curve are tangent to it.The critical points P and Q can be obtained by some special lines as shown in Figure 3, which are as follows.
L 1 is the iso-momentum line, which is defined as (21).Along this line, the sum of the momentum and static pressure is constant, and the equation corresponds to the case when the momentum equation ( 18) is equal to zero.
L 2 is defined from the case when ( 17) is equal to zero; therefore, the points on the line satisfy the following equation: L 3 is the iso-stagnation-enthalpy line, which is defined as (23) and is also deduced by supposing (18) and ( 17) to be zero.Along this line, the stagnation enthalpy is constant, which is the sum of enthalpy and kinetic energy.
L 4 is the sonic line, which is defined as (24).On this line, Ma equals 1. Region in the phase diagram above this line is subsonic and below is supersonic.
The critical points P and Q are the intersection of the three lines: L 1 , L 2 , and L 3 .The direction field is influenced by parameters such as adiabatic coefficient , inlet Mach number, and Prandtl number Pr.The direction field shown in Figures 3(a) and 3(b) corresponds to two cases with different Ma and other identical parameters, such as the adiabatic coefficient  being 1.4 and Pr being 3/4.As is shown, when Ma decreases, P and Q will approach the line L 4 , and in the limit, when Ma = 1, the three lines L 1 , L 2 , and L 3 will be tangent with each other, and the tangent point is on L 4 .

Analytical Solution for Pr = 3/4.
The solution of compressing wave is composed of two parts: the inner solution and the outer solution: the inner solution corresponds to the transition region; the gradient of it is large and the thickness is thin; the outer solution corresponds to the two steady states, which has gentle gradient and the dimension is larger, and the relationship between the two steady states can be obtained by traditional Rankine-Hugoniot relationship formula and Prandtl formula.The complete solution of a compressing wave must take into account the relationship between the two parts.However, for an arbitrary Pr, the analytical solution of the inner solution is hard to derive.But for the special case with Pr being 3/4, the analytical solution can be obtained.

Inner Solution.
As is shown in Figure 3, it can be found that the tangent of each point on the curve L 3 coincides with its direction vector.Since L 3 itself passes through the two steady-state points of P and Q, the curve L 3 coincides with the T-w solution curve of the compressing wave, and ( 23) is a form of its analytical solution.According to this equation, T can be eliminated and solution about w can be first obtained as follows: This equation is an implicit solution, in which w 1 and w 2 denote the velocity of the two steady states, and the independent variable  is expressed as a function of the dependent variable w.As for the temperature T, its solution can be obtained according to (23) based on the velocity solution (25), as is shown below:

Outer Solution.
In the inner solution, the integral constants are still unknown, such as m, N, M, w 1 , and w 2 ; to get the complete solution of a compressing wave, the outer solution is needed.In this regard, it is necessary to first determine the propagation velocity of the compressing wave c.Due to the external solution based on the fact that the continuum hypothesis is consistent with the shock relation obtained from the discontinuity plane model, the classical Prandtl formula is still applicable to the outer solution, and the velocity of the compressing wave c and piston velocity V satisfy the following equation: The velocity of the compressing wave can be obtained by solving the above algebraic equation as below, in which a stands for the sonic velocity.
Based on this velocity, other integral constants can be determined: the total energy M before and after compressing wave can be determined according to the boundary conditions ( 11)-( 14), and the result is as shown in the following equation: The mass flow rate per unit area m can be determined as follows: And the total pressure N can be determined as follows: Figure 4 shows the distribution curves of relative velocity and temperature for two specific cases.Figure 4(a) shows the distribution of relative velocity w when the piston velocity is 5 m/s and 10 m/s. Figure 4(b) shows the distribution of temperature T when the piston velocity is 5 m/s and 10 m/s.For the two cases, the initial temperature of the gas is 300 K, and the initial velocity is 350 m/s.In the two figures, the relative coordinate  is drawn in accordance with the transformed model shown in Figure 2, with the positive direction pointing to the left.It can be seen that when the piston velocity is larger, the velocity of the compressed region is smaller, and the temperature is higher.
In Figure 4, the thickness of the transition region is also illustrated, the definition of which is the distance of the two locations at which the velocity difference is 0.99 times the steady velocity difference, as is illustrated in the following equation: For the case of piston velocity being 10 m/s, the thickness  1 is 0.026 mm, while for the case of piston velocity being 5 m/s, the thickness  2 is 0.052 mm.In fact, as the piston velocity increases, the thickness will decrease drastically.26) are only applicable for the case of Pr = 3/4; for other Prandtl numbers, the analytical solution will be difficult to derive.So, the direction field can be very useful for analysis, and the internal solution can be obtained by numerical calculation.Figure 5 shows the direction field in the T-w diagram for Pr = 10 and Figure 6 is for Pr = 0.3.It can be found that when Pr > 3/4, the T-w solution line will shift to the line L 2 ; and when Pr < 3/4, the line shifts to L 1 .And mathematically, when thermal conductivity  tends to be zero, Pr will tend to be infinity; thus the T − w curve will be infinitely close to L 2 .Similarly, when  tends to be zero, Pr will tend to be zero, and the T-w curve will be infinitely close to L 1 .Therefore, in compressing wave process, heat conduction and viscosity cannot be neglected; it is just their existence that leads to the entropy generation across the compressing wave.

Thickness of Compressing Wave and the Propagation
Velocity.As the piston velocity increases, the thickness of the transition region  will decrease and the propagation velocity of the compressing wave c will increase, as is shown in Figure 7.In Figure 7, the abscissa represents the piston velocity, which is drawn in the form of logarithmic coordinate in order to show the variation tendency of  and c at low piston velocity.The left ordinate represents the thickness of the transition region, which is also drawn in logarithmic form.
The parameters used for calculation in Figure 7 are as in Table 1.
It can be found that when the piston velocity V is large enough, the thickness of compressing wave will decrease drastically to the magnitude of the 0.1 m, and the propagation velocity will increase to the magnitude of 10 3 m/s, which is far larger than the sonic velocity a, and the compressing wave will become a shock wave.
It is noticeable that the above results are obtained from the viewpoint of mathematics, while for actual situations the applicability of continuum hypothesis needs to be considered.For a small disturbance velocity V, Ma calculated by c will not be very large and the thickness of compressing wave will be much larger than the molecule mean free path, which is of the magnitude of 10 −6 m, and, in this case, the continuous solution of the compressing wave is reasonable.As the disturbance velocity V increases, Ma will increase and the thickness of compressing wave will gradually reach the order of the molecule mean free path.For example, when  = 280 m/s, Ma = 1.6; then  will be 10 times of the molecule mean free path; the applicability of the continuum hypothesis needs to be further discussed.

Relationship between Entropy Generation and Temperature.
From the analytical solution, the entropy generation process in the compressing wave can be calculated easily as (33) which is taken as the dimensionless form based on the universal gas constant R g , and T * and w * are the dimensionless temperature and relative velocity, respectively, which are based on the parameters T 1 and w 1 in the unaffected sated, as is illustrated in (34).
The relationship between the entropy generation and the temperature across compressing wave for the case of Pr = 3/4 under different Ma is shown in Figure 8.In the figure, three lines of entropy generation and temperature relationship were drawn, and the fourth line illustrated the relationship between the net entropy generation and temperature across compressing waves, which can be obtained by the Rankine-Hugoniot relation formula.It can be observed that the entropy generation does not increase monotonously with the increase of temperature but increases first and then decreases to the definite value calculated by discontinuity plane model.For given  and Ma, the net entropy generation is definite, but different Pr corresponds to different curves of Δs-T.If the Pr is large enough, it can be monotonous, while for small Pr, the nonmonotony of Δs-T relationship is severe.

Conclusions
Based on the continuum assumptions, distribution of the parameters such as temperature, velocity, and entropy generation across the compressing wave was calculated analytically and the following conclusions were reached.
(1) As the governing equations of compressing wave can be turned into two autonomous equations, the phase diagram is an effective tool for the analysis of the compressing waves, and by means of it, the special case when Pr = 3/4 was discovered, in which the iso-stagnation-enthalpy characteristic line was just tangent to a set of vectors, which is just the relationship curve between the parameters of temperature and velocity, so the analytical solution can be obtained.As for other cases, the vector field in the phase diagram can reflect the guideline of the solution.
(2) The solution of compressing wave is composed of two parts: the inner solution and the outer solution.And the inner solution is greatly influenced by Pr; when Pr = 3/4, the analytical inner solution can be obtained, in which the specific stagnation enthalpy is constant, which is the sum of the kinetic energy and enthalpy per unit mass.The solution of velocity has the form of implicit function, and the solution of temperature can be solved based on velocity solution.
(3) As the disturbance velocity increases, the thickness of compressing wave will decrease, and the propagation velocity of the compressing wave will increase.When the disturbance velocity is large enough, the thickness of the transition region will reach the magnitude of molecule mean free path, and the propagation velocity will be much larger than sonic velocity, and the compressing wave will become shock wave.

Figure 1 :
Figure 1: Propagation of a one-dimension compressing wave driven by a moving piston.

Figure 2 :
Figure 2: Transformed compressing wave model in the new coordinate.

Figure 4 :Figure 5 :
Figure 4: Distribution of relative velocity  and temperature  in the shock wave.

Figure 6 :
Figure 6: Direction field in the phase diagram for Pr = 0.3.

Figure 7 :
Figure 7: Thickness of transition region and propagation velocity of compressing wave for different piston velocity.

Figure 8 :
Figure 8: Relationship of entropy generation and temperature across compressing wave for different Ma.

Table 1 :
Calculating parameters for case shown in Figure7.