Shockwave Compression and Joule-Thomson Expansion

Structurally-stable atomistic one-dimensional shockwaves have long been simulated by injecting fresh cool particles and extracting old hot particles at opposite ends of a simulation box. The resulting shock profiles demonstrate tensor temperature, with the longitudinal temperature exceeding the transverse, and Maxwell's delayed response, with stress lagging strainrate and heat flux lagging temperature gradient. Here this same geometry, supplemented by a short-ranged external"plug"field, is used to simulate steady Joule-Kelvin throttling flow of hot dense fluid through a porous plug, producing a dilute and cooler product fluid.

Stationary One-Dimensional Shockwaves-Shockwaves are arguably farther from equilibrium than are any other readily available states of a nonequilibrium fluid [1][2][3][4][5][6][7][8][9] . In just a few collision times, or mean free paths, the shock transforms cold equilibrium fluid (or solid) into a hot compressed state [1][2][3][4][5] . Laboratory shockwaves at a few terapascals can compress condensed matter as much as threefold, to densities and pressures far greater than those at the center of the earth 6 . Because the shock transformation is a steady small-scale continuous process, converting kinetic energy to internal energy without any external heating, steadystate shockwave structures can be replicated with computer simulations [2][3][4][5][7][8][9] . The inset of Figure 1 shows an interior snapshot of a typical simulation, with cold particles entering at the left and hot ones exiting to the right. The corresponding density profile snapshot using Lucy's weight function for the spatial averaging 7-11 is the smooth curve. The shockwidth can be estimated from the maximum slope. It is just a few atomic spacings. Steadystate profiles generated in this way are fully consistent with the transient profiles generated with [1] shrinking periodic boundaries or [2] headon collisions of two similar blocks of cold material 5,7-9 .
Both experiments and simulations show that initially-sinusoidal shockfronts soon become planar. Steady shockwaves are accurately one-dimensional 3,5,7 . Accordingly, the mass, momentum, and energy fluxes ( in the x direction, the propagation direction ) are all constant 2 in the comoving coordinate frame of Figure 1, the frame moving with the shockwave 1 : ρu ; P xx + ρu 2 ; (ρu)[ e + (P xx /ρ) + (u 2 /2) ] + Q x ; All Three Fluxes Constant .
Here ρ(x) and u(x) are the mass density and the flow velocity, P xx (x) is the pressure-tensor component in the propagation direction. e(x) is the internal energy per unit mass and Q x is the heat flux vector, measuring the conductive flow of heat in the comoving frame.
The cold entrance velocity is +u s (the "shock velocity") and the hot exit velocity is +(u s − u p ) (where u p is the "particle velocity") in the hot fluid. Away from the shockfront the cold and hot pressure and energy have their thermodynamic equilibrium values : Eliminating u s and u p from the three constant-flux equations gives the "Hugoniot Equation" or "shock adiabat", ∆e = P ∆v , where P is the mean pressure, [ P cold + P hot ]/2 , and ∆v is the overall change in volume per unit mass, (1/ρ) cold − (1/ρ) hot . Though there is no external heating there is heat flow within the shockwave structure. For weak shockwaves it is given The limiting values of the energy flux divided by the mass flux far from the shockwave are equal : In shockwaves the inflow is supersonic so that the kinetic energy cannot be ignored. Choosing the initial thermodynamic state along with the particle velocity determines the shock velocity as well as the pressure and energy of the resulting "hot" state.
Although such a potential was perfectly satisfactory for the shockwave simulations of twofold compression it suggests the possibility of poor behavior at high density, where the force is a decreasing function of compression. Accordingly we compared results with a modified pair potential for which the force remains constant, with its maximum value F max at separations less than r max = (1/7) = 0.377964473 : φ max (r < r max ) = (6/7) 4 + F max (r − r max ) ; F max (r < r max ) = 8(6/7) 3 (1/7) . Joule-Thomson profiles including this φ pair precaution weren't significantly changed from those with the unmodified potential.
Corresponding time-averaged density and velocity profiles are shown in Figure 3, along with the (necessarily constant) mass flux, ρu . Just as in our shock work the one-dimensional grid-profile averages were all computed using Lucy's one-dimensional smooth-particle weight function 7-11 , with h = 3 : With an input speed of 0.5, which quickly accelerates to 0.62, the velocity speeds up to 1.25 on passing through the plug potential. Straightforward Runge-Kutta simulation converges relatively simply and quickly to a flow satisfying slight modifications of the conservation relations which hold for shockwave simulations.
The longitudinal momentum flux drops at the barrier because the barrier force removes momentum. The overall flux drop exactly matches (F barrier /L y ) , with Mass, momentum, and energy fluxes are shown in Figure 3.
The energy flux is particularly interesting. Adding the contributions of pair interactions (ẋ i +ẋ j 2 )x ij F x ij to the "convective flux" gives perfect agreement between the entrance and exit flows. These contributions can be divided equally between Particles i and j. Alternatively, they can be velocity-weighted: (ẋ i 2 )x ij F x ij for i and (ẋ j 2 )x ij F x ij for j. The effect of this choice on the energy flux is insignificant, of order 0.001. In shockwaves the total pressure-tensor component P xx includes the ρkT xx which is absent in our Joule-Thomson flux. The derivations for these two slightly different expressions for the energy flux are both familiar textbook fare 15 . The reason for the difference is interesting. The x component of the purely kinetic part of the energy flux (excluding the contributions from φ and F ) involves local sums cubic in the velocity components. In the equilibrium case the cubic sum can be expressed in terms of the stream velocity and the deviations from it, which can in turn be expressed in terms of temperature : (v x /2)(v 2 x + v 2 y ) = (1/2) (u + δv x ) 3 + u(δv y ) 2 = (1/2)u 3 + (3/2)ukT xx + (1/2)ukT yy .

6
The resulting "extra" ρukT xx can, if desired, be combined with the potential part of uP xx so as to agree with the continuum energy-flux expression. Far from equilibrium this simplification does not hold and the full cubic kinetic-theory sums must be evaluated.
In general, it is interesting to note that the hot and cold momentum fluxes don't match in the Joule-Kelvin experiment though they do in the shockwave. ( A series of Joule-Thomson states can be generated by using several plug barriers rather than just one. Accordingly, we believe that they will, like shockwaves, provide a useful source of computer-experimental constitutive information for flow states far from equilibrium and help in choosing the optimum weight function for correlating microscopic and macroscopic flow descriptions.