The Entropy Density Behavior across a Plane Shock Wave

Abstract

Entropy density behavior poses many problems when we study non-equilibrium situations. In particular, the local equilibrium hypothesis (LEH) has played a very important role and is taken for granted in non-equilibrium problems, no matter how extreme they are. In this paper we would like to calculate the Boltzmann entropy balance equation for a plane shock wave and show its performance for Grad’s 13-moment approximation and the Navier–Stokes–Fourier equations. In fact, we calculate the correction for the LEH in Grad’s case and discuss its properties.

1. Introduction

Although the entropy concept has been considered for about 150 years, there are issues that have pervaded until the first quarter of the XXI century; a crucial one was presented succinctly by Hoover and Hoover in the preface to the first edition of their book entitled “Time Reversibility, Computer Simulation, Algorithms, Chaos” [1]:

Today a small army of physicists, chemists, mathematicians and engineers has joined forces for a renewed attack on a classical problem, the “irreversibility paradox”.

The irreversibility paradox is also known as Loschmidt’s paradox, since he pointed out that Newton’s laws are invariant under time reversal while, Boltzmann’s $\mathcal{H}$ function is not. Later on, Zermelo used Poincaré’s recurrence theorem to provide further criticisms of Boltzmann’s probabilistic interpretation of the second law, which led to a continuing exchange between Boltzmann and Zermelo [2,3]. The more recent discoveries (chaos, Lyapunov exponents, thermostats and fractals) have played a role in discussions of the irreversibility paradox in the XX and XXI centuries [1].

There is a comment by von Neumann to Shannon in relation to a question raised by the latter about what would be an appropriate term for a formula that he derived; the following answer by von Neumann to Shanon was communicated to Myron Tribus by Shannon [4]:

“You should call it entropy, for two reasons. In the first place your uncertainty function has been used in statistical mechanics under that name, so it already has a name. In the second place, and more important, no one really knows what entropy really is, so in a debate you will always have the advantage.”

The previous joke is important because it shows that entropy is not a simple concept. In fact, its meaning must be given within a clear context which settles the frame for its correct interpretation.

The concept of entropy may be attributed to Claussius (1879) [5], however, as it is usual for scientists to use the previous intuition of other authors. In the case of entropy, Carnot’s ideas [6] were essential to the development of the entropy concept, providing Clausius with relevant insights about it [7,8]. The origin of the concept by Carnot and Clausius was in the phenomenological field of classical thermodynamics, meaning that it is given in terms of macroscopic variables. In the second half of the XIX century, Boltzmann [9] adopted the nascent atomic hypothesis, an idea that at such time had detractors, to develop an expression for entropy which is engraved on his tombstone [7]

$${\mathcal{S}}_{{}_W}=k_{B}lnW,$$

(1)

where ${\mathcal{S}}_{{}_W}$ is the entropy, $k_B$ is the Boltzmann constant and W the number of microstates compatible with the conditions of an isolated macroscopic system. Furthermore, he constructed the Boltzmann kinetic equation and defined his famous $(\mathcal{H})$ functional in terms of the distribution function $f(\mathbf{c},t)$, for homogeneous systems. It is simply generalized to consider the non-homogeneous case in which $f(\mathbf{c},\mathbf{x},t)$, then

$$\mathcal{H}(\mathbf{x},t)=\int f(\mathbf{c},\mathbf{x},t)\mathrm{L}nf(\mathbf{c},\mathbf{x},t)d\mathbf{c},$$

(2)

where $f(\mathbf{c},\mathbf{x},t)d\mathbf{c}d\mathbf{x}$ gives us the number of particles which at time t are in the volume element $d\mathbf{c}d\mathbf{x}$. It should be mentioned that the distribution function is now defined in the so called $\mu$ space $(\mathbf{c},\mathbf{x})$ and it must be a solution of the Boltzmann equation. In the case of equilibrium, the distribution function becomes Maxwell’s velocity distribution [3,10], and it gives the negative of the equilibrium entropy for an ideal gas, and as a consequence, it is natural to define the local entropy as

$${\mathcal{S}}_{\mathcal{H}}(\mathbf{x},t)=-k_B\mathcal{H}(\mathbf{x},t).$$

(3)

The entropy definition given in Equation (3) is valid for non-equilibrium situations as far as the $\mathcal{H}$ functional is well defined.

We do not pretend to give a review on the origin of the concept and its development; the interested reader will find critical insights in references [1,3,4,9,10,11,12,13] among others.

In this work our objective is more modest and amounts to a discussion of the entropy as given by Boltzmann’s $\mathcal{H}$ entropy $({\mathcal{S}}_{\mathcal{H}})$ for a shock wave and its relation with the local equilibrium hypothesis. As we will see, even in this simple case some problems arise.

The structure of this work is as follows. After this section we discuss the shock-wave problem in Section 2. Section 3 deals with the distribution function for Grad’s 13-moment approximation and the one provided by the Chapman–Enskog method at the Navier–Stokes–Fourier (NSF) level. In Section 4 we deal with the entropy density, and in Section 5 we provide numerical results for several profiles. We end by giving some final remarks in Section 6.

2. The Shock-Wave Structure

Let us consider a dilute gas in which a strong perturbation is produced so that a shock wave propagates with a constant velocity along the positive x direction. In the reference frame traveling with the shock wave we will have a stationary problem. Under these conditions, the usual mass, momentum and total energy conservation equations reduce to the case where the mass, linear momentum and total energy fluxes become constants. Then

