Fluctuation-Dissipation Relations for the Shock Compression of Hydrodynamic Solids

I. INTRODUCTION

It has not been well recognized that von Neumann considered a fluctuation-dissipation relation for hydrodynamic shocks¹ in a technical report prior to the publication of the famous paper (with R. D. Richtmyer) that introduced the artificial viscosity to regularize shocks in hydrodynamic flows.² In the report he defined the dissipation ("degradation of energy" in his term) by the variant fluctuations of both kinetic and internal energies and calculated their magnitudes based on harmonic fluctuation oscillations.

Recently D. Grady has published a series of papers that are concerned with the fluctuation-dissipation relation for steady shock transitions in solid materials ${}^{3,4}$. In contrast to von Neumann, Grady considered a variant velocity correlation function to define energy dissipation with focus on kinetic energy. He then hypothesized an exponential correlation function and derived a dissipated energy-time relationship to characterize the shock transition. The Grady derivation of the energy-time relationship provides a foundational support for his shock invariant concept in the manner of classical action principle. The often referenced fourth power relationship between shock pressure and strain rate ${}^5$ is a direct derivative of this invariant relationship. For reference a summary of their approach is described in Section 2.

Study of fluctuations in the shock compression of solid materials is obviously not limited to the works of von Neumann and Grady. As briefly reviewed in a paper, noteworthy examples of modeling fluctuations are those of M. Baer, his associates and Y. Meshcherikov to name a few. Baer's work was concerned with the multidimensional simulation of shock wave propagation in heterogeneous media by use of CTH code. In contrast Meshcheryakov's work was concerned with the multiscale mechanics of shock compression processes with focus on the meso-macro exchange of momentum and energy. Interestingly, however, there is no discussion of fluctuation-energy relations in their works.

Other types of fluctuation modeling include the discrete element modeling of shock propagation in two dimensional anisotropic crystals ${}^{10}$ and a turbulence modeling of entropy change ${}^{11}$. Among the works on fluctuations in shock waves, unique is the work of L. Margolin. He investigated the meaning of artificial viscosity and interpreted it as a physical phenomenon based on the evolution of integral averages of the fluid solution (Navier-Stokes equations) over finite length scales ${}^{12}$ . His approach is reminiscent of von Neumann's work that made a distinction between average flow and the associated fluctuations to explain the entropy increase across a shock in adiabatic flow.

What separates the work of von Neumann and Grady from others is the idea that fluctuation is the dissipation, which is inherent to shock processes. Von Neumann emphasizes a viewpoint¹ that "oscillations must develop as a shock crosses a mass point, and they have a good physical significance and represent the thermal agitation caused by the degradation of energy through the shock." The idea of inherent oscillation across a shock is basically similar to that of Grady approach. But in the Grady approach fluctuation is defined as a random motion about the average and is only concerned with kinetic energy whereas von Neumann considered the variance of both kinetic and internal energies assuming harmonic motions for the fluctuations of field variables.

Experimental observations of fluctuations in shock waves are macroscopic and based primarily on interferometric measurements of particle velocity of emerging shock wave at rear sample surfaces. Additionally, existing data are mostly visual. So far as the authors are aware, the data on fluctuation have not been quantitatively analyzed and related to dissipated energies. Measurements of stress and density fluctuations do exist¹³, but they are rare and are again visual as the interferometric velocity data.

There are three-fold goals of this paper. First is to extend the von Neumann approach by removing the limitation of harmonic oscillations and calculating the variance (dissipation) based on an assumed probability distribution function. The key mathematical tool is the Luss formula for Taylor series expansion of a function of a random variable to calculate the variance ${}^{14}$. One advantage of the approach based on a probability distribution function is that it enables the investigation of fluctuations in terms of either underlying stochastic heterogenous meso-structures or shock-induced fluctuations in a probabilistic framework. Second is to examine the missing component of the Grady model that considered kinetic energy only. Third is application of the extended model to open a new perspective into looking at subjects such as shock-induced energy localization caused by underlying heterogeneous deformation mechanisms or heterogeneous microstructure or both. In this paper shock initiation of explosives that results in detonation is chosen as a test case. It is thought that probabilistic fluctuation-dissipation relations are well suited for modeling shock initiation of chemical reaction resulting from localized energy distribution. In the next section a brief description of the von Neumann and Grady models is provided for subsequent reference prior to the extension of the von Neumann model.

II. VON NEUMANN AND GRADY MODELS OF FLUCTUATION AND DISSIPATION

2.1 Von Neumann model

There are three key elements of the von Neumann model. First is the decomposition of internal energy into volumetric and entropic parts as shown in Eq. (1).

Then the total specific degraded energy ${\mathrm{\Delta }}_s$ (per unit mass) across a shock is given by Eq. (2).

where $p=-\frac{d{U}_{\ast }}{d\nu }$ and it is assumed for simplicity that $p_o=0$ . It is noted that $p$ is isentropic compression pressure. Second key-element is the discretized total energy of the point mass system that approximates the original hyperbolic conservation equations and is given in Eq. (3).

