Theory and simulation of shock waves: Entropy production and energy conversion

We have considered a shock wave as a surface of discontinuity and computed the entropy production using nonequilibrium thermodynamics for surfaces. The results from this method, which we call the “Gibbs excess method” (GEM), were compared with results from three alternative methods, all based on the entropy balance in the shock-front region, but with different assumptions about local equilibrium. Nonequilibrium molecular dynamics (NEMD) simulations were used to simulate a thermal blast in a one-component gas consisting of particles interacting with the Lennard-Jones/spline potential. This provided data for the theoretical analysis. Two cases were studied, a weak shock with Mach number $M\approx 2$ and a strong shock with $M\approx 6$ and with a Prandtl number of the gas $\text{Pr}\approx 1.4$ in both cases. The four theoretical methods gave consistent results for the time-dependent surface excess entropy production for both Mach numbers. The internal energy was found to deviate only slightly from equilibrium values in the shock front. The pressure profile was found to be consistent with the Navier-Stokes equations. The entropy production in the weak and strong shocks were approximately proportional to the square of the Mach number and decayed with time at approximately the same relative rate. In both cases, some 97% of the total entropy production in the gas occurred in the shock wave. The GEM showed that most of the shock's kinetic energy was converted reversibly into enthalpy and entropy, and a small amount was dissipated as produced entropy. The shock waves traveled at almost constant speed, and we found that the overpressure determined from NEMD simulations agreed well with the Rankine-Hugoniot conditions for steady-state shocks.

Figure 1

Illustration of the entropy density, ${\rho }_{\text{s}}$, (solid line) and the entropy production, ${\sigma }_{\text{s}}$, (dashed line) around a shockwave moving from left to right with velocity $v^{\text{s}}$. The dash-dot line is a macroscopic illustration of the shockwave with discontinuous fluid properties.

Figure 2

Gibbs' equal-area construction for determination of shock-front position, $\ell$, by requiring that the surface excess mass density ${\rho }^{\text{s}}=0$ (the two shaded areas in the insert are equal). The black circles are results from the NEMD simulations for the strong shock as described in Sec. 5. The vertical dashed line shows the position of the equimolar surface, and the dashed-dot lines are least-squares fit to the bulk data. The uncertainties in the NEMD data, determined as three standard errors, are shown as vertical bars.

Figure 3

MD cell layout. The bright orange regions at the ends show where the energy was added during the heat pulse. The aspect ratio was different from that shown in the illustration.

Figure 4

Shock-front position as determined with the Gibbs' equal-area construction as function of time. The dotted lines are third-order polynomial fits.

Figure 5

Wave speed shown as Mach number as function of time for the two cases shown in Fig. 4.

Figure 6

Profile of $\mathsf{\text{T}}=\frac{1}{3}({\mathsf{\text{T}}}_{xx}+{\mathsf{\text{T}}}_{yy}+{\mathsf{\text{T}}}_{zz})$ for the strong shock at time $t^{*}=600$. The insert shows that the normal and tangential components of the temperature tensor are different in the shock front, but equal immediately behind the shock. The uncertainties are three standard errors. The vertical dashed line shows the position of the equimolar surface.

Figure 7

Particle speed distributions for the weak and strong shocks. The data for the weak shock were recorded at $x^{*}=3434$ and $t^{*}=1000$, and for the strong shock at $x^{*}=6204$ and $t^{*}=600$. The weak shock shows a perfect Maxwell-Boltzmann distribution, whereas the strong shock does not. The uncertainties are three standard errors.

Figure 8

Specific internal energy as function of $x$ in the front region of the strong shock at $t^{*}=600$. The dots show the total internal energy determined by NEMD, $u$(total), and the squares show the configurational (potential) part of it, $u_p$. The black and white squares show data from the nonequilibrium and equilibrium simulations, respectively. Note that the configurational contributions (referring to the right axis) are so small that they do not visibly separate the kinetic contributions to the internal energy from the total in $u$. The error bars are three standard errors, the errors for the equilibrium results are smaller than the symbol size. The vertical dashed line shows the position of the Gibbs equimolar surface.

Figure 9

Contributions to the pressure for $5900\le x^{*}\le 6400$ at $t^{*}=600$ for the strong shock. The vertical dashed line shows the position of the Gibbs equimolar surface at $x^{*}=6154$. (a) The black dots and triangles represent the total nonequilibrium pressure, $\frac{1}{3}\text{Tr}\left(\mathsf{\text{P}}\right)$, and its kinetic contributions, $p_{\text{ig}}$, respectively. The squares show the configurational nonequilibrium and equilibrium contributions. Note that these contributions refer to the axis on the right. (b) The $x$-component of the pressure tensor (${\mathsf{\text{P}}}_{xx}$) and the viscous pressure computed from the NEMD data (${\mathrm{\Pi }}_{xx}^{\text{MD}}$) and Navier-Stokes (${\mathrm{\Pi }}_{xx}^{\text{NS}})$.

Figure 10

Data used in the BBM for the strong shock at $t^{*}=600$. Uncertainties have not been determined, but are probably similar to the scatter on the left of the wave. (a) The functions $\frac{\partial }{\partial t}{\rho }_{\text{s}}(x,t)$ and $\frac{\partial }{\partial x}{J}_{\text{s}}(x,t)$ and their sum ${\sigma }_{\text{s}}(x,t)$. The vertical dashed line is the position of the Gibbs equimolar surface. (b) The integrand ${\sigma }_{\text{s}}(x,t)$ and the fitted Gaussian. We also show results from the LIT method used by Velasco and Uribe [30], i.e., the Gibbs equation used locally for each control volume in the simulation (see text).

Figure 11

(a) Surface excess entropy density and internal energy density as function of time for the strong shock. (b) Surface excess internal energy density vs surface excess entropy density. The slope was determined to $T^{\text{s}*}=3.9$.

Figure 12

(a) Plots of ${\sigma }_q$ and ${\sigma }_j$ for the strong shock as determined from Eqs. (32) and (33), respectively, at $t^{*}=600$. The vertical dashed line shows the position of the Gibbs equimolar surface. (b) The different terms in the parentheses of Eq. (33) as function of $x$ at $t^{*}=600$. The line is the sum of the four terms.

Figure 13

Surface excess entropy production as function of time for the strong (a) and weak (b) shocks, with Mach numbers $\approx 6$ and $\approx 2$, respectively. The lines show data from the four methods used and the line with the circles are exponential functions fitted to all the four data sets.

Figure 14

Peak overpressure as function of time for the strong and weak shocks. The dashed lines represent the RH conditions [Eq. (58)], and the solid lines are results from the NEMD simulations.