$$\begin{array}{c}\hfill \rho u=c_1,\rho u^2+P_{xx}=c_2,\end{array}$$

(4)

$$\begin{array}{c}\hfill \rho \left(\frac{u^2}{2}+e_{int}\right)u+P_{xx}u+q_x=c_3,\end{array}$$

(5)

where $\rho$ is the mass density, u the velocity, $P_{xx}=p+p_{xx}$ is the $xx$ component of the pressure tensor, p is the pressure, T the temperature, $p_{xx}$ the $xx$ component of the viscous tensor, $e_{int}$ the internal energy density and $q_x$ the heat flux. The quantities $c_1,c_2,c_3$ are constants to be determined. We will consider an ideal monatomic gas in such a way that the pressure is related to the density and temperature through the ideal gas equation of state: $p=\rho k_{B}T$. In this case, the internal energy is given as $e_{int}=\frac{3}{2}\frac{k_{B}T}{m}$. The thermodynamic equilibrium points are characterized by $P_{xx}=p,q_x=0$, and we take the so-called upflow point as a reference $({\rho }_0,u_0)$, which corresponds to the cold and supersonic part of the shock. Now we define reduced variables as

$$v=\frac{u}{u_0},{\rho }^{*}=\frac{\rho }{{\rho }_0},\tau =\frac{k_{B}T}{m{u}_0^2},s=\frac{x}{\lambda },\lambda =\frac{4}{3}\frac{{\mu }_0}{{\rho }_0{u}_0},\eta =\frac{\mu }{{\mu }_0},\mathbb{P}=\frac{p_{xx}}{{\rho }_0{u}_0^2},\mathbb{Q}=\frac{q_x}{{\rho }_0{u}_0^3},$$

(6)

where the subscript 0 means the values of the quantities evaluated at the upflow. Equations (4) and (5), written in terms of the reduced variables, are given as

$${\rho }^{*}v=1,\mathbb{P}+\frac{\tau }{v}+v-1-{\tau }_0=0,v\mathbb{P}+\mathbb{Q}+\frac{1}{2}\left(5\tau -5{\tau }_0+v^2-1\right)=0;$$

(7)

their solution when $\mathbb{P}=0$ and $\mathbb{Q}=0$ are called the Rankine–Hugoniot jump conditions and determine the equilibrium points coordinates. One of them is the upflow, which is taken as the reference, and the downflow is characterized as the subsonic and hot part of the shock:

$$\begin{array}{c}\hfill \mathrm{Upflow}:v_0=1,\tau ={\tau }_0,\end{array}$$

(8)

$$\begin{array}{c}\hfill \mathrm{Downflow}:v_1=\frac{1+5{\tau }_0}{4},{\tau }_1=\frac{3+14{\tau }_0-5{\tau }_0^2}{16},\end{array}$$

(9)

where ${\tau }_0=\frac{3}{5M^2}$ and $M=\frac{u_0}{c_0}$ is the Mach number and $c_0$ the speed of sound in the fluid calculated at upflow.

Now let us calculate the entropy density change occurring between the equilibrium points taking the upflow and downflow coordinates. Since we are dealing with an ideal gas, the entropy change is well defined between the upflow and downflow coordinates and we obtain

$$\frac{m\Delta \mathcal{S}}{k_B}=\frac{m({\mathcal{S}}_{down}-{\mathcal{S}}_{up})}{k_B}=\mathrm{L}n\left[{\left(\frac{{\tau }_1}{{\tau }_0}\right)}^{3/2}\left(\frac{v_1}{v_0}\right)\right]>0.$$

(10)

As an immediate consequence, it is clear that we are dealing with an irreversible process occurring between the equilibrium points. This means that some dissipative processes take place and the shock wave must have a structure caused by them [14,15,16,17]. Now the problem turns out to be a search for the behavior of the dissipative effects present in the problem. They are represented by the viscous tensor and the heat flux due to the fact that both have already appeared in the conservation Equations (4) and (5). To go further, we recall that the system is a dilute gas, where the ideal gas equation of state and the internal energy written in terms of local variables are valid for the stationary shock wave.

Here we will have two lanes of research: First, the phenomenological scheme based on the linear irreversible thermodynamics, which starts with the fundamental entropy relation written for local variables. This means that the set of macroscopic variables describing the system satisfies the “local equilibrium hypothesis”(LEH). Then, the thermodynamic equilibrium relations are valid for local variables, the entropy density balance equation has the common entropy flux and the entropy production is the product of fluxes and thermodynamic forces [18].

The other scheme goes through the kinetic theory of gases based on the Boltzmann kinetic equation and the corresponding distribution function (df). In this case, the entropy density is directly related with the $\mathcal{H}$ Boltzmann functional as written in Equation (3) [9]. The application of the Maxwell transport equation allows for the writing of an entropy density balance equation, which in some particular cases is compatible with the LEH. Its compatibility depends on the approximation taken for the distribution function.

3. The Distribution Function