where $x_i$ denotes the coordinate of $i$-th point mass and $\nu ={\nu }_i=x_i-x_{i-1}$ . Third key-component is the interpretation of "oscillations" of two energies as "degraded energy", representing thermic agitation induced by crossing of a shock wave. He then assumes that the oscillations are harmonic and shows for the kinetic energy that

where $x_{osc}$ denotes the velocity amplitude of harmonic oscillations and $V_{osc}$ its velocity. Likewise, the degraded energy due to the oscillation of internal energy $U_{\ast }$ is given by Eq. (5). It is noted that flow is assumed to be adiabatic and S in Eq. (1) is constant.

where $c^2=\frac{dp}{d\rho }$ and $\rho =\frac{1}{\upsilon }$ . $c$ is local sound speed. So, the total degraded energy in Eq. (2) is now given by Eq. (6).

It is noted for later reference that Eq.(4) can be put in a form similar to that of Eq.(5). That is,

where $\bar{x}=0$ for the harmonic oscillations.

2.2 Grady model

Grady model is only concerned with the fluctuation of kinetic energy and starts with the decomposition of particle velocity into average field and its variant fluctuation as shown in Eq. (8) ${}^4$ .

Dissipated energy is related, as shown in Eq. (9), to the expected value of the square of the variant velocity on a plane (x) whose normal is collinear to the planar shock wave direction. The arrow denotes a vector quantity, but vectorial integrals were not discussed. Instead, integrated quantities such as standard deviation are deployed as shown in Eq. (9). That is, actual vectorial average calculation was not made explicit.

where $\sigma (t)$ is the standard deviation identified with kinetic dissipation. Secondly a fluctuation-dissipation relation shown in Eq.(10) is introduced as "a measure of the temporal correlation of kinetic dissipation and identified physically as "the acoustic phonon viscosity within the shock transition, which is said to be integral to the structuring of the time history of the wave through the shock transition."4

where $s=\frac{\rho x}{Z}$ and $Z=\rho c$ . Thus, variable $s$ is a normalized space coordinate $x$ scaled by the local sound speed $c$ . That is, $s=\frac{x}{c}$ . It is a measure of the domain of influence of sound propagation. He then introduces a quantity called frictional dissipation constant $\mathrm{\Gamma }$ within a time span of $\tau$ (shock rise time) such that

where $dm(=Zds)$ is said to be an element of areal mass. But there exists an ambiguity in the description of $\mathrm{\Gamma }$ . For example, it is found that $dm=Zds=Zd\left(\frac{\rho x}{Z}\right)$ . So, $Z(=\rho c)$ will not cancel unless it were constant during a span of shock transition time $\tau$ . So, the picture of $dm$ is not clear at least to the uninitiated to the analysis. Finally, an exponential variant velocity correlation function, Eq. (12), is introduced to perform the integration in Eq. (11).

\Gamma = $\frac{1}{2}$ \sigma^ {2} (\tau) Z \tau = E _ {k} (\tau) Z \tau \tag {13}

$\frac{\mathrm{\Gamma }}{Z}$ = E _ {k} (\tau) \tau = constant \tag {14}

However, the question of whether $Z(=\rho c)$ in Eq. (11) is constant within the time span $(\tau )$ of shock transition (or shock rise time) is an unanswered question. But, if $Z$ were assumed constant, then we can surmise upon comparison of Eq. (11) with Eq. (10) that

So, $D(\tau )$ is the quantity Grady calls "shock invariant." It may be of interest to close this section by noting that Eq. (14) can be generalized based on energy law observations of the comminution of solids as shown in Eq. (16)¹⁵.

where $\eta$ and $\beta$ are parameters representing a possible fractal dimension of time and a scaling factor respectively. Appearance of parameter $\beta$ suggests the existence of a critical time such that $\beta s={\tau }_c^{\eta }$ . Detail is referred to the paper¹⁵.

III. AN EXTENSION OF VON NEUMANN AND GRADY MODELS

Although specific details of the two models are different, the fundamental idea of relating fluctuation to dissipation is same in that shock dissipation is defined as the variance of field variables such as kinetic and internal energies. An illustration of their equivalence is given in Appendix. However, neither of the two models discusses probabilistic underpinning of the variance of field quantities. Therefore, the purpose of this work is to fill this gap using a distribution function that is amenable to analytic manipulation. In so doing we shall follow and generalize the von Neumann model and consider (1) both the kinetic and internal energies, considering their independence, and (2) show a consequence of treating only the kinetic energy fluctuation.

The chosen probability density function is a $\beta$ -function as shown in Eq. (17).

f (x) = C (x - a) ^ {\alpha} (b - x) ^ {\beta} \tag {17}

where $a$ and $b$ are fixed parameters, , and $C$ is a normalizing constant. Expectation (average) and variance of variable $x$ are given in Eqs. (18) and (19) respectively.

