Nonequilibrium thermodynamics of surfaces captures the energy conversions in a shockwave

The local entropy production in a shock wave was analysed in the framework of non-equilibrium thermodynamics (NET) of surfaces. We show that the thermodynamic state variables in the shock front are equal to their equilibrium values, despite lack of global equilibrium in the dense gas. This observation was used to develop a theory for the entropy production in a shock wave using Gibbs' surface excess properties. The theoretical results were compared with a numerical evaluation of the entropy balance for the shock front and confirmed by non-equilibrium molecular dynamics (NEMD) simulations. The NET analysis shows that the dominant contribution to the entropy production is the dissipation of kinetic and compression energy. This opens the door to accurate representations of energy conversions in shock waves.


Abstract:
The local entropy production in a shock wave was analysed in the framework of non-equilibrium thermodynamics (NET) of surfaces. We show that the thermodynamic state variables in the shock front are equal to their equilibrium values, despite lack of global equilibrium in the dense gas. This observation was used to develop a theory for the entropy production in a shock wave using Gibbs' surface excess properties. The theoretical results were compared with a numerical evaluation of the entropy balance for the shock front and confirmed by non-equilibrium molecular dynamics (NEMD) simulations. The NET analysis shows that the dominant contribution to the entropy production is the dissipation of kinetic and compression energy. This opens the door to accurate representations of energy conversions in shock waves.
The basic theory for shock waves was developed in the late 19th century by Rankine and Hugoniot. 1,2 Important theoretical developments were made during and after the second world war, 3 and by use of computer simulations in more recent years, 4 so shock waves are now pretty well understood. There are, however, remaining questions, such as exactly how the kinetic and compression energy carried by a shock wave is dissipated or converted to other forms. This is an important question in fields like detonations, 5 material science, 6 and formation and collapse of bubbles, 7 to mention a few. Energy conversion is a topic of thermodynamics, and since shock waves are irreversible processes, more specifically non-equilibrium thermodynamics (NET). A shock front has a sharp gradient in the density, similar to a liquid-vapor surface. This has led us to consider the shock front as a surface and use NET for surfaces 8 as a tool to extract detailed information about the shock wave. We show that both energy dissipation and reversible conversion can be determined from this analysis and that it gives new information about energy conversion at the shock front. We start by deriving the governing equations for the Gibbs excess method applied to a shock wave.
The Gibbs excess method. We consider a planar shock wave moving in positive x-direction (from left to right). The shock front is treated as a discontinuity represented by excess variables for the surface (in excess of the bulk phase). This is similar to the typical treatment of e.g. vapor-liquid interfaces. For example, the surface excess mass density is defined by where superscript "s" denotes a surface excess property, is the position of the surface, Θ is the Heaviside step function, and x d < and x u > are positions in the bulk phases. The superscripts "d" and "u" denote the extrapolated values of ρ(x) from the bulk values on the downstream (left) and upstream (right) side of the shock. Eq. (1) is the Gibbs definition of excess densities. 9 Furthermore, we assume that thermodynamic relations between surface variables remain valid also when the system at large is out of equilibrium, as introduced by Bedeaux, Albano and Mazur. 10,11 Many theoretical and simulation studies have showed that this assumption applies to interfaces perturbed far beyond global equilibrium. 8,12,13 In the Gibbs excess method, one must define a dividing surface. We will do this by requiring that ρ s equals zero, which determines (t). The surface moves with a velocity given by The entropy density is represented as 8 where ρ s s is the surface excess entropy density. The balance equation for entropy for the planar shock wave is where ρs, Js, and σ are the entropy density, entropy flux, and entropy production per unit volume, respectively. The entropy flux in Eq. (4) is: where J q is the measurable heat flux, T is the temperature, and vcm,x is the local streaming velocity. Substituting Eq. (3) into Eq. (4), we obtain after some algebra the following balance equation for the surface excess entropy density, i.e. the entropy that is assigned to the shock front as represented by excess and extrapolated variables: where we have used the notation for the difference across the surface. Here, Js−v s ρs = J q /T + ρs(v−v s ) is the entropy flux in the surface frame of reference. The excess entropy density, ρ s s , is found by replacing the mass density with the entropy density in Eq. (1) and using the same value of as determined by the Gibbs construction for ρ s .
We introduce next the Gibbs equation applied to the excess densities of the surface, where ρ s u is the excess internal energy density, the surface temperature is defined as T s = ∂ρ s u /∂ρ s s at constant ρ s , and µ s is the specific Gibbs energy of the surface. Note that ρ s in Eq. (8) equals zero by construction. The ρ s u was found in the same way as ρ s s , i.e. by replacing ρ by ρu in Eq. (1). The local equilibrium hypothesis in the excess description amounts to assuming that Eq. (8) is valid. 13 Rearranging Eq. (8) gives Conservation of energy across the shock leads to the following balance equation for the excess internal energy density: (10) where Pxx = p + Πxx is the xx-component of the pressure tensor (including the viscous pressure component Πxx = −( 4 3 ηS+ηB) ∂v ∂x ). All properties in the brackets are bulk properties. By introducing Eq. (10) into Eq. (9) and comparing the result with the entropy balance, Eq. (6), we obtain the following expression for the excess entropy production, using the same bracket notation as in Eq. (7): where and In Eq. The speed of sound in the gas upstream of the shock was determined from independent equilibrium molecular dynamics simulations and found to be 1.298, which is very close to the ideal-gas value of 1.291. The blast caused the shock wave to travel at a slowly retarding supersonic speed with a Mach number of 2.1. This is a weak shock, on the borderline of the validity range of the local equilibrium condition. 14,15 The question of local equilibrium. Shock waves are nonequilibrium structures. For instance, the velocity distribution and the kinetic temperature in the shock front is anisotropic. 16,17 However, many studies have confirmed that the classical local equilibrium hypothesis 8 holds when the interfacial properties are described by Gibbs excess variables. 8,12,13 In agreement with previous results, 17,18 we also found that the local kinetic temperature was anisotropic in the shock front. On the other hand, we found that the Boltzmann H-function based on the particle speeds from the NEMD data was consistent with a state of local equilibrium. This is illustrated in Fig. 1a, based on the speed of 35,996 particles (total from 20 runs) that were in the control volume of thickness ∆x * , centred at the shock wave front at x * = 3420, at t * = 1000. The fitted Maxwell-Boltzmann distribution gave a temperature T * = 1.9 ± 0.1, in good agreement with the kinetic temperature T * = 1.87 ± 0.03 (uncertainties given as three standard errors of the mean). The separate particle velocities in x-, y-and z-directions confirmed the equilibrium longitudinal and transverse distribu- tions and the corresponding local kinetic temperatures. This is illustrated in Fig. 1b. We conclude from this that the nonequilibrium entropy determined from the H-function agrees with the equilibrium value within the estimated uncertainty.
Direct numerical evaluation of the entropy balance equation. In order to verify the validity of Eqs. (11) -(13) and the Gibbs excess method, we computed the local entropy production by direct numerical evaluation of the entropy balance equation, Eq. (4), over the surface region. The only assumption behind this method is that the local properties are determined by the equation of state as discussed above. The first term in Eq. (4), (∂ρs/∂t) was evaluated by numerical differentiation of the data from NEMD using a five-point method.
The heat flux J q in Eq. (5) was computed directly from the NEMD simulations. Although the heat flux has a sharp peak in the shock front, it contributes at most only 3 % to Js in the present case. The second term in Eq. (4), (∂Js/∂x), was computed by numerical differentiation of the entropy flux profiles. The two terms on the left-hand side of Eq. (4) are large and of opposite sign around the shock front. The uncertainty in the sum at the left-hand-side of Eq. (4) is therefore large. Nevertheless, the entropy production shown in Fig. 2 displays a distinct positive peak around the shock front, in agreement with the second law of thermodynamics.
Finally, σ(x) was integrated from x * = 3000 to x * = 3800 with a simple trapezoidal rule. Noise on both sides of the peak gives positive (shaded green in Fig. 2) and negative (shaded red) contributions to the integral, which cancel out to zero. The non-zero contribution to the spacial integral of σ(x, t) is from the peak centred on the shock front. The integral determines the surface excess entropy production per cross sectional area, σ s . An estimate at t * = 1000 gave an excess entropy production in the shock wave σ s * = 0.007 ± 0.002.
Numerical evaluation of Eqs. (11) -(13). The dominant contribution to the excess entropy production comes from σj. The term σq is small because the heat flux in the bulk phases is small (slightly negative downstream of the front and zero upstream). Fig. 3a shows the profiles of the four terms in the bracket in Eq. (13) at t * = 1000. The viscous pressure term varies little over the shock front. The difference between the extrapolated values is practically zero.
The kinetic energy term includes the center-of-mass velocity relative to the shock wave velocity. This relative velocity is larger upstream than downstream, so the difference defined by the bracket is positive. The specific Gibbs energy is the difference between the specific enthalpy and the product of temperature and specific entropy, µ = h − T s. Both h and s increase when the shock wave passes, but the entropy term increases more than the enthalpy, so the total effect is a decrease in the specific Gibbs energy. The term (T − T s )ρs/ρ increases because both (T − T s ) and ρs/ρ increase from upstream to downstream. The mass flux is constant across the shock front because mass is conserved, and therefore equal to the upstream value, j = −ρv s . In total, the term σj is positive. Hence, for the propagating shock examined in this work, the overall picture is that kinetic energy and chemical energy are partially converted to entropy across the shock front.
The two contributions to the entropy production and their sum are shown in Fig. 3b. The relatively small contribution from heat conduction and viscous dissipation is a consequence of the low density of the gas, and we expect these terms to be larger in fluids with higher densities. Eqs. (11) -(13) provide a tool to quantify this.
The total excess entropy production as given by Eq. (11) is a difference between properties extrapolated to the surface position. This extrapolation is illustrated by the horizontal lines in Fig. 3b and the difference is illustrated by the double arrow. We emphasize that the values in the shock-front region have no significance in this context, only the extrapolated values are relevant. We found that the total excess entropy production per unit surface area was σ s * = 0.009 ± 0.001, which compares well with the value for the excess entropy production based on Eq. (4).
Conclusions. We have presented a new method to analyze the entropy production in a propagating shock wave by use of non-equilibrium thermodynamics for surfaces, using surface excess variables. The only assumptions behind this method is that the local properties are determined by the equation of state and that the Gibbs equation holds for the surface excess properties. These assumptions have been found valid for surfaces, and so also in the present case. A numerical evaluation was made with data from non-equilibrium molecular dynamics simulations of a weak shock wave. Within the  accuracy of the simulations, this method gave the same surface excess entropy production as a direct numerical evaluation of the local entropy balance in the shock-front region. The new method is a powerful tool for analysis of energy conversions in shock waves because it quantifies the different contributions to the excess entropy production. A consistent representation of dissipation in shocks is of key importance for the dynamic description of shock waves in a variety of fields.