The kinetic theory approach to the study of non-equilibrium behavior for a dilute gas is based on the Boltzmann equation describing the distribution function $f(\mathbf{c},\mathbf{x},t)$ [10,19]. The kinetic scheme is not so easy to develop due to its structure and the lack of detailed knowledge about the intermolecular interaction between particles. As a consequence, some approximate methods to deal with it have been studied. Here we will restrict ourselves to the Chapman–Enskog (CE) and Grad’s moment methods; the first one is a perturbative method where the Knudsen number defined as the ratio of the mean free path in the gas and the macroscopic lenght play a very important role. On the other hand, Grad’s method is a cumulant expansion in terms of the so called peculiar velocity $\mathbf{C}=\mathbf{c}-\mathbf{u}(\mathbf{x},t)$, where $\mathbf{u}(\mathbf{x},t)$ corresponds to the average velocity in the gas [20] and $\mathbf{c}$ is the atomic velocity. Both have their own limitations. The CE method is a perturbative one, and consequently the Knudsen number must be small $(Kn\ll 1)$ to obtain a reasonable approximation. On the other hand, Grad’s method does not contain a smallness parameter; instead, a closure hypothesis is involved. Here we will consider the first order $Kn$ number expansion in the CE method which drives to the Navier–Stokes–Fourier (NSF) constitutive equations. This means that we will take the viscous tensor and the heat flux proportional to the corresponding thermodynamic forces, i.e., the velocity and temperature gradients. In contrast, Grad’s method is taken with a 13-moments closure hypothesis.

The local Maxwell distribution function is taken to start the writing of an approximate solution of the Boltzmann equation:

$$f^{(0)}=n{\left(\frac{m}{2\pi k_{B}T}\right)}^{3/2}exp\left(-\frac{m{C}^2}{2k_{B}T}\right),$$

(11)

where $n(\mathbf{x},t)$ is the local number density so that $\rho (\mathbf{x},t)=n(\mathbf{x},t)m$, and $T(\mathbf{x},t)$ is the local temperature and $C=\parallel \mathbf{C}\parallel$. The first order $(Kn)$ in the CE approximation allows us to write the df as

$$f^{CE}=f^{(0)}\left[1+\frac{m{p}_{ij}}{2p{k}_{B}T}{C}_i{C}_j-\frac{m{q}_i}{p{k}_{B}T}{C}_i\left(1-\frac{m{C}^2}{5k_{B}T}\right)\right],q_i=-\kappa {\nabla }_{i}T,p_{ij}=-2\mu \overline{({\nabla }_i{u}_j)},$$

(12)

where $\kappa$ is the thermal conductivity and $\mu$ the shear viscosity, and they depend on the intermolecular potential. We recall that the bulk viscosity vanishes for a dilute gas and $\overline{({\nabla }_i{u}_j)}$ means the traceless $ij$ components of the deformation rate tensor.

The Grad’s distribution function is written up to the thirteen moments approximation, called (G13). In this case we will have a set of 13 coupled and nonlinear equations containing the whole set of variables:

$$f^{\mathrm{G}13}=f^{(0)}\left[1+\frac{m{p}_{ij}}{2p{k}_{B}T}{C}_i{C}_j-\frac{m{q}_i}{p{k}_{B}T}{C}_i\left(1-\frac{m{C}^2}{5k_{B}T}\right)\right],$$

(13)

where the heat flux and the viscous tensor must be determined according to their transport equations, obtained directly from the Boltzmann equation. Both dfs (12) and (13) share their structure in the sense that they represent deviations of the local Maxwellian df, although they are qualitatively different due to the fact that the CE df is closed according to the first order $Kn$ number. G13 has an arbitrary closure in 13 moments, leaving the viscous tensor and the heat flux to be determined by the equations of motion. In the particular case of the plane shock wave traveling along the x-direction, they can be written as

$$f=f^{(0)}\left[1+\frac{m{p}_{xx}}{2p{k}_{B}T}\left(C_x^2-\frac{1}{2}(c_y^2+c_z^2)\right)-\frac{m{q}_x}{p{k}_{B}T}{C}_x\left(1-\frac{m{C}^2}{5k_{B}T}\right)\right],$$

(14)

which can be written in terms of the dimensionless velocity components:

$$C_{r1}=\frac{(c_x-u_0)}{u_0},c_{r2}=\frac{c_y}{u_0},c_{r3}=\frac{c_z}{u_0}.$$

(15)

Then, the dfs are written as

$$f^{(0)}=\frac{n}{u_0^3}{\left(\frac{1}{2\pi \tau }\right)}^{3/2}exp\left(-\frac{C_r^2}{2\tau }\right),$$

(16)

$$f=\frac{n}{u_0^3}{\left(\frac{1}{2\pi \tau }\right)}^{3/2}exp\left(-\frac{C_r^2}{2\tau }\right)\left[1+\frac{v\mathbb{P}}{2{\tau }^2}\left(C_{r1}^2-\frac{1}{2}(c_{r2}^2+c_{r3}^2)\right)-\frac{v\mathbb{Q}{C}_{r1}}{{\tau }^2}\left(\frac{C_r^2}{5\tau }-1\right)\right],$$

(17)

where we have taken ${\rho }^{*}v=1$ according to (7) valid for a steady shock wave. The df in Equation (17) can be written as

$$f=f^{(0)}\left(1+\varphi \right),$$

(18)

where $\varphi$ represents the deviation from the local Maxwellian df. It should be noticed that in the CE df, the deviation represented is a first order in the $Kn$ number, so $\varphi \sim \mathcal{O}(Kn)$, whereas in G13 there is not a smallness parameter.

