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). $$ U = U_{*}(v) + U_{*}(S) \tag{1} $$ Then the total specific degraded energy $\Delta_{s}$ (per unit mass) across a shock is given by Eq. (2). $$ \Delta_{s} = \Delta U_{*} = U_{*1} - U_{*0} = \frac {p _ {1}}{2} \left(v _ {1} - v _ {0}\right) - \int_ {v _ {0}} ^ {v _ {1}} p d v \tag{2} $$ where $p = -\frac{dU_{*}}{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). $$ \frac {1}{2} \sum_{i} \left(\frac {d x _ {i}}{d t}\right) ^ {2} + \sum_{i} U_{*}(x_{i} - x_{i-1}) \tag{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 $$ \overline {{\frac {1}{2} \left(V_{osc}\right) ^ {2}}} = \Delta_{s} (\text{kinetic}) = \frac {1}{4} \left(\dot{x}_{osc}\right) ^ {2} \tag{4} $$ 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_{*}$ is given by Eq. (5). It is noted that flow is assumed to be adiabatic and S in Eq. (1) is constant. $$ \Delta_{s} \left(U _ {*}\right) = \overline {{U _ {*} \left(x _ {i} - x _ {i - 1}\right)}} - U _ {*} \left(\overline {{x _ {i} - x _ {i - 1}}}\right) = \frac {1}{2} \frac {c ^ {2} \overline {{\left(\mathrm {v} _ {\text{osc}}\right) ^ {2}}}}{\mathrm {v} ^ {2}} = \frac {1}{4} \frac {c ^ {2} \overline {{\left(\mathrm {v} _ {\text{osc}} ^ {\text{am}}\right) ^ {2}}}}{\mathrm {v} ^ {2}} \tag{5} $$ 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). $$ \Delta_{s} (total) = \frac {1}{4} \left[ \frac {c ^ {2} \overline {{\left(\mathrm {v}_{\mathrm{osc}} ^ {a m}\right) ^ {2}}}}{\mathrm {u} ^ {2}} + \left(x_{osc} ^ {\cdot am}\right) ^ {2} \right] = \frac {p _ {1}}{2} \left(\mathrm {v} _ {1} - \mathrm {v} _ {0}\right) - \int_ {\mathrm {v} _ {0}} ^ {\mathrm {v} _ {1}} p d \mathrm {v} \tag{6} $$ It is noted for later reference that Eq.(4) can be put in a form similar to that of Eq.(5). That is, $$ \Delta_{s} (\text{kinetic}) = \overline {{\frac {1}{2} \dot{x}^{2}}} - \frac {1}{2} \left(\bar {x}\right) ^ {2} \tag{7} $$ where $\overline{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$ . $$ V (x, t) = \mu (t) + \vec {u} (x, t) \tag {8} $$ 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. $$ \text {Kinetic dissipation} = E _ {k} (t) = \langle u ^ {2} (x, t) \rangle / 2 = \sigma^ {2} (t) / 2 \tag {9} $$ 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 $$ D (t) = \int_ {0} ^ {\infty} \langle \vec {u} (0, t) \vec {u} (s, t) d s = \int_ {0} ^ {\infty} K (s, t) d s \tag {10} $$ 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 $\Gamma$ within a time span of $\tau$ (shock rise time) such that $$ \Gamma = \frac {1}{2} \int_ {0} ^ {\infty} \left\langle \vec {u} (0, \tau) \vec {u} (s, \tau) \right\rangle d m = \frac {1}{2} \int_ {0} ^ {\infty} K (s, \tau) Z d s \tag {11} $$ where $dm(= Zds)$ is said to be an element of areal mass. But there exists an ambiguity in the description of $\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). $$ K (s, t) = \langle \vec{u} (0 , \tau) ^ {2} \rangle e ^ {- s / \tau} = \sigma^ {2} (\tau) e ^ {- s / \tau} \tag {12} Last step in the Grady analysis is the integration of Eq.(11) with the assumption of constant $Z$ . The result is shown in Eq. (13). $$ \Gamma = \frac {1}{2} \sigma^ {2} (\tau) Z \tau = E _ {k} (\tau) Z \tau \tag {13} $$ Since $Z$ is assumed constant over the time span of $\tau$ (i.e. shock transition), Eq. (13) expresses an energy-time criterion called shock invariant at the terminal Hugoniot state. That is, $$ \frac {\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 $$ D (\tau) = \int_ {0} ^ {\infty} K (s, \tau) d s = \frac {\Gamma}{Z} = E _ {k} (\tau) \tau \tag {15} $$ 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)¹⁵. $$ \frac {\Gamma}{Z} = \frac {1}{\beta} E _ {k} (\tau) \tau^ {\eta} \tag {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, $\alpha > -1$ , $\beta > -1$ and $C$ is a normalizing constant. Expectation (average) and variance of variable $x$ are given in Eqs. (18) and (19) respectively. $$E [ x ] = a + \frac {\alpha + 1}{\alpha + \beta + 2} (b - a) \tag {18} $$ $$ \operatorname {v a r} [ x ] = E \left[ x ^ {2} \right] - (E [ x ]) ^ {2} = \frac {(b - a) ^ {2} (\alpha + 1) (\beta + 1)}{(\alpha + \beta + 2) ^ {2} (\alpha + \beta + 3)} \tag {19} $$ 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 $\overline{u}^2$ . Thus, it is of interest to compare it with the original artificial viscosity $q^{-}$ given in Eq. (20) $^2$ . $$ q = c _ {2} \rho \left| \frac {d u}{d x} \right| ^ {2} \Delta x ^ {2} \tag {20} $$ where $c_{2}$ is a dimensionless constant, $\rho$ is mass density, and $\Delta x$ is the size of a computational cell. Since $\frac{du}{dx}$ can be considered constant over a small computational cell $\Delta x$ , viscosity $q$ is proportional to the kinetic energy change $(\Delta u^{2})$ over the cell distance $\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 $\left(0, u_{\max}\right)$ . 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. $$ f (u) = C (u) ^ {\alpha} \left(u _ {\max } - u\right) ^ {\beta} \tag {21} $$ $$ E [ u ] = \bar {u} = \frac {\alpha + 1}{\alpha + \beta + 2} u _ {\text {m a x}} \tag {22} $$ $$ \operatorname {v a r} [ u ] = \frac {(\alpha + 1) (\beta + 1)}{(\alpha + \beta + 3) (\alpha + \beta + 2) ^ {2}} u _ {\max } ^ {2} \tag {23} $$ 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 $\overline{u}$ . That is, average state is considered to be the most probable state. Then, $$ u (\frac{\partial f}{\partial u} = 0) = \frac {\alpha}{(\alpha + \beta)} u _ {m a x} \tag {24} $$ Now, equating Eq. (24) to Eq. (22), it is found that $\alpha = \beta$ and that $$ \bar {u} = \frac {1}{2} u _ {\text {m a x}} \tag {25} $$ $$ v a r [ u ] = \frac {1}{(2 \alpha + 3)} (\bar {u}) ^ {2} \tag {26} $$ The maximum particle velocity of $2\overline{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\overline{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). $$ \Delta_ {s} (t o t a l p e r u n i t m a s s) = \frac {1}{3} g c _ {0} ^ {2} \bar {\varepsilon} ^ {3} + \dots \tag {27} $$ $$ {where} \overline {{\varepsilon}} = 1 - \frac {\overline {{v}}}{v _ {o}} = 1 - \frac {\rho_ {o}}{\rho} = \frac {\overline {{v}}}{v _ {o}}. $$ Then using the linear $U - \overline{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). $$ \Delta_ {s} (\text {total}, \text {unit mass}) = \frac {2 g}{3} \frac {\bar {u}}{c _ {o}} \left(\frac {\bar {u} ^ {2}}{2}\right) + \dots \tag {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 $$ \frac {1}{(2 \alpha + 3)} \frac {\left. \bar {u}\right) ^ {2}}{2} = \frac {1}{2} \left[ \frac {2 g}{3} \frac {\bar {u}}{c _ {o}} \left(\frac {\bar {u} ^ {2}}{2}\right) \right] \tag {29} $$ Rearranging the terms, one finds an estimate for $\alpha$ as given in Eq. (30). $$ \alpha = \frac {3}{2} \left(\frac {c _ {o}}{g \bar {u}} - 1\right) \tag {30} $$ Eq. (30) is valid only for $\overline{u} \leq \frac{c_o}{g}$ within the confine of model assumptions and approximations. Physically, the limiting equality signifies the emergence of a uniform distribution at $\overline{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. $$ \alpha = \frac {3}{2} \left(\frac {c _ {o}}{2 g \bar {u}} - 1\right) \tag {31} $$ 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. $$ Al 6061: \rho_{0}=2710 \frac{kg}{m^{3}}; U = 5.35 + 1.34 u \left(\frac{km}{sec}\right) $$ $$ PBX9501: \rho_{0}=1830 \frac{kg}{m^{3}}; U = 2.95 + 1.50 u \left(\frac{km}{sec}\right) $$ where $U$ and $u$ are shock and particle velocity respectively and $u = \overline{u}$ . Then, the respective values of the threshold pressure are $$ P (Al6061) = 116 \text{GPa} \text{and} P (PBX9501) = 21.2 \text{GPa} $$ 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, $$ U _ {*} (\varepsilon) (p e r u n i t v o l u m e) \simeq \int_ {0} ^ {\varepsilon} p _ {h} d \varepsilon = \int_ {0} ^ {\varepsilon} \frac {\rho_ {o} c _ {o} ^ {2} \varepsilon}{(1 - g \varepsilon) ^ {2}} d \varepsilon = \frac {1}{2} \rho_ {o} c _ {o} ^ {2} \varepsilon^ {2} + \frac {2}{3} \rho_ {o} c _ {o} ^ {2} g \varepsilon^ {3} + \dots \tag {32} $$ Then, using the Luss approximation, it is found that $$ \Delta_ {s} (i n t e r n a l e n e r g y) = \overline {{U _ {*} (\varepsilon)}} - U _ {*} (\overline {{\varepsilon}}) = \frac {1}{2} \sigma_ {\varepsilon} ^ {2} c _ {o} ^ {2} (1 + 4 g \overline {{\varepsilon}} + \dots) \tag {33} $$ where $\sigma_{\varepsilon}^{2}$ is the variance of $\varepsilon$ as defined by the density distribution function given in Eq. (34). $$ f (\varepsilon) = K (\varepsilon) ^ {a} \left(\varepsilon_ {\max } - \varepsilon\right) ^ {b} \tag {34} $$ The average and variance are given in Eqs. (35) and (36) respectively. $$ E [ \varepsilon ] = \bar {\varepsilon} = \frac {a + 1}{a + b + 2} \varepsilon_ {m a x} \tag {35} $$ $$ \operatorname {v a r} [ \varepsilon ] = \frac {(a + 1) (b + 1)}{(a + b + 3) (a b + 2) ^ {2}} \varepsilon_ {\text {m a x}} ^ {2} \tag {36} $$ Then the requirement that the internal energy dissipation is half of the total dissipated energy, yields $$ \sigma_ {\varepsilon} ^ {2} = \frac {1}{3} g \bar {\varepsilon} ^ {3} \tag {37} $$ Now upon substitution of Eqs. (35) and (36) into Eq.(37), we find that $$ \varepsilon_ {m a x} = \frac {(b + 1) (a + b + 2)}{(a + b + 3) (a + 1) ^ {2}} \left(\frac {3}{g}\right) \tag {38} $$ Similarly, the requirement that the maximum density distribution is located at the average strain yields the results that are similar to those for the particle velocity and are given in Eqs. (39) and (40). $$ \varepsilon_ {m a x} = \frac {a + b}{a} \bar {\varepsilon} \tag {39} $$ $$ \bar {\varepsilon} = \bar {u} / U = \bar {u} / \left(c _ {o} + g \bar {u}\right), \tag {40} $$ Finally, combining Eqs. (38)-(40), one finds $$ a (= b) = \frac {3}{2} \left(\frac {1}{g \bar {\varepsilon}} - 1\right) \tag {41} $$ Again, uniform distribution of strain fluctuation emerges when the average strain $\bar{\varepsilon}$ reaches the value of $1 / g$ . A representative value of $g$ in the LASL Handbook is about 1.5. Then Eq. (41) yields a very high relative strain of $\frac{2}{3}$ for the uniform distribution of strain fluctuation. Again, no data are available to evaluate the threshold value, but it is a reasonably high value for the appearance of a uniform distribution. Since $\overline{\varepsilon}$ is related to $\overline{u}$ by the definition: $\frac{\overline{\varepsilon} = \overline{u}}{U = \overline{u}} / (c_{o} + g\overline{u})$ , coefficients $a$ and $\alpha$ are interrelated by Eq. (42). $$ a = \alpha + \frac {3}{2} \tag {42} $$ At present no data are available to evaluate the prediction. But Eq. (42) will be useful to estimate the coefficient $\alpha$ from $\alpha$ that can be determined by interferometric particle velocity measurements. # IV. AN APPLICATION OF PROBABILISTIC FLUCTUATION-DISSIPATION RELATIONS It is shown in the previous section that with a minimum number of input data fluctuation-dissipation relations for hydrodynamic shock compression can be developed based on an underlying density distribution function of either particle velocity or relative compressive strain or both. An interesting result is the appearance of an upper limit of the dissipated energy for materials exhibiting the linear U-u relation as given in Eq. (43), which also corresponds to the emergence of a uniformly distributed fluctuation and is obtained by combining Eqs. (28) and (29). $$ \max . \Delta_ {s} (\text {unit mass}) = \frac {1}{3} \left(\frac {- 2}{2}\right) \tag {43} $$ Understanding the meaning of this maximum dissipation fraction for the materials obeying the linear U-u relation is left as a future problem. The goal of this section is to explore a possible application of the fluctuation-dissipation relations based on the beta density distribution function. A chosen test case is the initiation of detonation for which Eqs. (28) and (33) are interpreted to represent the sum of distributed "hot spots" (high dissipated-energy spots). For example, since $\varepsilon_{max} = 2\bar{\varepsilon}$ , max. $\Delta_s = 8\Delta_s(\bar{\varepsilon})$ . Then, if we use an approximation for the temperature change through use of the equation: $C_v\Delta T = \Delta_s$ where $C_v$ is the heat capacity, corresponding maximum hot-spot temperature will be eight times the average Hugoniot temperature. This magnification is definitely sufficient to initiate localized chemical reactions even at relatively low shock pressures. For a specific illustration of application of the distributed internal energy, a simple, but commonly used reaction model is chosen to calculate reaction-completion time as preparation for the next section where reaction time is reinterpreted probabilistically by use of the fluctuation-dissipation relation. It is noted that reaction-completion time is thought to be able to represent several situations. It could be the time necessary to form critical initiation nuclei, time to detonation, or rise time of detonation wave front, depending on the interpretation and form of the kinetic equation $^{18,19}$ . In this paper reaction-completion time is considered to be the time-to-detonation in Pop-plots $^{17}$ . # 4.1 Reaction-completion time and modeling threshold conditions Goal of this section is to calculate the resection-completion time based on the kinetic equation given in Eq. (44) as preparation for application of the fluctuation-dissipation relations. This equation is an analog of standard temperature-based reaction time calculation $^{19}$ and resembles the use of variable $Z$ in the CREST model for which $Z$ denotes the entropy-dependent part of internal energy function ( $e$ ) similar to the von Neumann decomposition in Eq.(1). That is, $e = e_0(\upsilon) + f(\upsilon)Z(S)$ where $S$ is entropy $^{20}$ . $$ \dot {x} = v \exp \left(- \frac {E _ {a}}{E _ {o} + \Delta_ {s} + x \Delta q}\right) \tag {44} $$ where $x$ is fraction of mass reacted, $x$ its time rate of change, $\nu$ frequency factor, $\Delta q$ energy release due to reaction, $E_{a}$ activation energy, $E_{o}$ isentropic compression energy at inert Hugoniot state. In integrating Eq. (44), we assume that the release of $\Delta q$ has no influence on $E_{o}$ and $\Delta_{s}$ . Then, the time to completion $(\tau)$ is given in Eq. (45). $$ \tau = \frac {1}{v} \frac {\left(E _ {a} + \Delta_ {s}\right) ^ {2}}{E _ {a} \Delta q} \exp \left(\frac {E _ {a}}{E _ {o} + \Delta_ {s}}\right) \tag {45} $$ There are many phenomenological models of threshold conditions such as shock initiation and Pop-plots that make use of the time to reaction-completion $^{19,21}$ . In this paper a commonly used power function given in Eq.(46) is chosen to illustrate an application of Eq.(45). $$ P ^ {n} \tau = \text {constant} \tag {46} $$ where $P$ is sustained inert shock pressure. Then, the power coefficient $n$ can be evaluated and compared with experimental values through use of Eq. (47). $$ n = - \frac {P}{\tau} \frac {d \tau}{d P} \tag {47} $$ In order to calculate Eq. (47) by use of using Eq. (45), some approximations that are commensurate to the level of approximation for calculating dissipated energy in Eq. (27), are introduced. These approximations are $$ u = \frac {c _ {o} \varepsilon}{1 - g \varepsilon} \simeq c _ {o} \varepsilon (1 + g \varepsilon) \tag {48} $$ $$ P = \rho_ {o} U u \simeq \rho_ {o} c _ {o} ^ {2} \varepsilon \left(1 + 2 g \varepsilon + g ^ {2} \varepsilon^ {2} + \dots\right) \simeq \rho_ {o} c _ {o} ^ {2} \varepsilon \tag {49} $$ Then upon combining Eq. (49) with Eq. (27), it is found that $$ P = \left(\frac {3 \rho_ {o} ^ {3} c _ {o} ^ {4}}{g}\right) ^ {1 / 3} \Delta_ {s} ^ {1 / 3} \tag {50} $$ Finally, it may be shown after some algebra that $$ n = - \frac {d l n \tau}{d l n P} = - \frac {P}{\tau} \frac {d \tau / d \Delta_ {s}}{d P / d \Delta_ {s}} = \frac {3 \Delta_ {s}}{\left(E _ {o} + \Delta_ {s}\right)} \left[ \frac {E _ {a}}{\left(E _ {o} + \Delta_ {s}\right)} - 2 \right] = \frac {3 \Delta_ {s}}{E _ {h}} \left[ \frac {E _ {a}}{E _ {h}} - 2 \right] \tag {51} $$ where it is noted that the sum of the isentropic energy and the dissipated energy is approximately the Hugoniot energy even when the former is approximated by use of the Hugoniot curve. Also, Eq. (51) can be replaced by Eq. (52) where the equipartition of energy is used at the Hugoniot state, i.e., $E_{h} = \frac{u^{2}}{2}$ $$ n = \frac {6 \Delta_ {s}}{u ^ {2}} \left(\frac {2 E _ {a}}{u ^ {2}} - 2\right) \simeq 4 g \left(\frac {u}{c _ {o}}\right) \left(\frac {E _ {a}}{u ^ {2}} - 1\right) \tag {52} $$ where it is used that $\varepsilon = \frac{u}{U} \simeq \frac{u}{c_o}$ to the first order, commensurate to the approximation used to get Eq. (27). As expected, $n$ is a function of particle velocity and depends on factors that include the approximations involved. Thus, after a spot check on Eq. (52) in the next section, the effect of particle velocity distribution on $n$ is examined in the succeeding section in stipulation that experimentally measured values of $n$ correspond better to average $n$ than specific functional values based on Eq. (52). Again, data from LASL Explosive Handbook¹⁷ are used for the spot check on the function in Eq. (52). # 4.2 A spot-check on Eq.(52) Since LASL data specify the power coefficients for select ranges of impact particle velocity, Eq.(52) may be best checked by inverting it to get the velocity for any given $n$ as shown in Eq. (53), which has a single positive root given in Eq. (54) $$ u ^ {2} + \left(\frac {c _ {o}}{4 g} n\right) u - E _ {a} = 0 \tag {53} $$ $$ u = \frac {1}{2} \left(\frac {c _ {o}}{4 g} n\right) \left[ - 1 + \sqrt {1 + 4 E _ {a} \left(\frac {4 g}{n c _ {o}}\right) ^ {2}} \right] \tag {54} $$ In testing Eq. (54) with LASL data, the following conversion formula is used: $1 \text{joule} = kg(m/sec)^2$ . Four materials are chosen arbitrarily. They are PETN, TATB, TNT, and PBX 9501. Experimental data are grouped by the initial density. # PETN LASL data: $$ E _ {a} = 4 7. 0 kcal = 0. 7 8 9 \left(\frac {k m}{s e c}\right) ^ {2}; \text {molecular weight} = 3 1 6. 1 5 $$ (a) $\rho_{o} = 1.72\frac{g}{cc}; n = 2.0; U = 1.83 + 3.45u\left(\frac{km}{sec}\right)$ ; $$ 0. 1 < u < 0. 7 \left(\frac {k m}{s e c}\right); 1. 7 < P < 3. 9 (G P a) $$ (b) $\rho_{o} = 1.77\frac{g}{cc}; n = 4.545; U = 2.87 + 1.69u\left(\frac{km}{sec}\right)$ ; $$ 0. 5 < u < 1. 5 \left(\frac {k m}{s e c}\right); 1. 7 < P < 2. 5 4 (G P a) $$ where the range of $P$ is for Pop-plot measurements. A scatter is obvious, but the substitution of the data into Eq. (54) yields (a) $u(n = 2.0) = 0.765\frac{km}{sec}; P = 5.88GPa$ (b) $u(n = 4.55) = 0.347\frac{km}{sec}; P = 2.12(GPA)$ Predicted value of Pop-plot pressure for case (a) is outside of the experimental range. But the second prediction is close to the median value of the experimental range. At present it is difficult to interpret the discrepancy in case (a) due to the scatter in the experimental data and their sensitivity to minor differences in microstructure and processing parameters. TATB (one matching density for shock and Pop-plot) LASL data: $$ E _ {a} = 5 9. 9 kcal = 0. 9 7 0 1 \left(\frac {k m}{s e c}\right) ^ {2}; molecular number = 2 5 8. 1 8 $$ (a) $\rho_{o} = 1.876\frac{g}{cc}; n = 2.778; U = 2.037 + 2.497u\left(\frac{km}{sec}\right)$ ; $$ 0. 4 8 < u < 1. 5 4 \left(\frac {k m}{s e c}\right), 3. 2 7 < P < 5. 6 4 (G P a) $$ (b) $\rho_{o} = 1.876g / cc$ ; $n = 2.778$ ; $U = 1.663 + 2.827u\left(\frac{km}{sec}\right)$ ; $$ 0. 4 8 < u < 1. 5 4, 5. 9 3 < P < 1 6. 5 (G P a) $$ (c) $\rho_{o} = 1.876g / cc$ ; $n = 2.778$ ; $U = 1.46 + 3.68u\left(\frac{km}{sec}\right)$ ; $$ 0. 0 < u < 0. 4 8 \left(\frac {k m}{s e c}\right), 1 1. 4 < P < 1 6. 2 2 (G P a) $$ Data yield predicted values as follows. (a) $u(n = 2.778) = 0.7419\left(\frac{km}{sec}\right)$ ; $P = 5.413$ (GPa) (b) $u(n = 2.778) = 0.7996\left(\frac{km}{sec}\right)$ ; $P = 5.89$ (GPa) (c) $u(n = 2.778) = 1.188\left(\frac{km}{sec}\right)$ ; $P = 13.0$ (GPa) Predicted values fall within the ranges of Pop-plot velocity and pressure measurements except case (b). But its deviation is only about $0.6\%$ . # TNT LASL Data: $$ E _ {a} = 3 4. 4 \text {kcal} = 0. 6 3 3 7 \left(\frac {\text {k m}}{\text {s e c}}\right) ^ {2}; \text {m o l a r w e i g h t} = 2 2 7. 1 3. $$ $$ \rho_ {o} = 1. 6 2 - 1. 6 3 4 g / c c; n = 3. 2 2; U = 2. 5 7 + 1. 8 8 u \left(\frac {k m}{s e c}\right); 0. 0 < u < 2. 0 $$ $$ 0. 0 < u < 2. 0 \frac {(k m}{s e c)}, 9. 1 7 < P < 1 7. 1 (G P a) $$ Corresponding predicted values are $$ u (n = 3. 2 2) = 0. 4 1 8 \frac {(k m}{s e c}), P = 2. 2 9 (G P a) $$ The value of pressure is well below the experimental pressure range and is an exception among the four cases examined. PBX 9501 LASL data: $$ \begin{array}{l} E _ {a} = 4 0. 1 kcal = 0. 5 8 5 9 \left(\frac {k m}{s e c}\right) ^ {2}; molecular weight = 2 8 6. 3 6 \\ (a) \rho_ {o} = 1. 8 3 3 g / c c; n = 1. 8 8 7; U = 2. 5 0 1 + 2. 2 6 1 u \left(\frac {k m}{s e c}\right); \\ 0. 0 7 < u < 0. 8 \left(\frac {k m}{s e c}\right), 2. 3 8 < P (G P a) < 7. 3 2 \\ (b) \rho_ {o} = 1. 8 4 4 g / c c; n = 2. 2 2 2; U = 2. 9 5 3 + 1. 5 0 7 u \left(\frac {k m}{s e c}\right); \\ 0. 1 < u < 0. 9 \left(\frac {k m}{s e c}\right), 2. 4 7 < P < 7. 2 1 (G P a) \\ \end{array} $$ Although initial densities are slightly different ( $\sim 0.6\%$ ), two cases are treated as the same material. That is, the difference in cases (a) and (b) is viewed as a scatter for the same material. Then upon substitution of relevant numbers, predicted velocities and pressures are given as follows. $$ \begin{array}{l} \mathrm {(a)} u (n = 1. 8 8 7) = 0. 5 4 4 \left(\frac {k m}{s e c}\right); P = 3. 7 3 (G P a) \\ \left(\mathrm {b}\right) u (n = 2. 2 2 2) = 0. 3 9 5 \left(\frac {k m}{s e c}\right); P = 2. 5 8 (\mathrm {G P a}) \\ \end{array} $$ These values are well within the experimental ranges of Pop-plot measurements. Overall, it is concluded that the results of the above-described spot check indicate that the experimentally measured power coefficients are best modeled by considering distributed fluctuations of field variables. That is, measured values of $n$ are averaged $n$ over the particle velocity distribution. # 4.3. Influence of distributed velocity field on $n$ Test calculations in section 4.2 where it is seen that for the select materials the values of $n$ were found to be within or close to the experimental ranges (with the exception of TNT), gave a hint that averaged $n$ over a distributed velocity field would offer a more realistic estimate of experimentally measured values. Therefore, by combining Eqs. (21) and (52), this average is defined as shown in Eq. (55). $$ \bar {n} = k \int_ {\text {lower bound}} ^ {\text {upper bound}} 4 g \left(\frac {u}{c _ {o}}\right) \left(\frac {E _ {a}}{u ^ {2}} - 1\right) (u) ^ {\alpha} \left(u _ {\max } - u\right) ^ {\beta} d u \tag {55} $$ where $k$ is the normalizing constant for the distribution function. Since Eq. (55) is not easily amenable to get general analytical solutions, a special case is considered using the previously considered distribution function for which $\alpha = \beta = 1$ . But here they are not tied to loading strength and considered to be an assumed function. Then, it is found that $$ \bar {n} = \frac {4 k g}{c _ {o}} \int_ {u _ {\min }} ^ {u _ {\max }} \left(E _ {a} - u ^ {2}\right) \left(u _ {\max } - u\right) d u \tag {56} $$ where $k = \frac{6}{u_{max}^3}$ for the distribution range of $(0, u_{max})$ . Also, it may be seen now that since $0 < n$ , $u_{max} = \sqrt{E_a}$ . Then due to the symmetry of the distribution function, $u_{min}$ is $\sqrt{E_a} / 2$ , which is the average particle velocity, indicating that averaging must be performed over the velocity range where dissipated energy exceeds the average value. Then, the integral of Eq. (56) yields the desired average $n (= \overline{n})$ as shown in eq.(57). $$ \bar {n} = \frac {g}{c _ {o}} \sqrt {E _ {a}} \left(1 0 - \frac {1 6 0 8}{1 9 2}\right) \simeq \frac {1 . 6 3 g}{c _ {o}} \sqrt {E _ {a}} \tag {57} $$ Again, testing of this equation is done using the samples used earlier for the spot check. The results are listed below. The values of $\overline{n}$ in the parentheses are experimental values in the Handbook for comparison. $$ \begin{array}{l} \underline {{\mathrm {P E T N}}} \cdot E _ {a} = 0. 7 8 9 \left(\frac {k m}{s e c}\right) ^ {2}, c _ {o} = 1. 8 3 \left(\frac {k m}{s e c}\right), g = 3. 4 5. \overline {{n}} = 2. 7 3 (2. 0) \\ \underline {{\mathrm {T A T B}}} \cdot E _ {a} = 0. 9 7 0 1 \left(\frac {k m}{s e c}\right) ^ {2}, c _ {o} = 2. 0 3 7 \left(\frac {k m}{s e c}\right), g = 2. 4 9 7. \overline {{n}} = 1. 9 7 (2. 7 7) \\ \underline {{\mathrm {T N T}}} \cdot E _ {a} = 0. 6 6 3 \left(\frac {k m}{s e c}\right) ^ {2}, c _ {o} = 2. 2 5 7 \left(\frac {k m}{s e c}\right), g = 1. 8 8. \overline {{n}} = 0. 9 6 4 (3. 2 2) \\ \underline {{\mathrm {P B X}}} 9 5 0 1. E _ {a} = 0. 5 8 5 9 \left(\frac {k m}{s e c}\right) ^ {2}, c _ {o} = 2. 9 5 3 \left(\frac {k m}{s e c}\right), g = 2. 2 6 1. \overline {{n}} = 1. 1 3 (1. 8 8 7) \\ \end{array} $$ With the exception of TNT again, predictions are in a reasonable agreement with the experimental values, considering the approximations made in the analysis and the use of an untested, hypothetical distribution function. They can be improved upon once the statistical data become available. At least closeness of the predictions to experimental values indicates a possible new approach to look at Pop-plot data based on the probabilistic fluctuation-dissipation relations. # V. CONCLUSIONS Fluctuation-Dissipation (F-D) relations proposed by J. von Neumann and D. Grady for shock wave propagation in hydrodynamic solids were extended by introducing underlying distribution functions for both kinetic and internal energies. Due to the lack of experimental data, a beta function is hypothesized for analytical convenience to develop and evaluate the probabilistic F-D relations. One application of the probabilistic F-D relations is investigated in prediction of Pop-plot power coefficients. Predicted coefficient values were found to be in a reasonable agreement with experimentally measured values. The results are thought to open a new perspective to analyze stochastic heterogeneous shock-induced energy localization and its significance through use of probabilistic F-D relations. # ACKNOWLEDGMENTS First author (Y.H.) wishes to acknowledge (1) partial support of Sandia National Laboratories (SNL) under Contract SPO #2432815-4 for the work reported in this paper and (2) help of Ms. C. Rountree (librarian at the Air Force Research Lab. Eglin AFB, Florida) in literature search. SNL is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Security Administration under Contract No. DE-NA0003525. # APPENDIX. EQUIVALENCE OF VON NEUMANN AND GRADY MODELS Von Neumann defined degraded energy as the variance of both K.E. and I.E. For example, if we let $v$ to be the particle velocity, then the degraded kinetic energy is given by Eq. (58). $$ \text {Degraded} \mathrm {K.E.} = \left\langle v ^ {2} \right\rangle - \left\langle v \right\rangle^ {2} \tag {58} $$ If we decompose $v$ as Grady did in such a way that $v = u + \delta$ , then von Neumann's degraded energy is given by Eq. (59). $$ Degraded K. E. = \left\langle (u + \delta) ^ {2} \right\rangle - \left\langle u + \delta \right\rangle^ {2} \tag {59} $$ Expanding the bracket and assuming that $\langle \delta \rangle = 0$ , $\langle u^2\rangle = \langle u\rangle^2$ , and $u$ and $\delta$ are independent, we obtain Eq. (60) that shows the equivalence of the two models. $$ \text {Degraded} \quad \mathrm {K.E.} = \langle v ^ {2} \rangle - \langle v \rangle^ {2} = \langle u ^ {2} \rangle - 2 \langle u \rangle \langle \delta \rangle + \langle \delta^ {2} \rangle - \langle u \rangle^ {2} = \langle \delta^ {2} \rangle \tag {60} $$

Generating HTML Viewer...