In a physical space as discussed in the next section, parameters $a$ and $b$ represent lower and upper limits of variable $x$ respectively. This means, for example, in the case of particle velocity $\mathbf{u}$ , dissipated energy due to fluctuation is proportional to ${\bar{u}}^2$ . Thus, it is of interest to compare it with the original artificial viscosity $q^{-}$ given in Eq. (20) ${}^2$ .

where $c_2$ is a dimensionless constant, $\rho$ is mass density, and $\mathrm{\Delta }x$ is the size of a computational cell. Since $\frac{du}{dx}$ can be considered constant over a small computational cell $\mathrm{\Delta }x$ , viscosity $q$ is proportional to the kinetic energy change $(\mathrm{\Delta }{u}^2)$ over the cell distance $\mathrm{\Delta }x$ . Thus, as done in the analysis of frictional force on the Brownian motion using Langevin equation, the artificial viscosity may be attributed to frictional force on the fluctuating mass motion $u$ . But the demonstration is beyond the scope of this paper. We shall now consider the kinetic and internal energy dissipations separately though use of Eqs. (17) - (19).

3.1 Kinetic energy fluctuation and dissipation

For the velocity fluctuation we shall assume a domain of particle-velocity magnitude to be finite and defined as $(0,u_{max})$ . Then the probability density distribution of particle-velocity magnitude is given in Eq. (21). The corresponding expectation and variant values are given in Eqns. (22) and (23) respectively.

One of the power coefficients (either $\alpha$ or $\beta$ ) and $u_{max}$ can be determined from the requirement that the maximum density distribution occurs at the average particle velocity $\bar{u}$ . That is, average state is considered to be the most probable state. Then,

Now, equating Eq. (24) to Eq. (22), it is found that $\alpha =\beta$ and that

The maximum particle velocity of $2\bar{u}$ is reasonable considering the fact that the reflection of a shock wave at solid free surface doubles the magnitude of particle velocity in a good approximation.

The determination of power coefficients $\alpha$ and $\beta$ requires experimental data that are not available at the moment. But a theoretical estimate can be made of their values by considering the overall dissipated energy that can be calculated based, for example, on a linear shock velocity-particle velocity relationship (i.e. $U=c_o+g\bar{u}$ ) or a thermodynamic argument. That is, it is known¹⁶ that for the linear shock-particle velocity model, the overall dissipated energy is given by Eq. (27).

Then using the linear $U-\bar{u}$ equation as a test case, Eq. (27) can be transformed into the dissipated energy in terms of the particle velocity to the same degree of approximation as shown in Eq. (28).

Now equating Eq. (27) and (28), and noticing that the dissipated kinetic energy based on the variance is only one half of the total dissipated energy in Eq. (28) (more later on this subject), one obtains

Rearranging the terms, one finds an estimate for $\alpha$ as given in Eq. (30).

Eq. (30) is valid only for $\bar{u}\le \frac{c_o}{g}$ within the confine of model assumptions and approximations. Physically, the limiting equality signifies the emergence of a uniform distribution at $\bar{u}=\frac{c_o}{g}$ , i.e., $\alpha =\beta =0$ and Eq. (19) = constant. If the dissipated kinetic energy is the only source to the total energy dissipation as assumed in the Grady analysis, then Eq. (30) needs to be modified to Eq. (31), indicating the emergence of the uniform distribution at a much lower loading.

At present there is no data to test Eq. (30) or Eq. (31). But it may be of interest to see the level of pressure needed to reach the uniform distribution of kinetic energy fluctuation. As a test case two materials are selected from the LASL Handbook¹⁷. They are Aluminum Alloy 6061 and PBX 9501 (unreacted) and have the following properties respectively.

where $U$ and $u$ are shock and particle velocity respectively and $u=\bar{u}$ . Then, the respective values of the threshold pressure are

These values are reasonable levels of shock pressure to expect the onset of uniform distribution of particle velocity fluctuation. It would be of interest to evaluate the prediction experimentally and numerically. They are well within the existing technical capability.

3.2 Internal energy fluctuation and dissipation

Same procedure as that used for the kinetic energy fluctuation will be followed for calculating the fluctuation and dissipation of internal energy. The difference is the calculation of the isentropic internal energy as defined by von Neumann. In this paper the isentropic energy will be approximated, following ${\mathrm{Grady}}^5$ , by integration of the Hugoniot curve $(p_h)$ based again on the linear U-u relation $(U=c_o+gu)$ as a test case. That is,

Then, using the Luss approximation, it is found that

where ${\sigma }_{\epsilon }^2$ is the variance of $\epsilon$ as defined by the density distribution function given in Eq. (34).

The average and variance are given in Eqs. (35) and (36) respectively.

Then the requirement that the internal energy dissipation is half of the total dissipated energy, yields

Now upon substitution of Eqs. (35) and (36) into Eq.(37), we find that