All the kinetic calculations will be based on the df written as in Equation (18), hence we should study its behavior as a function of the velocity for some Mach number values. To proceed, we define the transversal speed $C_{rt}=\pm \sqrt{C_{r2}^2+C_{r3}^2}$ and we give some examples where we calculate the dfs for different Mach number values usually taken to calculate the shock-wave structure. First of all we should notice that the heat flux and the viscous tensor values must be given. However, they are functions of the dimensionless position s along the shock wave, which means that we must calculate them in some position across the shock. We have chosen the upflow, downflow and the center points. The center is the taken with the condition that the normalized density profile satisfies the relation $\frac{\rho (s)-{\rho }_0}{{\rho }_1-{\rho }_0}=\frac{1}{2}$. Furthermore, we need the corresponding values for $\mathbb{P},\mathbb{Q}$, which are calculated with the shock-wave NSF solution for a given Mach number. It should be noted that the density profiles depend on the Prandtl number and the viscosity model taken to solve NSF equations. In particular, we have taken the soft sphere temperature dependence in such a way that $\eta (s)={\left(\frac{\tau (s)}{{\tau }_0}\right)}^{\sigma }$. For Argon, we have fitted the value of the viscosity index ($\sigma$) so that the solutions to the NSF equations reproduce the experimental normalized density profiles provided by Alsmeyer [21,22]. The previous results have shown that the NSF equations are in good agreement with Alsmeyer’s normalized density profiles if the viscosity is enhanced; this procedure holds for Mach numbers in the range $(1.55\le M\le 9)$ provided that the viscosity index is fitted for each Mach number [22]. From the fitted NSF numerical solution we have determined the values $xx$ component of the viscous pressure tensor and the heat flux needed to evaluate the distribution function given by Equation (17). Notice that a similar procedure can in principle be carried out for Grad’s approximation, but only for Mach numbers lower than 1.65, since for higher Mach numbers, Grad’s approximation does not provide a shock-wave solution [20].

Figure 1 shows the behavior of the reduced distribution function $g=u_0^{3}f/n_0$ as a function of the velocity components $C_{r1},C_{rt}$ for Ar at $M=1.55$.

The dfs calculated at the upflow and downflow are Gaussian dfs in the corresponding equilibrium points. The df named as center describes the situation when the normalized density equals $1/2$. This point has been considered as representative of the region where the gradients in the normalized density and temperature profiles have their most important values. It has been observed that in such a region, the df may have negative values [23,24]. In Figure 1, negative values in the df are not apparent; however, in Figure 2 we provide two intersections of the distribution function with the planes $C{r}_t=0$ and $C{r}_t=2$; negative values of the distribution function are clearly exhibited in the last case.

For larger Mach numbers, the situation drastically changes, as shown in Figure 3.

On the other hand, we must notice that the calculations make sense only if the $\varphi$ expansion is convergent, as usually it is argued when $\varphi \ll 1$.

4. The Entropy Density

The entropy density is defined as

$$\rho (\mathbf{x},t)\mathcal{S}(\mathbf{x},t)=-k_B{\int }_{-\infty }^{\infty }f(\mathbf{C},\mathbf{x},t)\mathrm{L}n\left(\frac{u_0^{3}f(\mathbf{C},\mathbf{x},t)}{n_0}\right)d{\mathbf{C}}_r,$$

(19)

where $\mathcal{S}(\mathbf{x},t)$ is the specific entropy. The calculations are performed with the CE and G13 dfs as written in Equation (18). Notice that the expression in (19) contains an adimensional quantity in the logarithmic function. This was taken to consider as a reference the numerical density $n_0$, which in the case of the shock wave will be the upflow value, as for all other variables in the system. Let us define $g(\mathbf{C},\mathbf{x},t)=\frac{u_0^{3}f(\mathbf{C},\mathbf{x},t)}{n_0},g^{(0)}(\mathbf{C},\mathbf{x},t)=\frac{u_0^{3}f(\mathbf{C},\mathbf{x},t)}{n_0}$ to shorten the notation.

$$\begin{array}{ccc}\hfill & & \rho \mathcal{S}=-k_B\int f\mathrm{L}ngd\mathbf{C}=-k_B\int f^{(0)}(1+\varphi )\left[\mathrm{L}n{g}^{(0)}+\mathrm{L}n(1+\varphi )\right]d\mathbf{C}\hfill \\ & & =-k_B\int f^{(0)}\left[\mathrm{L}n{g}^{(0)}+\varphi -\frac{{\varphi }^2}{2}\right]d\mathbf{C}-k_B\int f^{(0)}\varphi \left[\mathrm{L}n{g}^{(0)}+\varphi -\frac{1}{2}{\varphi }^2\right]d\mathbf{C}+\cdots \hfill \\ & & =\underset{\rho {\mathcal{S}}^{loc}}{\underbrace{-k_B\int f^{(0)}\mathrm{L}n{g}^{(0)}}}d\mathbf{C}\underset{\mathrm{vanishes}}{\underbrace{-k_B\int f^{(0)}\left(1+\mathrm{L}n{g}^{(0)}\right)\varphi d\mathbf{C}}}\underset{\mathrm{non}\mathrm{trivial}}{\underbrace{-\frac{k_B}{2}\int f^{(0)}{\varphi }^{2}d\mathbf{C}}}+\cdots \hfill \\ & & \hfill \end{array}$$

(20)

The direct substitution of Equation (16) gives the expression for the local specific entropy, which in fact coincides with the equilibrium entropy where the equilibrium values are changed by their local counterparts,

$$\frac{\rho {\mathcal{S}}_{local}}{n{k}_B}=\frac{m{\mathcal{S}}_{local}}{k_B}=\left[\mathrm{L}n\left(v{\tau }^{3/2}\right)+\mathrm{a}\mathrm{constant}\right].$$

(21)

It should be noticed that the change of the local entropy between two points eliminates the constants, and we can then calculate the change in upflow–downflow in the shock wave. Furthermore, a different reference point can be chosen to measure the specific entropy; for example, we can take the entropy change between the upflow and any coordinate x before the downflow.

The second term in Equation (20) vanishes due to the compatibility conditions. The last term gives the first non trivial contribution consistent with the approximation ${\varphi }^2$ we have considered. The corresponding calculation is direct, although somewhat cumbersome; then

$$\frac{\rho \Delta \mathcal{S}}{n{k}_B}=\frac{m\Delta \mathcal{S}}{k_B}=\frac{\rho \Delta {\mathcal{S}}_{local}}{n{k}_B}-\frac{v^2}{{\tau }^2}\left(\frac{3{\mathbb{P}}^2}{8}+\frac{{\mathbb{Q}}^2}{5\tau }\right)+...$$

(22)

Equation (22) shows that the entropy density has two completely different terms; the first one comes from the local contribution in the df. In contrast, the second part contains quadratic contributions in the dissipative effects. If we consider the CE expansion up to first order in the $Kn$ number we notice that the local part is independent of the Knudsen number approximations. However, the terms containing the dissipative fluxes are second order $Kn$ contributions. It is not the case when we consider the description in the Grad’s moments approximation; then the second term is the simplest non-trivial contribution to the entropy density.

The entropy flux will be defined as follows:

$${\mathbf{J}}_{\mathcal{S}}=-k_B\int \mathbf{C}f\mathrm{L}n\left(\frac{u_0^{3}f}{n_0}\right)d\mathbf{C},$$

(23)

where we only need the x component, which can be written as follows:

$$\begin{array}{ccc}\hfill & & \frac{J_{\mathcal{S}x}}{k_B{u}_0^4}=-\int C_{rx}f\mathrm{L}ngd{\mathbf{C}}_r=-\int C_{rx}{f}^{(0)}(1+\varphi )\left[\mathrm{L}n{g}^{(0)}+\mathrm{L}n(1+\varphi )\right]d{\mathbf{C}}_r\hfill \\ & & =-\int C_{rx}{f}^{(0)}\left[\mathrm{L}n{g}^{(0)}+\varphi -\frac{{\varphi }^2}{2}+...\right]d{\mathbf{C}}_r=-\int C_{rx}{f}^{(0)}\varphi \left[\mathrm{L}n{g}^{(0)}+\varphi -\frac{1}{2}{\varphi }^2+...\right]d{\mathbf{C}}_r\hfill \\ & & \underbrace{=-\int C_{rx}{f}^{(0)}\mathrm{L}n{g}^{(0)}d{\mathbf{C}}_r}\underbrace{-\int C_{rx}{f}^{(0)}\mathrm{L}n{g}^{(0)}\varphi d{\mathbf{C}}_r}\underbrace{-\frac{1}{2}\int C_{rx}{f}^{(0)}{\varphi }^{2}d{\mathbf{C}}_r}+...\hfill \\ & & vanishesnon-trivialnon-trivial\hfill \end{array}$$

(24)

Now we calculate the terms in Equation (24) and define the dimensionless entropy flux in the x direction giving the following result:

$$\frac{J_{\mathcal{S}x}}{n_0{k}_B{u}_0}=j_s=\frac{\mathbb{Q}}{\tau }\left(1-\frac{2}{5}\frac{v\mathbb{P}}{\tau }+...\right),$$

(25)

valid up to the same approximation as the entropy density. It must be noted that Equation (22) is obtained from the local distribution function labeled as $f^{(0)}$, and the result is consistent with the entropy density in equilibrium evaluated with the local variables. Furthermore, due to the compatibility conditions, the entropy density can be calculated with f instead of $f^{(0)}$, and both results coincide. On the other hand, the first term in the entropy flux structure coincides with its usual expression in linear irreversible thermodynamics (LIT) [18] and it contains the heat flux over the temperature $\frac{\mathbb{Q}}{\tau }$. This contribution is the first non-trivial term in the entropy flux and, in the CE approximation, is a first order $Kn$ contribution. The second term corresponds to a correction which contains the product of dissipative fluxes, being a second order in the Knudsen number.

4.1. Navier–Stokes–Fourier (NSF)

To begin the analysis we will consider the first order $Kn$ number in the df according to the CE method. In this case the constitutive equations for the dissipative constributions are given be the usual Navier–Newton and Fourier equations [10,19,23]. The NSF set of equations are then obtained directly by the direct substitution of the constitutive equations for the viscous tensor and the heat flow, as specified in Equation (12), all of them written in terms of the dimensionless variables defined in Equation (6). Then

$$\begin{array}{c}\hfill \frac{\tau }{v}+v-\eta v^{\prime }={\tau }_0+1,\end{array}$$

(26)

$$\begin{array}{c}\hfill -v\eta v^{\prime }-\frac{15\eta }{8Pr}{\tau }^{\prime }+\frac{1}{2}(5\tau +v^2)=\frac{1}{2}(5{\tau }_0+1).\end{array}$$

(27)

where $v^{\prime },{\tau }^{\prime }$ mean the velocity and temperature derivatives with respect to s, $Pr=\frac{c_p{\eta }_0}{{\kappa }_0}=\frac{2}{3}$ is the Prandtl number and the reduced viscosity $\eta$ is modeled through a power of the temperature:

$$\eta =\frac{\mu }{{\mu }_0}={\left(\frac{\tau }{{\tau }_0}\right)}^{\sigma }.$$

(28)

The viscosity index ($\sigma$) can be determined by fitting experimental viscosity values or by using calculated viscosities obtained from formulas of the kinetic theory of gases with the interaction potential obtained from ab initio calculations.

Now we take the local equilibrium entropy density balance calculated as follows:

$$\begin{array}{c}\hfill \frac{d}{dx}\left(\rho {\mathcal{S}}_{local}u+J_{\mathcal{S}x}^{{}^{loc}}\right)=\frac{n_0{k}_B{u}_0}{\lambda }\frac{d}{ds}\left[Ln(v{\tau }^{3/2})+\frac{\mathbb{Q}}{\tau }\right]=\\ \hfill \frac{n_0{k}_B{u}_0}{\lambda }\left(-\frac{\mathbb{Q}{\tau }^{\prime }}{{\tau }^2}+\frac{{\mathbb{Q}}^{\prime }}{\tau }+\frac{{\tau }^{3/2}{v}^{\prime }+(3/2)v\sqrt{\tau }{\tau }^{\prime }}{v{\tau }^{3/2}}\right).\end{array}$$

(29)

According to the conservation laws valid in the shock wave (see Equation (7)), it is possible to write

$$\begin{array}{c}\hfill {\mathbb{Q}}^{\prime }=-\frac{3{\tau }^{\prime }}{2}-\left[1+{\tau }_0-v\right]v^{\prime },\mathbb{P}=\frac{(1+{\tau }_0)v-v^2-\tau }{v},\\ \hfill \frac{d}{dx}\left(\rho {\mathcal{S}}_{local}u+J_{\mathcal{S}x}^{{}^{loc}}\right)=\frac{n_0{k}_B{u}_0}{\lambda }\left(-\mathbb{Q}\frac{{\tau }^{\prime }}{{\tau }^2}+\frac{\tau +v^2-v(1+{\tau }_0)}{v\tau }{v}^{\prime }\right)=\end{array}$$

(30)

$$\begin{array}{c}\hfill \frac{n_0{k}_B{u}_0}{\lambda }\left(-\frac{{\tau }^{\prime }}{{\tau }^2}\mathbb{Q}-\frac{v^{\prime }}{\tau }\mathbb{P}\right).\end{array}$$

(31)

Therefore, the reduced entropy balance equation is

$$\frac{d}{ds}\left({\mathcal{S}}_{local}^{*}+j_{{}_{\mathcal{S}x}}^{loc}\right)={\Sigma }_{\mathcal{S}}^{*}{|}_{local},$$

(32)

where

$${\mathcal{S}}_{local}^{*}=\frac{m{\mathcal{S}}_{local}}{k_B},j_{\mathcal{S}x}^{loc}=\frac{\mathbb{Q}}{\tau },{\Sigma }_{\mathcal{S}}^{*}{|}_{local}=\frac{\lambda }{n_0{k}_B{u}_0}{\Sigma }_{\mathcal{S}}=-\frac{\mathbb{Q}{\tau }^{\prime }}{{\tau }^2}-\frac{\mathbb{P}{v}^{\prime }}{\tau }.$$

(33)

The entropy production density in its dimensionless expression is given as

$${\Sigma }_{\mathcal{S}}^{*}{|}_{local}=-\frac{1}{\tau }\left(\frac{\mathbb{Q}{\tau }^{\prime }}{\tau }+\mathbb{P}{v}^{\prime }\right).$$

(34)

We recall that the thermodynamic forces associated with the heat flux and viscous tensor are given as $X_{\tau }=\frac{d}{ds}\left(\frac{1}{\tau }\right),X_v=-\frac{1}{\tau }\frac{dv}{ds}$, respectively. The expression given in Equation (32) corresponds to the entropy density balance where the local equilibrium hypothesis is taken as it is in the usual linear irreversible thermodynamics approach [18]. In particular, the entropy production is the product of thermodynamic forces by fluxes.

$${\Sigma }_{\mathcal{S}}^{*}{|}_{local}=X_{\tau }\mathbb{Q}+X_v\mathbb{P}.$$

(35)

To calculate the numerical values we need the solution for the local variables $v(s),\tau (s)$ and the corresponding fluxes $\mathbb{P},\mathbb{Q}$. In this case we must take the NSF equations as written in Equations (26) and (27) to calculate the temperature and velocity profiles for a given Mach number and the viscosity index for a particular substance.

4.2. The Grad’s 13-Moment Approximation G13

Now we will work with G13 approximation, and as a first step we write the equations to be solved in this approximation. To do that we take the conservation equations, Equation (7) with the G13 equations, to determine the behavior of the viscous tensor and the heat flux, namely, [20,25],

$$\begin{array}{c}\hfill \frac{3}{4}\frac{d}{ds}(v\mathbb{P})+\frac{2}{5}\frac{d\mathbb{Q}}{ds}+\mathbb{P}\left(\frac{\tau }{\eta v}+v^{\prime }\right)+\frac{\tau }{v}{v}^{\prime }=0,\end{array}$$

(36)

$$\begin{array}{c}\hfill v{\mathbb{Q}}^{\prime }+(\tau -v\mathbb{P})\frac{d\mathbb{P}}{ds}+\frac{5\tau }{2v}{\tau }^{\prime }+\left[\frac{7}{2}{\tau }^{\prime }-v\left(\frac{d}{ds}\frac{\tau }{v}\right)\right]\mathbb{P}+\left(\frac{16}{5}{v}^{\prime }+\frac{8\tau }{9\eta v}\right)\mathbb{Q}=0.\end{array}$$

(37)

The set of Equations (7), (36) and (37) can be solved numerically to obtain the corresponding profiles with the same model for viscosity, as in the NSF case.

As a second step, the G13 df is taken to calculate the entropy density, which will contain the already calculated local part and some additional terms coming from the non-local contributions; both are shown in Equations (20)–(22). It should be noted that the non-local contributions are nonlinear terms in the viscous tensor and the heat flux. Furthermore, the entropy flux will have the local part and the non-local terms, as shown in Equation (25), where the non-local terms are bilinear in the viscous tensor and the heat flux. Now we will take the higher order terms, $O({\varphi }^2)$, and calculate their derivative with respect to the dimensionless position; hence, using the results from the Appendix A, we obtain that the entropy balance equation is then written as follows:

$$\begin{array}{c}\hfill \frac{d}{ds}\left(\rho S_{total}u+J_{{}_{\mathcal{S}x}}^{total}\right)=\frac{\lambda }{n_0{k}_B{u}_0}\{\frac{1}{\tau }\left(-\frac{\mathbb{Q}{\tau }^{\prime }}{{\tau }^2}-\frac{\mathbb{P}{v}^{\prime }}{\tau }\right)\\ \hfill +\frac{1}{4\tau }\left(-\frac{7v^{\prime }}{5}+\frac{3v^2{v}^{\prime }}{\tau }+\frac{27v{\tau }^{\prime }}{5\tau }\right)\mathbb{P}+\frac{1}{5\tau }\left(\frac{2{\tau }^{\prime }}{\tau }+\frac{4v{v}^{\prime }}{\tau }-\frac{2v^{\prime }}{v}+\frac{3v^2{\tau }^{\prime }}{{\tau }^2}\right)\mathbb{Q}\\ \hfill +\frac{v}{4{\tau }^2}\left(\frac{3v{\tau }^{\prime }}{\tau }-\frac{7v^{\prime }}{5}\right){\mathbb{P}}^2+\frac{1}{5{\tau }^3}\left(\frac{3v^2{\tau }^{\prime }}{\tau }-2v{v}^{\prime }\right){\mathbb{Q}}^2+\frac{2}{5{\tau }^2}\left(-v^{\prime }+\frac{v^2{v}^{\prime }}{\tau }+\frac{2v{\tau }^{\prime }}{\tau }\right)\mathbb{P}\mathbb{Q}\},\end{array}$$

(38)

where $\rho S_{total}u$ is given in Equation (22) and the entropy flux $J_s^{total}$ is written in Equation (25).

A close view of the entropy balance Equation (38) shows us that it contains the usual local part and several terms which have a different structure, in the sense that they are not the products of fluxes and thermodynamic forces as it happened in the NSF case (see Equation (35)). In fact, we know that the thermodynamics forces and their corresponding fluxes definitions are not unique due to the existence of Meixner transformations [26]. Hence, we wonder if it is possible to redefine the thermodynamic forces and recover the mentioned structure; however, it is easily seen that there are some kind of crossed effects which prevent such efforts. Furthermore, there are additional nonlinear contributions which make the problem worse.

5. The Local Equilibrium Hypothesis (LEH) and Beyond

In order to go a step further in the analysis of the calculations performed in the last section, we will specify as clearly as possible the content and implications of the local equilibrium hypothesis. First of all, it is important that it has been the cornerstone in the development of linear irreversible thermodynamics [18] and a lot of irreversible studies and generalizations [27,28,29,30,31,32,33,34]. It takes the classical thermodynamics relations to establish that the local variables describing a situation out of thermodynamic equilibrium can be written in the same way as they stand in thermodynamic equilibrium. This means that the equation of state, the caloric equation as well as the fundamental $Td\mathcal{S}$ relation are valid in local equilibrium. The line of thinking behind it focus on some situations which are not too far from true equilibrium, no matter that we do not have a quantitative criterion to be sure about such a fact. The NSF set of equations has kinetic support with the CE expansion in the Knudsen number $Kn$ as we noticed before; the entropy density calculation must then be taken only with what we called the local term due to the fact that quadratic terms in the dissipative fluxes $\mathbb{P},\mathbb{Q}$ become of the order $K{n}^2$. The entropy flux at order $Kn$ is just proportional to $\mathbb{Q}/\tau$, and the entropy production up to the same order corresponds to the product of fluxes and the thermodynamic forces, as identified in the literature [18,35]. It should be mentioned that according to the LEH, the entropy production density must be positive definite for any position along the shock wave. This condition is completely fulfilled with the NSF set of equations. Sometimes it is taken to extract conclusions about the transport coefficients, shear viscosity, thermal conductivity, diffusion coefficients, etc., when they are needed. In this context we see that the NSF set of equations is completely consistent with the LEH up to the first order in the Knudsen number.

Concerning the G13 equations, we must recall that their kinetic support is given through an arbitrary closure hypothesis and they do not have a smallness parameter to give them a systematic way to characterize the kind of approximation they represent. With the calculations carried out before, we have a starting point to analyze the results. First of all, the calculations were performed with the ${\varphi }^2$ contributions. The NSF calculations have the same structure; however, the Knudsen number order has given us a way to classify each term and neglect $K{n}^2$ contributions. The lack of the smallness parameter in G13 avoids such a classification and we must consider the complete expression. Some terms are consistent with LEH, but some others are not. In this case we can say that we have non-local contributions which are not negligible. Obviously, the entropy production density has a different structure than the one consistent with LEH.

It should be pointed out that the G13 equations may have a phenomenological interpretation through the extended irreversible thermodynamic (EIT) description. EIT considers the fluxes as relevant and independent variables and as constituting the basis to write a generalization of the $Td\mathcal{S}$ fundamental relation [34]. The steps followed by EIT drive an entropy production with the same structure as we have found for the G13 set of equations [36]. This means that the LEH is not obeyed in the EIT formulation. A similar quotation can be made for Burnett equations, which rise as an approximate solution for the Boltzmann equation at the second order in the $Kn$ number [10]. The Burnett’s equation’s origin is completely kinetic and derived from the CE method; however they have a thermodynamic context, as shown in Reference [37]. In both cases, G13 and Burnett, among other sets of equations, are not compatible with LEH, although they correspond to a thermodynamic context.

The numerical results will be given us a clear idea about the relevance of the non-local contributions.

6. Numerical Results

Let us study the performance of NSF and G13 in relation to the entropy density behavior. In order to compare the models the normalized variables are defined as follows:

$${\rho }_n(s)=\frac{1/v(s)-1/v_0}{1/v_1-1/v_0},{\tau }_n(s)=\frac{\tau (s)-{\tau }_0}{{\tau }_1-{\tau }_0}.$$

(39)

Figure 4 shows the normalized and temperature profiles for $He$ at $M=1.59$, and to reduce the variables we have taken the following:

$$\begin{array}{ccc}\hfill u_0& =& 1183.3884\frac{\mathrm{mt}}{\sec },{\rho }_0=1.93\times {10}^{-8}\frac{\mathrm{g}r}{{\mathrm{cm}}^3},k_B=1.380649\times {10}^{-23}\frac{\mathrm{J}}{\mathrm{K}},\hfill \\ \hfill m& =& 6.6465\times {10}^{-27}\mathrm{kg},{\mu }_0=13.0310\times {10}^{-6}\mathrm{Pa}·\mathrm{s}\Rightarrow \hfill \\ \hfill \lambda & \approx & 0.7607\mathrm{mm},{\tau }_0=\frac{3}{5M^2}\Rightarrow {\tau }_0=2000/8427\mathrm{for}M=1.59.\hfill \end{array}$$

(40)

In addition, we used the soft sphere model, given by Equation (28), to fit the values of the transport coefficients of ${}^4$He reported by Hurly [38] using ab initio values for the interatomic potential. The resulting value for the viscosity index, obtained from the fit, is $\sigma =0.6716$ [25]. With the previous information, the solutions to the NSF differential equations, given by Equations (26) and (27), or to Grad’s equations, given by Equations (36) and (37), were determined by electing the normalized density profile with the value 1/2 for $s\approx 18.3$ [25,39].

The normalized density and temperature profiles and their comparison with experimental data for G13 and NSF have been discussed in the literature [25,39], as shown in Figure 4. The results represented here are shown to demonstrate that they have the usual behavior in both NSF and G13 models. The viscous tensor and the heat flux are shown in Figure 5 for both models. It is clear that they vanish in the equilibrium points, as it should be.

Figure 6 shows the entropy density change profile in the cases corresponding to the NSF local expression; the local contribution in G13 model and the complete entropy density change for the G13 case are calculated with the solution for the G13 equations. The comparison between the local and the total entropy change calculated with the G13 model solutions shows how relevant the non-local terms become. It shows in a clear way that the non-local contributions in the G13 approximation are not negligible. A conspicuous feature of the local entropy density profile is that it exhibits a maximum, as shown in Figure 6. The entropy discussed in the work by Margolin et al. [17] is monotonic, which is expected outside of equilibrium according to some researchers. However, the results of numerically solving the Boltzmann equation by Malkov et al. [40] report non-monotonic entropy profiles. For other theories in which the entropy in a shock wave is discussed, the interested reader is referred to the relevant bibliography [41,42].

In Figure 7, the usual entropy flux for the NSF model is shown in contrast with the local contribution in the G13 case where a big difference is evident. In this case, the non-local contribution in the G13 approximation becomes the leading part.

The entropy production is presented and compared in Figure 8, where it is shown that the non-local contribution plays a very significative role in the total entropy production. In this case, it is seen that its difference with the total entropy production is somewhat small.

7. Concluding Remarks

There is no doubt about the entropy concept difficulties, much of them caused because an incomplete context is given when it is applied to a particular problem. Even more, the local equilibrium hypothesis must be applied carefully, since it has somewhat blurry limitations. In particular, sometimes it is extended indiscriminately, taking as a basis the second law of thermodynamics for an isolated system, without questioning its validity under the particular conditions imposed by the problem. The work presented here has tried to show explicitly the contrast between the entropy calculations based on the NSF and the Grad13 models, within the frame of their application to study the shock-wave structure in dilute gases. The NSF model is consistent with the LEH, and the Grad13 moments do not perform in the same way. In fact, the calculation shown allows a quantitative estimate of how different the entropy production is in this case. The differences come from the numerical solutions for the variables $(v(s),\tau (s))$ in each model; however, the difference in the structure is more important, as it is shown in Equation (38). In particular, the last three terms are second order in the fluxes, and due to the lack of a smallness parameter, their order of magnitude is not known. The entropy production in the G13 model does not have the structure required by the LHE. Moreover, its non-local contributions cannot be neglected, at least in this particular example. It must be said that Grad’s moment method is just one model to show the problem set in this work; for example, the Burnett set constitutes another set which is not consistent with LEH.

It should be noted that the shock-wave problem with the NSF or G13 equations can be solved independently of the definition of entropy used. For a glimpse of the vast literature available on shock waves in dilute gases, the reader can refer to the bibliography cited in reference [43], which also provides a few references to shock waves in dense gases. For the latter case, the interested reader may take a look at chapter 6 of the book by Hoover and Hoover [1].