Shear shock evolution in incompressible soft solids

TL;DR

This paper investigates how shear waves evolve into shocks in incompressible soft solids, deriving formulas and phase maps to predict shock formation and applying findings to brain tissue to inform injury prevention.

Contribution

It provides closed-form expressions for shock formation distance and non-dimensional phase maps, revealing the influence of loading and material parameters on shear shock evolution.

Findings

Shock formation is sensitive to amplitude, shape, and ramp time of loading.

Smaller brains are less susceptible to shear shocks under realistic loadings.

Non-dimensional maps can guide protective structure design to prevent shocks.

Abstract

Nonlinear evolution of shear waves into shocks in incompressible elastic materials is investigated using the framework of large deformation elastodynamics, for a family of loadings and commonly used hyperelastic material models. Closed form expressions for the shock formation distance are derived and used to construct non-dimensional phase maps that determine regimes in which a shock can be realized. These maps reveal the sensitivity of shock evolution to the amplitude, shape, and ramp time of the loading, and to the elastic material parameters. In light of a recent study (Espindola et al., 2017), which hypothesizes that shear shock formation could play a signicant role in Traumatic Brain Injury (TBI), application to brain tissue is considered and it is shown that the size matters in TBI research. Namely, for realistic loadings, smaller brains are less susceptible to formation of shear shocks. Furthermore, given the observed sensitivity to the imparted waveform and the constitutive properties, it is suggested that the non-dimensional maps can guide the design of protective structures by determining the combination of loading parameters, material dimensions, and elastic properties that can avoid shock formation.

Equations

$$x_1=X_1+u(X_2,t),x_2=X_2,x_3=X_3$$

$$\gamma (X_2,t)=\text{partialderivative}{x}_1{X}_2=\text{partialderivative}u(X_2,t)X_2$$

$$v(X_2,t)=\text{partialderivative}u(X_2,t)t\text{where}v(0,t)=V(t)$$

$$\mathbf{F}=\text{partialderivative}\mathbf{x}\mathbf{X}=\left[\begin{array}{ccc}1& \gamma & 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],\mathbf{B}=\mathbf{F}\cdot {\mathbf{F}}^T=\left[\begin{array}{ccc}{\gamma }^2+1& \gamma & 0\\ \gamma & 1& 0\\ 0& 0& 1\end{array}\right]$$

$$W=W(I_1,I_2)$$

$$\mathbf{P}=\mathbf{T}\mathbf{F}^{-T}\quad,\quad{\mathbf{T}}=-p~{}{\mathbf{{\mathit{{\mathbf{1}}}}}}+2\big{[}(W_{1}+I_{1}W_{2})\mathbf{B}-W_{2}\mathbf{B}^{2}\big{]}$$

$$\mathbf{P}(X_2,t)=\left[\begin{array}{ccc}2W_1+4W_2-p& 2\gamma (W_1+W_2)& 0\\ \gamma (-2W_2+p)& 2W_1+4W_2-p& 0\\ 0& 0& 2W_1+4W_2+2{\gamma }^2{W}_2-p\end{array}\right]$$

$${\int }_{{\mathcal{D}}_0}{\rho }_0\dot{\mathbf{v}}\mathrm{d}{\mathrm{V}}_0={\int }_{{\partial \mathcal{D}}_0}\mathbf{P}{\mathbf{n}}_0\mathrm{d}{\mathrm{A}}_0$$

$$\mathrm{Div}\mathbf{P}={\rho }_0\dot{\mathbf{v}}$$

$$\text{partialderivative}\tau X_2$$

$$\text{partialderivative}{P}_{22}{X}_2$$

$$\text{partialderivative}[2]ut=c^2\text{partialderivative}[2]u{X}_2\text{where}c(\gamma )=\sqrt{\frac{1}{{\rho }_0}\text{partialderivative}\tau \gamma }$$

$$c_s=\underset{\gamma \to 0}{\lim }c(\gamma )=\sqrt{\frac{\mu }{{\rho }_0}}$$

$$\text{partialderivative}vt$$

$$\text{partialderivative}\gamma t$$

$${\tau }^{\prime \prime }(|\gamma |)\frac{\partial |\gamma |}{\partial X_2}<0$$

$$Q(\gamma )=\int {\text{displaylimits}}_0^{\gamma }c(\gamma )\text{d}\gamma ,\frac{\partial Q}{\partial t}=c(\gamma )\frac{\partial \gamma }{\partial t},\frac{\partial Q}{\partial X_2}=c(\gamma )\frac{\partial \gamma }{\partial X_2}$$

$$v+Q=f(\beta )$$

$$\text{where}\text{differential}\beta =0\text{ along }\text{derivative}{X}_{2}t=-c$$

$$v(X_2\ge 0,t=0)=0,\gamma (X_2\ge 0,t=0)=0,v(X_2=0,t>0)=V(t)$$

$$\hat{v}(\alpha )$$

$$\hat{Q}(\alpha )$$

$$X_2=\hat{c}(\alpha )(t-\alpha )$$

$$\alpha =\alpha (X_2,t)$$

$$t_{\infty }^{\alpha }=\alpha +\frac{\hat{c}}{{\hat{c}}^{\prime }},t_{\infty }^{\alpha }>\alpha$$

$$X_{\infty }^{\alpha }=\hat{c}(t_{\infty }^{\alpha }-\alpha )=\frac{{\hat{c}}^2}{{\hat{c}}^{\prime }}$$

$$V^{\prime }(\alpha )=-c\text{derivative}\hat{\gamma }\alpha \text{i.e}\text{derivative}Vt=-c\text{derivative}\gamma (0,t)t$$

$$t_{\infty }^{\alpha }$$

$$X_{\infty }^{\alpha }$$

$$v_{X_2}^{\alpha }=\text{partialderivative}\hat{v}(\alpha (X_2,t))X_2={\hat{v}}^{\prime }\text{partialderivative}\alpha X_2,v_{X_2}^{\alpha }(X_2)=\frac{V^{\prime }(\alpha )}{{\hat{c}}^{\prime }{X}_2/\hat{c}-\hat{c}},{\tilde{v}}_{X_2}^{\alpha }(t)=\frac{V^{\prime }(\alpha )}{{\hat{c}}^{\prime }(t-\alpha )-\hat{c}}$$

$${\gamma }_{X_2}^{\alpha }=\text{partialderivative}\hat{\gamma }(\alpha ((X_2,t))X_2={\hat{\gamma }}^{\prime }\text{partialderivative}\alpha X_2=-\frac{V^{\prime }(\alpha )}{\hat{c}}\text{partialderivative}\alpha X_2=-\frac{v_{X_2}^{\alpha }}{\hat{c}}$$

$$a^{\alpha }=\text{partialderivative}\hat{v}(\alpha (X_2,t))t={\hat{v}}^{\prime }\text{partialderivative}\alpha t,a^{\alpha }(X_2)=\frac{V^{\prime }(\alpha )}{1-({\hat{c}}^{\prime }/{\hat{c}}^2)X_2},{\tilde{a}}^{\alpha }(t)=\frac{V^{\prime }(\alpha )}{1-({\hat{c}}^{\prime }/\hat{c})(t-\alpha )}$$

$$t^{\text{shock}}=\underset{\alpha \in [0,\infty )}{\min }{t}_{\infty }^{\alpha }$$

Full text

††journal: Journal of the Mechanics and Physics of Solids

1 Introduction

Solids are capable of transmitting mechanical shear waves, which can evolve into shocks depending on the nonlinearity of the shear response. Generating such shear shocks in stiff materials, such as metals, would require extreme strain rates at high stress and would lead to catastrophic damage and failure by other competing processes before shock wave phenomena can be observed (Marchand and Duffy, 1988; Mercier and Molinari, 1998). By contrast, even weak impacts can generate shear shocks within small distances in soft solids, owing to their low shear moduli and significantly higher nonlinearity in shear (Catheline et al., 2003). Many soft solids including biological tissues typically have a shear modulus that is orders of magnitude smaller than the bulk modulus and can hence be assumed to be incompressible. This paper concerns the nonlinear evolution of shear waves into shear shocks in such materials.

Nonlinear shear waves have been studied extensively by both the continuum mechanics and nonlinear acoustics communities. Propagation of finite amplitude plane shear waves and shear shocks were investigated by Chu (1964, 1967), in incompressible isotropic elastic materials. Therein, the condition for shear shock formation was established and the shear loading problem of an incompressible elastic half space was studied using the method of characteristics. Collins (1966, 1967) extended the investigation of wave propagation to materials with transverse isotropy and allowed for transverse displacements in two directions. Davison (1966) studied nonlinear shear waves and propagation of shock waves that are generated by impact loadings in compressible hyperelastic materials wherein accounting for coupled longitudinal motion is necessitated by the compressibility. There, an instantaneous constant loading is considered whereby a shock would immediately form at the loading surface, as opposed to evolution of a smooth loading waveform into shocks as considered by Chu (1964). Recently, Ziv and Shmuel (2019) applied the formulation developed by Davison (1966) to specific compressible hyperelastic material models. Aboudi and Benveniste (1973, 1974) developed a finite-difference based scheme to study such impact induced nonlinear waves. Solutions to the wave equations for compressible hyperelastic materials was provided by Destrade and Saccomandi (2005), including the possibility of dissipation, extending the pioneering works of Carroll (1967, 1974, 1977a, 1977b, 1979).

The nonlinear evolution of shear waves has also been extensively studied in acoustics, where it is routine to perform expansions of the strain energy density function, which is essentially an assumption of weak nonlinearity. The paraxial approximation of small but finite111The shear response would be linear for infinitesimal displacements. wave amplitudes is also employed to obtain reduced wave equations. A nonlinear parabolic wave equation for shear wave beams in isotropic solids was first derived by Zabolotskaya (1986) accounting for nonlinearity, viscous dissipation and diffraction. In the absence of diffraction, the equation reduces to the Modified Burgers equation (MBE) which is similar to the Burgers equation except the quadratic nonlinearity term is replaced by a cubic one. The MBE was studied in detail by Lee-Bapty and Crighton (1987). An alternative expansion of the strain energy density, suitable for application to soft solids, was provided by Hamilton et al. (2004) and used to derive equations that describe nonlinear propagation of plane shear waves for different polarizations (Zabolotskaya et al., 2004). The plane wave model was extended to account for diffraction by Wochner et al. (2008). Recently, Destrade et al. (2019) extended the formulation to any isotropic incompressible solid without having to rely on expansions of the strain energy density.

While the theoretical investigation of shear shocks in solids dates back to the 60’s, the first observation of shear shocks in an elastic medium was only reported recently (Catheline et al., 2003). This was achieved by application of a transient elastography technique to measure the displacements induced by transverse vibrations applied on one end of a gelatin phantom sample. A more recent study by Espíndola et al. (2017) employed a similar method and demonstrated the spontaneous evolution of smooth shear waves into shear shocks in a porcine brain and reported acceleration magnifications of up to a factor of $8.5$. Hence, it was suggested that shear shock waves could be an unappreciated damage mechanism that could play a significant role in traumatic brain injury. To reliably study shear shock generation in the brain and other soft solids, it is vital to consider realistic loadings and the possibility of large deformations. However, both the experimental studies discussed, apply harmonic loading, and make use of reduced wave equations and strain energy density expansions that apply for small amplitude deformations and weak nonlinearity.

In light of these recent studies, in this work we utilize a large deformation elastodynamics framework to study shear shock evolution using realistic constitutive relations that can capture the strongly nonlinear response of soft and biological materials, with the goal of answering the fundamental question: Can shear shocks be induced in soft materials, such as brain tissue, within the length of the impacted object, when subjected to realistic loadings at large deformations? To answer the question, we conduct a parametric analysis using three different constitutive relations that capture a wide range of incompressible material response, and a family of loading waveforms.

The paper is organized as follows: In the following section, we begin by defining our problem setting and deriving the governing equations. Then, in Sec. 3, we solve the boundary value problem by employing the method of characteristics, as in Chu (1964), and derive general expressions needed to evaluate the shock formation distance. Subsequently, in Sec. 4, we apply these expressions to analyze various stress responses, and arrive at closed form expressions for the shock formation distance. These expressions are used to construct non-dimensional phase maps that determine regimes in which a shock can exist and illuminate the effect of material properties and loading scenarios on shock evolution. A parallel supplementary investigation of the distance taken for realization of a given finite acceleration magnification is also carried out. Finally, in Sec. 5, we demonstrate application of our results to the problem of shear impact of the brain to answer the question posed earlier. We provide some concluding remarks in Sec. 6.

2 Nonlinear Shear Wave Equation and Method of Characteristics

Consider a homogeneous and isotropic semi-infinite medium whose undeformed stress-free configuration is described by the Lagrangian coordinates $\mathbf{X}=(X_{1},X_{2},X_{3})$, with $-\infty<X_{1},X_{3}<\infty$, and $X_{2}\geq 0$, as shown in Fig. 1(a). The body is subjected to shearing motion by imposition of a continuous time dependent shearing velocity $V(t)$, on its surface $X_{2}=0$, along $X_{1}$ as shown in Fig. 1(b). Using symmetry considerations and incompressibility, the mapping between the Lagrangian coordinates of a material point, $\mathbf{X}$, and its current coordinates, $\mathbf{x}$, at time $t$, is necessarily given by

[TABLE]

where $u(0,t)=\int_{0}^{t}V(t)\textrm{d}t$ is the shear displacement of the surface $X_{2}=0$. The shear strain $\gamma(X_{2},t)$ is given by

[TABLE]

and the particle velocity is given by

[TABLE]

The deformation gradient $\mathbf{F}$ and the left Cauchy Green deformation tensor $\mathbf{B}$ can then be written as

[TABLE]

and the invariants of $\mathbf{B}$ are given by $I_{1}=I_{2}=3+\gamma^{2}$ and $J=\textrm{det}(\mathbf{F})=1$.

In our analysis we consider a purely mechanical theory and restrict our attention to incompressible hyperelastic materials. The strain energy density per unit volume, $W$, can then be written in its most general form as

[TABLE]

Defining $W_{i}=\partialderivative{W}{I_{i}}$ for $i=1,2$, the first Piola-Kirchhoff stress tensor $\mathbf{P}$ and the Cauchy stress $\mathbf{T}$ are given by

[TABLE]

where $p=p(X_{2},t)$ is a pressure field that arises due to the incompressibility constraint and is determined by the boundary value problem. Substituting the deformation field (4) in (6), we arrive at a stress state of the form222Note that different authors (ex: Chu (1964), Horgan and Murphy (2011)) might report seemingly different longitudinal stress expressions for the simple shear problem, but they all differ only by a hydrostatic term which can be absorbed into the undetermined pressure field $p$.

[TABLE]

where the stress component $P_{12}(=T_{12})$ is denoted by $\tau(\gamma)(=2\gamma(W_{1}+W_{2}))$ and is referred to as the shear stress. Note that $\tau(\gamma)$ is an odd function.

2.1 Equation of Motion

While dilatational stresses are transmitted instantaneously in an incompressible elastic medium, shear stresses are transmitted at a finite velocity. As the surface $X_{2}=0$ is put to motion to generate deformation, a shear wave propagates into the undeformed material and its propagation is governed by balance of linear momentum, which, in absence of body forces, reads

[TABLE]

Here $\dot{\mathbf{v}}$ is the material acceleration, $\rho_{0}$ is the constant density, and integration is performed on a subregion $\mathcal{D}_{0}$ in the reference configuration. In the absence of shocks, assuming the smoothness of fields, this balance law can be localized as

[TABLE]

Substituting (7) into (9) gives us the following equations,

[TABLE]

While perfect incompressibility allows us to impose an arbitrary normal traction $\sigma_{n}=P_{22}$ that will be instantaneously equilibrated throughout the body (as seen from (10b)), in a compressible material, coupling between the longitudinal and shearing deformations might become significant (see for example Ziv and Shmuel (2019)).

Plugging (2) into (10a) gives us the nonlinear wave equation

[TABLE]

where in writing the expression for the wavespeed $c(\gamma)$ the tacit assumption has been made that $\tau^{\prime}(\gamma)>0$, to ensure hyperbolicity of the wave equation such that shear stresses are transmitted through waves. Note that $c(\gamma)$ is an even function. At the linear elastic limit shear waves are transmitted at the constant wavespeed $c_{s}$ given by

[TABLE]

where $\mu$ is the linear elastic shear modulus. In terms of strain and velocity, written in (2) and (3), equation (11) can be written as a system of equations,

[TABLE]

where (13b) is the compatibility condition.

It should be noted that, in practice, the simple shear deformation (4) is a non-trivial one to generate, especially in the large deformation settings, as it requires suitable tractions to be applied on the inclined surfaces, see Destrade et al. (2012), Horgan and Murphy (2011), and Horgan and Murphy (2012). Generating a dynamic simple shear deformation in finite dimensional blocks would be even more tedious, requiring time and spatially varying tractions on the inclined surfaces and would not correspond to any real life shearing scenarios. In the present study we consider shearing of an elastic half-space which can be considered an approximation of a physical shearing setting with large in-plane dimensions $X_{1}-X_{3}$, such that effects from the free surfaces can be neglected in the mid-section where the deformation can be assumed to be one of simple shear.

2.2 Shear shock condition

The shear wavespeed $c$ is a function of the local shear strain $\gamma(X_{2},t)$. Hence, if the local wavespeed increases along the wave propagation direction, a given wave form will spatially spread over time. On the other hand, if the wavespeed decreases along the propagation direction, the waveform would steepen spatially over time and can eventually develop a local discontinuity when a shock forms. Accordingly, for a shear wave travelling in the positive $X_{2}$ direction, shocks can form when $\partialderivative{c}{X_{2}}=\derivative{c}{\gamma}\partialderivative{\gamma}{X_{2}}<0$. Hence, using the expression for $c(\gamma)$ from (11) and the fact that $\tau(\gamma)$ is an odd function, we obtain the condition for shear shock formation as333Here we have used the definitions $\tau^{\prime\prime}(|x|)\equiv\evaluated{\tau^{\prime\prime}(\gamma)}_{\gamma=|x|}$, $\frac{\partial|\gamma|}{\partial X_{2}}\equiv\frac{\partial(\gamma\textrm{sign}(\gamma))}{\partial X_{2}}$, and we consider situations in which shear strain does not change sign. Also, whenever we discuss sign of $\tau^{\prime\prime}(\absolutevalue{\gamma})$ we consider non-zero shear strains.

[TABLE]

Formation of a shock thus depends on the loading program - loading or unloading the material, as defined by the sign of $\frac{\partial|\gamma|}{\partial X_{2}}$ and on the nonlinear shear response of the material through the sign of $\tau^{\prime\prime}(|\gamma|)$. A material with shear stress response such that $\tau^{\prime\prime}(|\gamma|)>0$ will thus evolve a smooth shear loading waveform into a shock when it is being sheared, and one with a shear stress response such that $\tau^{\prime\prime}(|\gamma|)<0$ will produce a shear shock when unloaded from a sheared state. Fig. 2 shows a loading shear waveform steepening into a shear shock versus spreading out, depending on the material constitutive shear response. If $\tau^{\prime\prime}(\gamma)=0$, that is if the material shear response is linear ($\tau=\mu\gamma$) and the wavespeed is constant, the waveform will be relayed at the constant speed $c_{s}$ without steepening or spreading out (dashed lines in Fig. 2).

The most commonly used neo-Hookean (Rivlin, 1948) and Mooney-Rivlin (Mooney, 1940) hyperelastic models both predict a linear shear stress response, in which case there would be no nonlinear shear wave evolution. Biological tissues such as the brain have a strain stiffening shear response (Mihai et al., 2017, 2015; Pogoda et al., 2014; Storm et al., 2005) and unsurprisingly commonly used hyperelastic models for soft solids and biological tissues such as the Gent model (Gent, 1996; Horgan, 2015) and the Fung model (Fung, 1993) predict stiffening shear response i.e $\tau^{\prime\prime}(|\gamma|)>0$ (see Mihai et al. (2015)). Thus soft tissues and other materials whose constitutive response is well captured by these models would allow for shock formation in shear loading. Hence, in this work, we focus our attention on shear loading (as opposed to unloading from a pre-sheared state). Note that the velocity $V$ and the strain $\gamma$ will be of opposite signs by virtue of their definition. Thus when the elastic half-space is sheared along positive $X_{1}$ direction (i.e $V>0$) it develops a negative shear strain.

While (14) gives us the condition under which a shear wave can nonlinearly evolve into a shock, it does not provide information about the length and time scales over which the nonlinear evolution of the smooth waveform into a shock happens, or how this evolution depends on the material nonlinearity and the waveform. To that end, we will analyze the transient evolution of shear waves up to shock formation. We will then apply the analysis to different constitutive stress responses.

2.3 Method of Characteristics

To investigate the transient evolution of shear waves we follow the formulation in Chu (1964) (results recapitulated here) which makes use of the method of characteristics, we first introduce the auxiliary function $Q(\gamma)$ and its derivatives

[TABLE]

The monotonicity of $Q$ ensures that a unique value corresponds to a given shear strain $\gamma$. Further, $Q(\gamma)$ is an odd function ($Q(-\gamma)=-Q(\gamma)$) as $c(\gamma)$ is an even function. Using the above definitions, we can reduce the system of equations (13) to two ordinary differential equations whose solution is

[TABLE]

The variables $\alpha$ and $\beta$ identify sets of characteristic curves described by $\derivative{X_{2}}{t}=\pm c$ along which characteristic relations $\differential{v}=\pm c\differential{\gamma}$ hold. The general solution (16) can be applied to our boundary value problem to solve for the transient shear wave evolution up to formation of shocks. Upon formation of discontinuities the assumption of smoothness used to localize the integral form of momentum balance (8) breaks down, and thus the method of characteristics solution will not hold in regions with shocks.

3 Boundary Value Problem

Summarising our boundary value problem we have the initial and boundary conditions

[TABLE]

where we consider $V(t)$ to be a continuous piece-wise differentiable function of $t$. Using the method of characteristics, the solution for the velocity and strain fields is444 There is no dependence on $\beta$ (from (16)) as the initial shear strain in the body is constant (and zero). (Chu, 1964),

[TABLE]

where fields written as a function of $\alpha$ have $\hat{()}$ accents and, as such, $\hat{Q}(\alpha)=Q(\hat{\gamma}(\alpha))$. According to (18), the velocity $v$ and shear strain $\gamma$ (or equivalently $Q$) are constant along characteristic lines identified by the characteristic variable $\alpha$. The equation of a characteristic line identified by $\alpha$ is (Chu, 1964)

[TABLE]

where $\hat{c}(\alpha)=c(\hat{\gamma}(\alpha))$ and the characteristics are labelled such that $\alpha=t$ on $X_{2}=0$ so that the characteristic $\alpha$ sets out from the loading surface at $t=\alpha$.

Equation (19) determines a relation

[TABLE]

which allows us to convert the solution field from functions of $\alpha$ to functions of $(X_{2},t)$. When characteristics meet, a shock is formed, $\alpha$ becomes multivalued at a given $(X_{2},t)$, and the solution (18) is invalid beyond this time as previously discussed. In this work we are only concerned with the transient nonlinear evolution of shear waves up to the first formation of a shock.

3.1 Shock formation

We now apply the formulas from the previous section to evaluate the length and time scales needed for the development of a shock discontinuity. Representative characteristic lines in the $t-X_{2}$ plane are shown in Fig. 3(a). The slope of a characteristic line setting out at the loading surface $X_{2}=0$ at time $t=\alpha$ is $1/\hat{c}(\alpha)$. Consider two characteristics described by some $\alpha>0$ and $\alpha+\differential{\alpha}>0$, that set out from the loading surface within an infinitesimal time separation $\differential{\alpha}$, with slopes $1/\hat{c}$ and $1/(\hat{c}+\differential{\hat{c}})$ respectively; if $\hat{c}^{\prime}>0$ they can meet and form a shock at a time

[TABLE]

and at the location (using (19))

[TABLE]

where, from here on, quantities with the superscript $\alpha$ are assosciated with the characteristic line $\alpha$. It can be verified that the shock condition $\hat{c}^{\prime}>0$ is consistent with the discussion in Sec. 2.2 (see A).

From (21) and (22), the time and distance taken for characteristics to meet depend inversely on $\hat{c}^{\prime}=\derivative{c}{\gamma}\hat{\gamma}^{\prime}$, where $\derivative{c}{\gamma}$ is smaller in magnitude for materials with weaker nonlinearity of shear response, and $\derivative{\hat{\gamma}}{\alpha}={\derivative{\gamma(0,t)}{t}}$ is the strain rate generated at the loading surface which can be related to the velocity loading rate by differentiating (18b),

[TABLE]

Accordingly, (21) and (22) can be written as

[TABLE]

Functions of $\alpha$ can be rewritten as functions of the strain and vice versa using (18b), i.e $\hat{F}(\alpha)=F(\hat{\gamma}(\alpha))$. Thus higher the loading rate ($|V^{\prime}|$) and/or stronger the material nonlinearity (higher $|\derivative{c}{\gamma}|$), the quicker the characteristics meet and over shorter distances.

The velocity gradient on a given characteristic $\alpha$, denoted by $v^{\alpha}_{X_{2}}$, can be found using (18a) and (19) as function of the time $t(>\alpha)$ or location $X_{2}$,

[TABLE]

Similarly, the strain gradient, $\gamma^{\alpha}_{X_{2}}$, can be found using (23) and (25), as

[TABLE]

Thus, using (21) or (22) in (25) and (26), it can be seen that when characteristics meet the magnitude of local spatial gradients of the velocity and strain become infinite. The local acceleration on a given characteristic $\alpha$, denoted by $a^{\alpha}$, can be also written as a function of $X_{2}$ or $t$ using (18a) and (19),

[TABLE]

and its magnitude is also seen to become infinite when characteristics meet (using (21) or (22) in (27)).

The earliest time at which characteristics meet is when the first shock forms and the method of characteristics solution will not hold at later times. Fig. 3(b) shows a set of representative charactersitics starting out at different times $\alpha_{i}$ meeting at different times $t_{\infty}^{\alpha_{i}}$, the earliest time at which characteristics meet, $t_{\infty}^{\alpha_{1}}$, is the time at which the first shock forms. Thus, in the general case, the time at which the first shock is formed is given by

[TABLE]

Using (24b), the location at which it forms, denoted by $X_{\textrm{shock}}$, can be found as

[TABLE]

The characteristics that start out at $t=\alpha_{\min}^{\infty}$ are the ones that meet earliest in time to form a shock. When $X^{\alpha}_{\infty}$ and $t^{\alpha}_{\infty}$ are smooth functions of $\alpha$, which is true when the material response and loading are smooth (from (24)), we can differentiate (22) and use (21) to write

[TABLE]

Thus the sign of $\derivative{X^{\alpha}_{\infty}}{\alpha}$ and $\derivative{t^{\alpha}_{\infty}}{\alpha}$ will be the same and $\alpha_{\min}^{\infty}$ will also be a minimizer of $X^{\alpha}_{\infty}$, that is,

[TABLE]

Hence, among the locations of intersection of characteristics, the location at which the first shock forms, $X_{\textrm{shock}}$, is also at the shortest distance from the loading surface, when the loading and material behaviour are smooth.

3.2 Acceleration magnification

While the spatial steepening of a waveform ultimately results in a shock, we would also like to quantify the nonlinear evolution of the shear wave as it progresses into the material. One way to do this is to analyse the magnification in acceleration along a characteristic with respect to the initial imposed acceleration. Physical systems are finite in length and we might be interested in knowing if a shock can form within a certain length for a given loading and if not what the maximum magnification in acceleration or field gradients would be, since higher velocity and strain gradients or equivalently higher accelerations, can cause damage in the material. Also, real materials are not perfectly elastic, they have some viscosity which has a spatial smoothening influence on the evolving waveform, hence a pure shock (discontinous fields) is never achieved, so it makes sense to talk in terms of acceleration/field gradient magnification which is a staple of nonlinear wave evolution irrespective of whether viscosity is included in the modelling.

The magnification in acceleration along a characteristic $\alpha$ with respect to the imposed acceleration at the loading surface ($a(X_{2}=0,t=\alpha)$), can be written as a function of progressing $t$ or $X_{2}$ using (27),

[TABLE]

Using (25), (26), and (32), it can be shown that the magnification in the velocity and strain gradients along a characteristic $\alpha$, with respect to their values at the loading surface, is equal to the magnification in acceleration, $M^{\alpha}$,

[TABLE]

Thus, henceforth, whenever we talk about magnification in acceleration, it is tantamount to talking about magnification in the velocity or strain gradients.

From (32), it can be seen that when $\hat{c}^{\prime}>0$ (converging characteristics in the $t-X_{2}$ plane as shown in Fig. 3(a)), the magnitudes of velocity gradient, strain gradient, and acceleration, along a characteristic $\alpha$, magnify as the wave progresses into the material and become infinite when characteristics meet and form a shock (seen by using (21) and (22) in (32)). On the other hand, when $\hat{c}^{\prime}<0$ (diverging characteristics in Fig. 3(a)), the magnitudes of acceleration and field gradient along the characteristic decrease, and the wave locally spreads out. If $\hat{c}^{\prime}=0$ (parallel characteristics in Fig. 3(a)) there is no magnification of fields. For a material with linear shear response $\hat{c}^{\prime}=\derivative{c}{\gamma}\hat{\gamma}^{\prime}=0$ everywhere and the imposed waveform is relayed with no nonlinear evolution at the constant linear elastic shear wavespeed ($v(X_{2},t)=V(t-X_{2}/c_{s})$ using (18a) and (19)). Since we are considering the loading of strain stiffening materials in this paper, the discussion hereon is for555Recall that $\hat{c}^{\prime}=\derivative{c}{\gamma}\hat{\gamma}^{\prime}$. At zero shear strain, $\derivative{c}{\gamma}$ can be zero, and thus $\hat{c}^{\prime}$ can be zero. $\hat{c}^{\prime}\geq 0$.

From (32), the time at which a given magnification in acceleration, $M$, is achieved along a characteristic $\alpha$, denoted by $t^{\alpha}_{M}$ , is given by

[TABLE]

and the location at which it is achieved is given using (19) and (34) as

[TABLE]

Note that for $\hat{c}^{\prime}\geq 0$ and $t<t^{\alpha}_{\infty}$, $M\in[1,\infty)$ using (32) and thus $\lambda_{M}\in[0,1)$. The earliest time at which an acceleration magnification $M$ is realized in the material, $t_{M}$, and the characteristic along which it is realized earliest, $\alpha_{\min}^{M}$ , are given by

[TABLE]

and the location at which the acceleration magnification $M$ is realized earliest, $X_{M}$ , is given by (using (35) and (36))

[TABLE]

Using (23) we can rewrite eqs. 34–37 as

[TABLE]

When the loading and material behaviour are smooth, we can differentiate (35) and use (34) to write

[TABLE]

Note that for a finite magnification $M$, we have $\lambda_{M}<1$ and thus, the signs of $\derivative{X^{\alpha}_{M}}{\alpha}$ and $\derivative{t^{\alpha}_{M}}{\alpha}$ need not be the same. Hence, in general, $\alpha^{M}_{\min}$ need not be the minimizer of $X^{\alpha}_{M}$ and thus, the location at which a given finite acceleration magnification is realized earliest in time need not be the shortest distance from the loading surface at which it is realized. The analysis for shock formation from the previous section is related to the acceleration magnification analysis by taking the limit $M\to\infty$, namely,

[TABLE]

From eqs. 25–27 and eqs. 32–33, it can be seen that a local discontinuity in the imposed waveform ($V^{\prime}\to\infty$) would mean infinite local acceleration and field gradients but not necessarily an infinite magnification. By basing our investigation on analysis of the acceleration magnification (32), we distinguish between two types of discontinuities: those that are evolved due to nonlinearity of the material; and those that are imposed by a discontinuity in the loading waveform $V(t)$. The discussion of shock waves in this paper refers to the former, which is associated with the intersection of characteristics, and in turn represents the limit of acceleration/field gradient magnification becoming infinite. The question may arise as to why one might want to draw the distinction between an imposed and a nonlinearly evolved discontinuity; in a physical system with even minute viscosity a high imposed local gradient would be immediately smoothened and a linear material will not try to steepen the smooth wave into a shock, whereas in a nonlinear material (with $\hat{c}^{\prime}<0$), the material would continuously try to steepen the waveform, increasing the strength of the field gradients/acceleration with time while the viscosity tries to dissipate them away (Bland, 1965).

3.3 Immediate shock formation

If the very first characteristics ($\alpha\to 0^{+}$) meet immediately, then a shock is instantaneously formed at $(t,X_{2})=(0,0)$ and there would be no need to perform the minimisation in (28) and (29). To see when this could happen we first find the time it takes for the very first characteristics to meet by setting $\alpha\to 0^{+}$ in (24a). Without loss of generality, we make the assumption $V\geq 0,\gamma\leq 0$ henceforth. When the initial shear strain is zero, we have $\hat{\psi}(0)=\psi(0)$, yielding

[TABLE]

where we have used $\eqref{charmeettime2}^{2,3}$ and (12). Since $\derivative{c}{\gamma}$ is bounded near zero shear strain, the only way for a shock to immediately form is when the initial loading rate, $V^{\prime}(0)$, is infinite. However, this condition alone is insufficient, since for common materials the slope $\derivative{c}{\gamma}\to 0$ may balance the loading singularity, leading to $t^{0}_{\infty}>0$, as seen by examining the above equation. The distance at which the very first characteristics meet is given, using (22), as

[TABLE]

The time taken for the very first characteristics to realize an acceleration magnification $M$, and the distance at which it is realized, are similarly given by setting $\alpha\to 0^{+}$ in (38a) and (38b)

[TABLE]

It can be seen from (41) and (43) that when the very first characteristics meet immediately ($t^{0}_{\infty}=0$), any acceleration magnification $M\geq 1$ is also immediately realized ($t^{0}_{M}=t_{M}=0$).

3.4 Solution for a general family of loadings

Our goal is to study the effect of material behaviour and the loading waveform on the shock evolution. To that end, we choose a general loading profile

[TABLE]

wherein the shearing velocity is ramped up in a continuous and monotonous fashion from zero at $t=0$ to a ramp velocity $V_{0}$ within a ramp time $t_{0}$, after which the velocity is held constant. Different values of the loading exponent $n>0$ reflect different ramping forms, as shown in Fig. 4, and in particular, $n=1$ corresponds to linear ramping where the acceleration is constant during the ramping phase. Note that as discussed earlier in Sec. 3.2, an infinite initial loading rate ($t=0$, $n<1$) is not an evolved shock.

For the loading (44), using (18a) and (18b), the velocity and strain fields can be seen to be constant for characteristics $\alpha>t_{0}$ and thus the characteristics are parallel and there is no magnification in acceleration along charactersitics $\alpha>t_{0}$ (setting $V^{\prime}(\alpha)=0$ in (38a) gives $t^{\alpha}_{M}\to\infty$ for $t>t_{0}$ and $M>1$). Thus we can restrict the minimisation in (36) and (38) to the ramping phase

[TABLE]

Henceforth, unless mentioned otherwise, we consider $\alpha\in[0,t_{0}]$. The loading is smooth in the interval of interest and hence, as discussed in Sec. 3.1, $X_{\textrm{shock}}$ is also the shortest distance from the loading surface where characteristics intersect (assuming smooth material response). Using (18b), we can specialize (38b) for the loading (44) as (see B.1)

[TABLE]

The time $t_{M}^{0}$ at which the very first characteristics attain acceleration magnification $M$, from (43), can be specialised for the loading (44) as (see B.1)

[TABLE]

If $p_{0}$ is the order of the Taylor expansion of $\derivative{c}{\gamma}$ for $\gamma\to 0^{-}$, i.e if

[TABLE]

where $c_{0}$ and $p_{0}$ are constants666$c_{0}>0$ for a shock to form and $p_{0}\geq 0$ since physically in the linear elastic limit we expect $|\derivative{c}{\gamma}|$ to be bounded., then using (48) in (47) yields

[TABLE]

We can define a critical loading exponent, denoted by $n_{c}(\leq 1)$,

[TABLE]

so that for $n<n_{c}$, we have $t^{0}_{M}=X_{M}=0$ (using (49) and $\eqref{tMzero}^{2}$), i.e there exists a material nonlinearity dependent loading exponent below which a shock forms immediately irrespective of the ramp velocity and ramp time. For $n>n_{c}$, we have $t^{0}_{M},X^{0}_{M}\to\infty$ (using (50) in (49), $\eqref{tMzero}^{2}$), that is the initial characteristics are parallel and there is no acceleration magnification along them (there will still be acceleration magnification and shock formation along later characteristics). For $n=n_{c}$, the values of $t^{0}_{M}$ and $X^{0}_{M}$ are finite and given by

[TABLE]

Essentially from (43), $t^{0}_{M}$ and $X^{0}_{M}$ are inversely proportional to the limit of the product ${V^{\prime}(\alpha)\derivative{c}{\gamma}}$ near zero strain and thus there is a competition between the loading rate and the material’s ability to resist nonlinear evolution decided by $|\derivative{c}{\gamma}|$ (the smaller it is the less nonlinear the shear response and the easier it is for the material to relay the waveform as a spatially smooth one). For $n<n_{c}$, the loading rate singularity immediately overpowers the material’s ability to resist shock formation (resisting characteristics intersecting / acceleration magnification limit from becoming infinite). The weaker the nonlinearity, the higher the $p_{0}$ (from (48)) and thus smaller the $n_{c}$ using (50), meaning more powerful/singular loadings ($n<n_{c}$) are required to immediately form shocks for materials with weaker nonlinearity.

An essential function required for the minimisation in (45) is the derivative $\derivative{t^{\alpha}_{M}}{\alpha}$, which for the loading (44) can be written as below (see B.2 for derivation),

[TABLE]

In the following sections we will apply the formulas developed here to investigate transient evolution of shear waves into shocks for materials with different constitutive responses. The Gent (Gent, 1996) and Ogden (Ogden, 1972) models are widely used for modelling the response of soft solids. We first use an exponential stress model which is not derived from any hyperelastic strain energy, but is one that greatly simplifies the analysis, and different features of the results therein carry over to the Gent and two parameter Ogden hyperelastic models (for parameter values that predict stiffening) analyzed subsequently. Finally we apply the results of our analysis to a physical problem of shear impact of the brain.

4 Analysis for different nonlinear shear constitutive responses

4.1 Exponential stress model

A simple function for the shear response that captures the nonlinear strain stiffening reported in experiments and allows for easy analytical tractability for our problem, is an exponential stress-strain response of the form

[TABLE]

where $k>0$ is a dimensionless parameter that quantifies the nonlinear stiffening. This function reduces to a linear elastic response in the limit of $\gamma\to 0$. However, it is not derived from any established hyperelastic strain energy function and is simply being used as a phenomenological model for analytical tractability; many features of the results carry over to the analysis for Gent and two parameter Ogden hyperelastic models in the subsequent subsections. For the response (53), we have $\tau^{\prime\prime}(|\gamma|)>0$ and hence shocks will be formed in loading. Using (11), the wavespeed is given by

[TABLE]

Fig. 5 shows plots of the stress response and the shear wavespeed for different values of the stiffening parameter $k$.

The derivative of the wavespeed is given by

[TABLE]

Thus from (48), we have $p_{0}=0$ and $c_{0}=c_{s}k/2$ for (55), and using (50) we get $n_{c}=1$. Thus, for $n<1$, we have $t_{\textrm{shock}}=t_{M}=0$ and $X_{\textrm{shock}}=X_{M}=0$. This is not an artifact of the initial loading rate being infinite for $n<1$ but instead represents the limit of acceleration magnification becoming infinite. In the subsequent sections for Gent and Ogden solids it will be shown that an infinite initial loading rate can still require a finite time and distance for the magnification in acceleration to increase or become infinite. Note that this deviation from conventional hyperelastic models essentially stems from the $\textrm{sign}(\gamma)$ factor in the stress response that allows for a quadratic term in the stress expansion whereas conventional energy expansions would only allow a cubic term. For $n\geq 1$ we have to perform the minimization in $\eqref{tm2}^{2}$ to find $\alpha_{\textrm{min}}^{M}$.

Further, we write,

[TABLE]

where we define a characteristic scaling velocity $V_{*}$ (not to be confused with velocity of a characterisitic). As mentioned earlier we assume $V\geq 0$, $\gamma\leq 0$. Using (18b), we have

[TABLE]

Using (57), we have for the exponential solid (refer C.1)

[TABLE]

Combining eqs (52), (58) and (18b) we have

[TABLE]

For $n=n_{c}=1$ we have $\derivative{t^{\alpha}_{M}}{\alpha}>0$ and thus $\alpha_{\textrm{min}}^{M}=0$. Acceleration magnification and shock formation thus initiate at the leading edge of the waveform for $n=n_{c}$ and using (37) and $\eqref{tm0_ncrit}^{2}$ we have

[TABLE]

where we define the loading mach number $m$($\geq 0$) and the characteristic length scale $L_{*}$ as

[TABLE]

Using (59) to perform the minimization in $\eqref{tm2}^{2}$ for $n>1$, we end up with (see C.1)

[TABLE]

where we have defined the the threshold mach number $m_{\textrm{th}}$ as

[TABLE]

The threshold mach number is an increasing function of both the loading exponent $n$ and the target acceleration magnification $M(=1-1/\lambda_{M})$. Additionally, $m_{\textrm{th}}\to 0$ as $n\to 1$. From (62), we see that for loadings with $n>1$, the earliest realization of target acceleration magnification or shocks can happen either at the the trailing edge of the waveform ($\alpha=t_{0}$) or in the ramping part ($0<\alpha<t_{0}$) depending on $m$ and $m_{\textrm{th}}$.

Using (62) in (38b), and the result (60), we have $X_{M}$, the location where an acceleration magnification of $M$ is realized earliest in the solid, non-dimensionalized by $L_{*}$ as (see C.1)

[TABLE]

The solution is continuous at $n\to 1^{+}$ and at $m\to m_{\textrm{th}}^{-}$ but is discontinuous at $n\to 1^{-}$ (see C.1). Setting $M\to\infty$, that is $\lambda_{M}=1$, in (64) gives the non-dimensionalized shock location,

[TABLE]

where $m_{\textrm{th}}=(n-1)/(n+1)$. The expressions in (64) and (65) provide us the location of first realization of given acceleration magnification and of first shock formation respectively.

Non-dimensional maps: The result in eq. (65) allows us to make non-dimensional phase maps, such as the representative one shown in Fig. 6(a) for a given loading exponent $n\geq n_{c}$, that would allow us to determine whether or not a shock would form within a given distance from the loading surface. The maps can serve as a guide in design problems. For example, they can be used to estimate the maximum length one could design for, while avoiding shock formation for a given material and loading777Neglecting reflections in a finite dimensional system.. Similarly they can be used to design loadings that will avoid shock formation given a material and a length scale, or to choose a material that will avoid shocks within a given length scale for a given loading. In Fig. 6(b), we use (65) to make such maps for different values of the loading exponent $n$, where we have the loading mach number $m=kV_{0}/(2c_{s})$ on the vertical axis and the distance for first shock formation non-dimensionalized by $L_{*}=c_{s}t_{0}$ on the horizontal axis. We can see that $X_{\textrm{shock}}/L_{*}$ reduces with increasing loading mach number for any $n$ (can be shown for any $n\geq 1$ from (65)). For $n>1(n_{c})$ there is a threshold mach number separating different functional forms of $X_{\textrm{shock}}$.

Effect of material parameters: We first fix the loading waveform and study the effect of material parameters on shock evolution. The material response is quantified by the nonlinearity parameter $k$, which influences the loading mach number $m$, and by the linear elastic shear modulus $\mu$, which appears through $c_{s}$ in both $m$ and $L_{*}$. Since $X_{\textrm{shock}}$ reduces with increasing $m$ and reducing $L_{*}$, higher material nonlinearity (higher $k$) and lower $\mu$ reduce the distance taken for first shock formation.

Effect of loading: Next, we fix the material properties and analyze the effect of the loading parameters. The loading mach number $m$ increases with $V_{0}$ and hence, a higher ramp velocity reduces $X_{\textrm{shock}}$ (for a fixed $n$ and $t_{0}$). The ramp time $t_{0}$ enters through $L_{*}=c_{s}t_{0}$, and thus $X_{\textrm{shock}}$ scales linearly with $t_{0}$ for given $V_{0}$ and $n$. Thus as $t_{0}\to 0$, we have $X_{\textrm{shock}}\to 0$ , i.e as the loading waveform approaches a discontinuity, a pure shock is immediately formed. The loading exponent $n$ seems to have the most profound and non-trivial effect on shock evolution. For example, it can be seen that going from $n\to 1^{-}$ to $n=1(n_{c})$ may be the difference between a shock immediately forming versus at a distance $L_{*}/m$ (see (65)). From Fig. 6(b), it can be seen that at a given loading mach number the trend in first shock location with respect to the loading exponent is non-monotonous. Also, the maps for different loading exponents intersect each other, meaning that there are different loading mach number regimes in which different loading exponents could delay shock formation. This could have implications for design of protective structures, if we know the mach regime we are operating in for a given problem, we can select a target waveform (by choosing $n$) that could delay the shock formation.

Acceleration magnification: Using (64) we can make non-dimensional maps as shown in Fig. 7, where we have the loading mach number on the vertical axis and the location of first realization of a given acceleration magnification $M$, non-dimensionalized by $L_{*}$ on the horizontal axis. From Figs. 7(a),(b), it can be seen that $X_{M}$ shows a similar qualitative dependence on the material and loading parameters as $X_{\textrm{shock}}$ and the commentary for shock formation carries over. However, when $n>1(n_{c})$ and $m>m_{\textrm{th}}$, the distance at which a given finite acceleration magnification is realized earliest in time need not be the shortest distance at which it is achieved (see C.1). In Figs. 7(c),(d), we have maps of ${X_{M}}/{L_{*}}$ for two different loading exponents and a range of increasing $M$, the limit $M\to\infty$ corresponds to ${X_{\textrm{shock}}}/{L_{*}}$. For a given material and loading it takes longer distances for first realization of higher acceleration magnifications, as one would expect. More importantly, it can be seen that the evolution from high acceleration magnifications ($M\sim 5$) to a shock wave ($M\to\infty$) happens over much shorter distances than the distance taken for realization of those high acceleration magnifications (starting from $M=1$). This is true for any $n\geq 1(n_{c})$ though only two loading exponents have been shown here. If the ramping is linear ($n=1$) so that the acceleration imposed at the loading surface is constant during ramping, maps like Figs. 7(c),(d) can also be used to predict the maximum acceleration realized within a length scale888Neglecting reflections in a finite dimensional system.. This might be useful for application to biological systems where the damage tolerance might be quantified in terms of maximum acceleration.

Post first shock formation: Note that we have restricted our discussion so far to the onset of shock formation and have not touched on the strength of the shock. The first shock that forms would be weak, with the jump in velocity and strain fields being small. Assuming a steady state pure shock ultimately propagates into the material, the weak initial shock would eventually grow into a stronger shock with a jump in velocity from [math] to $V_{0}$ across the shock, and this would happen over longer distances for materials with weaker nonlinearity of shear response. The strength of the steady state shock would be higher for higher ramp velocity $V_{0}$ and a consideration of the shock strength would be pertinent for materials with viscosity as it would be tougher to dissipate away larger discontinuities. Post the first shock formation we would have to apply the integral form of balance laws in a small volume surrounding the shock while using (16) in the smooth regions, however $f(\beta)$ (from (16)) will no longer be constant behind the shock and there will be information travelling along negative $X_{2}$ direction as well (see D for a relevant discussion). The problem becomes analytically intractable post shock initiation and would have to be studied numerically wherein numerical viscosity effects are inevitable.

4.2 Gent model

The popular and widely used Gent hyperelastic model (Gent, 1996) is based on the concept of limiting chain extensibility where the strain energy density function is designed to have a singularity when the first invariant of $\mathbf{B}$, reaches a limiting value $I_{m}=J_{m}+3$,

[TABLE]

For the simple shear deformation (4), this translates to limiting the magnitude of shear strain to a locking strain of $\gamma_{m}=\sqrt{J_{m}}$ as the shear stress keeps increasing (using (66) in (7)),

[TABLE]

Plots of the shear stress-strain curves for various values of $J_{m}$ are shown in Fig. 8(a), the smaller the $J_{m}$ the more nonlinear the shear response. When $J_{m}\to\infty$, the strain energy density function in (66) reduces to that of the neo-Hookean solid, yielding a linear shear stress-strain response. For finite $J_{m}$, we have $\tau^{\prime\prime}(|\gamma|)>0$ and hence, a solid whose shear response is well modelled by the Gent model will produce shear shocks in loading. We define a fractional shear strain $\delta$ which is the ratio of the shear strain and the locking strain $\gamma_{m}$,

[TABLE]

Once again, without loss of generality we assume $V\geq 0$ and thus $\gamma,\delta\leq 0$. Functions of the shear strain $f(\gamma)$ can be written as $\tilde{f}(\delta)=f(\delta\sqrt{J_{m}})$ where the tilde superscript denotes that the argument of the function is $\delta$. Using (11), the wavespeed for the stress response (67) is given by

[TABLE]

Here, as in (56), we have defined a stiffening parameter dependent characteristic scaling velocity $V_{*}$. Plots of the wavespeed as a function of the shear strain are showin in Fig 8(b) for various values of $J_{m}$. All of these curves collapse onto a single curve when plotted as a function of $\delta$ as shown in Fig. 9. We can also write (18b) in terms of the fractional shear strain,

[TABLE]

For the Gent stress response (67), we have

[TABLE]

where $\tilde{Q}$ is an odd function. We can divide (70) and (71) by $V_{*}$ to write

[TABLE]

where $|\tilde{q}(\delta)|$ can be thought of as the local mach number dependent on the local velocity field $V(\alpha(X_{2},t))$ (different from the loading mach number $m$ which depends on ramp velocity $V_{0}$ and is constant for a given loading). Equation (70) provides a mapping $\delta=\hat{\delta}(\alpha)$, allowing us to switch function arguments between $\delta$ or $\alpha$ i.e $\hat{F}(\alpha)=\tilde{F}(\hat{\delta}(\alpha))$. Thus, we can rewrite (46) as

[TABLE]

which can be evaluated further for the Gent solid as (see C.2)

[TABLE]

The loading mach number $m$ and characteristic length scale $L_{*}$ are defined in (61).

The derivative of the wavespeed (69) is given by

[TABLE]

Thus, from (48), $p_{0}=1,c_{0}={3c_{s}}/{J_{m}}$ for (75), and from (50) we have the value of critical loading exponent as $n_{c}=0.5$. Thus, for a Gent solid, even though the initial loading rate is infinite for $0.5<n<1$, the limit of acceleration magnification does not become infinite immediately. It does however for $n<0.5$. For loading exponent values $n\geq 0.5$, we have to perform the minimization (36) to find $\alpha_{\min}^{M}$ needed in (73). Writing the different functions in terms of fractional strain $\delta$ and then plugging them into (52) we end up with (see C.2)

[TABLE]

For $n=0.5(n_{c})$, it can be shown that for $0<|\delta|<1$, the function $\tilde{G}(\delta)$ is positive and hence, $\alpha^{M}_{\textrm{min}}=\delta^{M}_{\textrm{min}}=0$. Thus, acceleration magnification and shock formation initiates at the leading edge of the waveform once again for $n=n_{c}$, and using (37) and $\eqref{tm0_ncrit}^{2}$ we have

[TABLE]

For $n>0.5$, the plot of $\derivative{t^{\alpha}_{M}}{\alpha}$ vs $|\delta|$ is qualitatively similar to the plot in Fig. 10(a), with a zero, $\delta_{0}\in(0,1)$. That is, for $|\delta|<\delta_{0}$, the derivative $\derivative{t^{\alpha}_{M}}{\alpha}$ is negative, for $|\delta|>\delta_{0}$, it is positive and is zero at $|\delta|=\delta_{0}$. Thus if the loading mach number is high enough such that a fractional shear strain magnitude of $\delta_{0}$ is realized at the loading surface during the loading, then ${\delta_{\textrm{min}}^{M}}$ will be999Minus sign since $\delta_{0}>0$, ${\delta_{\textrm{min}}^{M}}<0$.$-\delta_{0}$. We can evaluate the threshold mach number $m_{\textrm{th}}$, which is the minimum value of the loading mach number required to realize a fractional strain magnitude of $\delta_{0}$ during loading, by setting $m_{\textrm{th}}=\tilde{q}(\delta_{0})$ using (72) (note that $\delta_{0}>0$). When $m<m_{\textrm{th}}$, the loading is not high enough to generate fractional strain magnitudes greater than or equal to $\delta_{0}$ and since $\derivative{t^{\alpha}_{M}}{\alpha}<0$ for $|\delta|<\delta_{0}$, the minimizer ${\delta_{\textrm{min}}^{M}}$ would simply be the largest (in magnitude) fractional strain realized. This would correspond to the strain state at the end of ramping when $|q|$ reaches $m$ in (72). Thus,

[TABLE]

Hence, once again, for $n>n_{c}$, depending on the loading mach number and threshold mach number, earliest realization of a shock or given acceleration magnification can happen at the trailing edge ($\hat{\delta}(t_{0})=-\tilde{q}^{-1}(m)$) or the middle of the ramping part of the waveform. Plots of $\delta_{0}$ as a function of the loading exponent $n$ for different values of $\lambda_{M}$ are shown in Fig. 10(b), it can be seen that $\delta_{0}\to 0$ as $n\to 0.5$. Also, $\delta_{0}$ is higher for higher acceleration magnification for a given loading exponent $n>0.5$, and for a given acceleration magnification, $\delta_{0}$ is higher for higher $n$.

Plugging (78) into (74), and using the result (77), we obtain the solution for non dimensionalized $X_{M}$ for Gent solids (see C.2),

[TABLE]

where $m_{\textrm{th}}=\tilde{q}(\delta_{0}(\lambda_{M},n))$ and $\delta_{0}(\lambda_{M},n)$ is the positive root of $\tilde{G}(\delta)$ in (76). Setting $\lambda_{M}=1$ gives the expressions for ${X_{\textrm{shock}}}/{L^{*}}$ for the incompressible Gent hyperelastic solid. Once again it can be verified that the only discontinuity in the solution space is in moving from $n\to n_{c}{}^{-}$ to $n_{c}$ ($0.5$ for Gent solid).

The results for a Gent solid are shown in Figs. 11-12 and can be seen to be exactly similar in nature to the results for the exponential stress model except the fact that $n_{c}$ is now $0.5$ instead of $1$ and thus shocks are immediately formed for loadings with loading exponent $n<0.5$ and the threshold mach number separating different functional forms of the non-dimensional evolution distances comes in the picture for $n>0.5$. Increasing material nonlinearity (lower $J_{m}$ instead of higher $k$) once again hastens shock evolution and other observations made for the results of the exponential model carry over.

4.3 Ogden model

Another popular hyperelastic energy function commonly used to model incompressible soft solids is the Ogden model which is a multiparameter model expressed in terms of the principal stretches $\lambda_{i}$, which are the square root of eigenvalues of $\mathbf{B}$, and material constants $\mu_{p}$ and $N_{p}$ where $p=1,2,3...,P$. For an incompressible material $\lambda_{3}=1/(\lambda_{1}\lambda_{2})$ and we have the strain energy density function (Ogden, 1972),

[TABLE]

where the linear elastic shear modulus is given by

[TABLE]

The Ogden model allows for a large number of parameters that can be chosen to fit different experimental stress-strain curves. Here we use the two parameter Ogden strain energy function for our analysis, i.e $P=1$, such that $N=N_{1}$ and $\mu=\mu_{1}N/2$ are the two independent material parameters that qualify the stress response of the material. The nonlinearity of the shear response is captured by the stiffening parameter $N$ and the shear response is given by (using (80) in (7))

[TABLE]

The sign of $N$ makes no difference to the shear stress response and hence we consider $N>0$ here without loss of generality. Plots of the shear stress as a function of the shear strain are shown in Fig. 13(a), for different values of the material parameter $N$. The value of $N=2$ yields the neo-Hookean shear response. We restrict our attention to materials with $N>2$ for which the response is stiffening. The shear wavespeed for (82) can be found using (11), as

[TABLE]

The large deformation Ogden model is not amenable to find closed form solutions using our analysis as the integration of $c(\gamma)$ is not analytically tractable (to find $Q(\gamma)$). With the goal of demonstrating similarity of results to the other stress models we make the assumption of moderate strains (such that $\gamma^{2}>>\gamma^{4}$). This allows us to use series expansions about $\gamma=0$, allowing for closed form expressions with the caveat of loss of accuracy at larger shear strains. Accordingly, the wavespeed (83) and its derivative can be expanded about $\gamma=0$ as

[TABLE]

Once again we assume $V\geq 0,\gamma\leq 0$ and from (48) we have $p_{0}=1,c_{0}=c_{s}(N^{2}-4)/8$ for (84). Therefore the critical loading exponent $n_{c}$ from (50) is once again $0.5$, similar to the Gent model. Thus for loading exponents below $0.5$, shocks are immediately formed for the stiffening Ogden solid (this result does not require the assumption of small/moderate strains). Further, $Q(\gamma)$ can be approximated as

[TABLE]

We define a modified strain $\delta=\gamma\sqrt{N^{2}-4}/4$, and functions of $\delta$ once again have a tilde accent. Functions of the shear strain $f(\gamma)$ can be written as $\tilde{f}(\delta)=f(4\delta/\sqrt{N^{2}-4})$ and thus the wavespeed from (84) becomes

[TABLE]

Plots of the wavespeed as a function of the modified strain for different material stiffening parameters $N>2$ are shown in Fig. 13(b), at moderate values of $\delta$ all the curves match well with the approximation (86). It can be seen that the more the stiffening (larger $N$), the larger the range of modified strain $\delta$ over which the approximation for wavespeed (86) is good101010At higher $N$ larger range of $\delta$ still corresponds to low/moderate shear strains..

We can use (70) and (73) for the Ogden solid as well, except that $\delta$ is now the modified strain and $\tilde{Q}$ is given by

[TABLE]

Equation (70) once again provides a mapping $\hat{F}(\alpha)=\tilde{F}(\hat{\delta}(\alpha))$. We can divide (70) and (87) by $V_{*}$ to write

[TABLE]

and $|\tilde{q}(\delta)|$ is once again the local mach number, different from the loading mach number $m$. We can specialize (73) for the Ogden solid as (see C.3)

[TABLE]

where $m$ and $L_{*}$ are defined in (61).

Evaluating all the different functions in terms of $\delta$ and substituting them in (52) we end up with (see C.3)

[TABLE]

For $n=0.5$, once again we can verify that $\derivative{t_{M}^{\alpha}}{\alpha}=\tilde{G}$ is positive for all $\delta$ and thus, $\alpha^{M}_{\textrm{min}}=\delta^{M}_{\textrm{min}}=0$. Hence, acceleration magnification and shock formation initiates at the leading edge of the waveform form for $n=n_{c}$, as before. For $n=0.5(n_{c})$, we have using (37) and $\eqref{tm0_ncrit}^{2}$

[TABLE]

For $n>0.5$ the plot of $\derivative{t^{\alpha}_{M}}{\alpha}$ vs $|\delta|$ resembles that in Fig. 10(a), with a zero, $\delta_{0}\in(0,1)$. The discussion of $\delta_{0}$ for the Gent solid holds for the Ogden solid as well except $\delta$ here refers to the modified strain instead of the fractional shear strain and Fig. 14(a) shows the plots of $\delta_{0}$ for an Ogden solid. Similarly (78) also applies for the Gent solid where $\tilde{q}$ is now given by (88). Plugging (78) into (89), and using the result (91) we obtain the solution for non dimensionalised $X_{M}$ at moderate strains for stiffening ($N>2$) Ogden solids (see C.3),

[TABLE]

Here $m_{\textrm{th}}=\tilde{q}(\delta_{0}(\lambda_{M},n))$ and $\delta_{0}(\lambda_{M},n)$ is the positive root of $\tilde{G}(\delta)$ in (90). Setting $\lambda_{M}=1$ gives the expressions for ${X_{\textrm{shock}}}/{L^{*}}$ for the incompressible stiffening Ogden solid. Once again it can be verified that the only discontinuity in the solution space is in moving from $n\to n_{c}{}^{-}$ to $n_{c}$ ($0.5$ for Ogden solid). The results for a stiffening Ogden solid (at moderate strains) are shown in Fig. 14(b) and Fig. 15, and are similar in nature to the results for the Gent model. A larger $N$ corresponds to higher material nonlinearity similar to a larger $k$ in the exponential stress model and smaller $J_{m}$ in the Gent model. All the observations made for the results of the Gent model also hold for the stiffening Ogden model when the strains are moderate ($\gamma^{2}>>\gamma^{4}$). The general nature of results might hold for large strains as well though closed form solutions were not attainable here.

5 Example application to shear shock formation in brain tissue

Having quantified shear shock evolution for different hyperelastic models in the previous section, we now apply our results to a physical problem of interest. Recently, Espíndola et al. (2017) demonstrated the evolution of smooth shear waves into shear shocks in a porcine brain, showing an acceleration amplification factor of up to $M\approx 8.5$ at the shock front; the magnification is not infinite as in our study, since the brain tissue has some viscosity. Nonetheless, it was seen in the previous section that high magnifications of that order are attained practically close to the shock location. In Espíndola et al. (2017), it was suggested that shear shock waves could be a previously unappreciated mechanism that could play a significant role in traumatic brain injuries (TBI). To investigate shear shock formation in the brain using our results, we idealize the brain as an isotropic hyperleastic body and neglect boundary effects and reflections. We also ignore viscosity and any spatial variation in properties. The Ogden model is the most commonly used hyperelastic model for constitutive studies of the brain and consequently we will use our results for the two parameter Ogden model.

The mechanical response of the human brain was carefully and systematically studied across different brain regions for various loading modes in Budday et al. (2017) and it was demonstrated that popular hyperelastic models are unable to accurately capture the stress response across different loading modes with a single set of parameters. However, we are primarily interested in the shear response of the brain and we make use of the parameters reported for the simple shear testing. For the two parameter Ogden shear response (82), averaged parameter values for $\mu$ in the range of 0.3-1.1\kPa and for $N$ in the range of $-26$ to $-21$ were reported for different regions of the brain. We take the density of the brain to be $\rho_{0}=1000\ \textrm{kg/}\textrm{m}^{3}$ which yields $c_{s}$ in the range of $0.5-1.1$ m/s. Assuming that a shear impact of the brain can be represented by the loading (44), we consider a linear ramping ($n=1$) from zero velocity to an impact velocity $V_{0}$, over an impact time of $t_{0}=100$ ms. Using these loading and material parameters in (92), we can make a plot of the shock formation distance as a function of the impact velocity as shown in Fig. 16. The two curves shown are for the material parameter values, from the range of reported values, that predict minimum and maximum shock distance. Thus, for a given impact we have a range of distances within which a shock is expected to form for an arbitrary brain sample. As discussed in the previous section, a brain that has a lower $\mu$ and higher nonlinearity parameter ($N$), evolves a shock over a shorter distance for a given impact. Since the Ogden parameter values were reported for tests conducted in the strain magnitude range of $0-0.2$ and since our results for the Ogden solid are more accurate for moderate strains we cap the curves at impact velocity values that would generate a shear strain of111111Maximum value of $V_{0}=Q(\gamma=0.5)$. $0.5$. Shock formation distance for any impact time can be obtained by remembering that the shock formation distance scales linearly with $t_{0}$.

We can readily use Fig. 16 in a design problem such as that of a helmet to avoid shear shocks in the brain. Since the average length of the human brain is $15$ cm, from Fig. 16, we find a conservative121212Viscosity would allow for higher impact velocities without shock formation. maximum velocity of $\sim 6$ cm/s that can be safely transferred to the brain after dissipation by the helmet and the skull, to avoid shear shocks. Also, the order of shock distances predicted by our analysis for the brain tissue is in agreement with the experimental values in Espíndola et al. (2017) where shocks (peak acceleration magnification) occur at distances of the order of $1$ cm for shear velocities of the order $1$ m/s. These results also suggest that experiments conducted on small brain samples in TBI research might not be able to capture shear shocks that could evolve in a real brain subjected to the same loading. Nevertheless, the non-dimensional maps can be used to rescale the loading waveform so that similar features of the shock evolution can be observed at smaller scales. Finally, it can be clearly seen that larger brains are more susceptible to shear shocks.

6 Conclusions

The nonlinear evolution of a loading shear wave into a shear shock was studied for incompressible strain stiffening hyperelastic solids using the method of characteristics in the framework of large deformation elastodynamics. A general family of loading waveforms was considered and closed form solutions for the distance of first shock formation or realization of a given acceleration magnification were obtained for three different stress models. The results were encapsulated in non-dimensional maps, formulated in terms of a loading mach number and the loading exponent. The maps can guide design selection of loading, and material parameters and dimensions, to avoid shock formation. The singularity of the loading required to immediately form shocks was quantified through a critical loading exponent. For weaker loadings a critical mach number, dependent on the loading waveform shape exists above and below which the functional forms of the evolution distances are different. Larger loadings over shorter times and/or higher material nonlinearity lead to shock evolution over shorter distances. The dependence of the shock evolution on the load profile shape was shown to be non-trivial and different shapes were seen to be better at delaying shock formation depending on the loading mach regime. We hope that our non-dimensional maps can guide design of protective structures and also that this work will serve as a stepping stone for future studies of nonlinear wave evolution in solids in the large deformation elastodynamics framework. Future directions include extension to coupled longitudinal and shear waves in compressible solids, inclusion of viscosity, and studying the evolution post shock initiation.

Appendix A Shock condition equivalence

We assume the shear strain does not change sign. From Sec. 3.1, we have the shock condition $\hat{c}^{\prime}>0$. Since $c(\gamma)$ is an even function, we have $c(\gamma)=c(|\gamma|)$. Thus we have $\hat{c}^{\prime}=\partialderivative{c}{|\gamma|}\derivative{|\hat{\gamma}|}{\alpha}$ where the sign of $\derivative{|\hat{\gamma}|}{\alpha}={\derivative{|\gamma(0,t)|}{t}}$ is decided by the loading program (positive for loading and negative for unloading) and $\partialderivative{c}{|\gamma|}$ has the same sign as $\tau^{\prime\prime}(|\gamma|)$. Thus a strain stiffening material ($\tau^{\prime\prime}(|\gamma|)>0$) would develop shocks under shear loading and a strain softening material ($\tau^{\prime\prime}(|\gamma|)<0$) would develop shear shocks in unloading. This is consistent with the discussion in Sec. 2.2.

Appendix B Solution for loading (44)

B.1 Derivation of $X_{M}$, ${t^{0}_{M}}$

We have from (43),

[TABLE]

For the loading (44), using (18b)

[TABLE]

Thus

[TABLE]

Substituting (95) in (38b), we obtain

[TABLE]

Substituting (95) in (93) and using $Q(\gamma)\to c_{s}\gamma$ as $\gamma\to 0$, we can write (we assume $V\geq 0,\gamma\leq 0$),

[TABLE]

B.2 Derivation of ${\derivative{t^{\alpha}_{M}}{\alpha}}$

We have from (38a)

[TABLE]

Thus

[TABLE]

For $\alpha\in[0,t_{0}]$ and the loading (44), using (18b) we can write

[TABLE]

We also have using (23),

[TABLE]

Thus the derivative $\derivative{t^{\alpha}_{M}}{\alpha}$ for $\alpha\in[0,t_{0}]$ is ,

[TABLE]

Using \psi(\gamma)=-{c^{2}}\big{(}{\derivative{c}{\gamma}}\big{)}^{-1}, we have ${R}(\gamma)=-{\derivative{\psi}{\gamma}}{c}^{-1}=2-c{\derivative[2]{c}{\gamma}}{(\derivative{c}{\gamma})^{-2}}$ and thus

[TABLE]

Appendix C Application to material models

C.1 Exponential model

In this section $\gamma<0$,

[TABLE]

From (57) we have $e^{-\frac{k}{2}\gamma}=1+{V}/{V_{*}}$, which along with (104) yields (58).

Evaluation of $\boldsymbol{\alpha_{\textrm{min}}^{M}}$ for $\boldsymbol{n>1}$: Differentiating (59) we have

[TABLE]

Substituting the loading (44) for $\alpha\leq t_{0}$ in (105) we obtain

[TABLE]

Thus $\derivative[2]{t^{\alpha}_{M}}{\alpha}>0$ for $n>1$ and $\alpha\leq t_{0}$, and thus we can find $\alpha_{\textrm{min}}^{M}$ using (59) and (44) as

[TABLE]

Equation (107) has a solution $\alpha_{\textrm{min}}^{M}\in[0,t_{0}]$ for $n>1$ only if the loading mach number $m$ is greater than or equal to $m_{\textrm{th}}$. For $m<m_{\textrm{th}}$, $\alpha_{\textrm{min}}^{M}=t_{0}$ since $\derivative{t^{\alpha}_{M}}{\alpha}<0$ for $\alpha$ $\in$ $[0,t_{0}]$ (Using $V/V_{*}<m<m_{\textrm{th}}$ in (59)). Using (107), for $m\geq m_{\textrm{th}}$,

[TABLE]

Evaluation of $\boldsymbol{X_{M}/L_{*}}$ for $\boldsymbol{n>1}$: Substituting (58) and the loading (44) for $t<t_{0}$ in (38b) we obtain

[TABLE]

From (62) and (44), for $n>1$ and $m<m_{\textrm{th}}$, we have $\alpha_{\min}^{M}=t_{0},V(\alpha_{\min}^{M})=V_{0}$, and thus (109) becomes

[TABLE]

From (62) and (44), for $n>1$ and $m\geq m_{\textrm{th}}$ we have \alpha_{\min}^{M}=t_{0}\bigg{(}\cfrac{m_{\textrm{th}}}{m}\bigg{)}^{\frac{1}{n}}, $V(\alpha_{\min}^{M})/V_{0}=m_{\textrm{th}}/m$, and thus (109) becomes

[TABLE]

where we have used $m_{\textrm{th}}=\frac{(n-1)\lambda_{M}}{n+\lambda_{M}}$.

**Solution continuity: ** We demonstrate the continuity of the solution (64) at $n\to 1^{+}$ and $m\to m_{\textrm{th}}^{-}$. We have $m_{\textrm{th}}=\frac{(n-1)\lambda_{M}}{n+\lambda_{M}}$ and thus for $n\to 1$, $m_{\textrm{th}}\to 0$. This means that only the $m\geq m_{\textrm{th}}$ case for $n>1$ is relevant for $n\to 1^{+}$ in (64). Setting $n\to 1^{+}$ for that case recovers the $n=1$ case,

[TABLE]

Setting $m=m_{\textrm{th}}$ in the $n\geq 1,m\geq m_{\textrm{th}}$ case in (64) yields

[TABLE]

Setting $m\to m_{\textrm{th}}^{-}$ in the $n\geq 1,m<m_{\textrm{th}}$ case in (64) recovers (113),

[TABLE]

**Shortest distance of realization of acceleration magnification: ** Let us denote the shortest distance from the loading surface at which a given magnification $M$ is realized by $X^{s}_{M}$. For $n<n_{c}$, $X_{M}$ and $X^{s}_{M}$ are both zero. For $n>n_{c},m<m_{\textrm{th}}$, the derivative $\derivative{t^{\alpha}_{M}}{\alpha}$ is always negative in the ramping interval (Using $V/V_{*}<m<m_{\textrm{th}}$ in (59)) and thus using (39), $\derivative{X^{\alpha}_{M}}{\alpha}$ is also negative during ramping. Thus the characteristic $\alpha=t_{0}$ will minimize both $t^{\alpha}_{M}$ and $X^{\alpha}_{M}$ and hence once again $X_{M}$ and $X^{s}_{M}$ are equal. For $n=n_{c}$, the derivative $\derivative{t^{\alpha}_{M}}{\alpha}$ is always positive and greater than one (from (59)) and hence, using (39) $\derivative{X^{\alpha}_{M}}{\alpha}$ is also always positive. Thus the characteristic $\alpha=0$ will minimize both $t^{\alpha}_{M}$ and $X^{\alpha}_{M}$, and thus $X_{M}$ and $X^{s}_{M}$ are equal once more. For the case $n>n_{c},m>m_{\textrm{th}}$, the distances $X_{M}$ and $X^{s}_{M}$ will not be the same. The same observations made here can be verified for the Gent and Ogden models as well where $n_{c}$ will be $0.5$ for them.

C.2 Gent model

Evaluating the different expressions for the Gent solid,

[TABLE]

Substituting the above expressions in (52) and using (72) we obtain ($\tilde{G}(\delta)=\tilde{G}(\hat{\delta}(\alpha))=\hat{G}(\alpha)$)

[TABLE]

Using (115) and (120) in (73) we have

[TABLE]

Using (78) in (126) we obtain

[TABLE]

C.3 Ogden model

Under the assumptions of moderate strains, for the Ogden model we have

[TABLE]

Substituting the above expressions in (103) we obtain

[TABLE]

Using (128) and (132) in (73) we have

[TABLE]

Plugging (78) in (136) we obtain

[TABLE]

Appendix D Post first shock formation

Consider the loading of an initially undeformed strain stiffening solid by the load (44) so that we have shocks. Let us assume that after a sufficiently long time from the first shock formation, a steady state pure shock propagates into the material, i.e there is a propagating velocity and strain discontinuity, with a velocity jump from [math] to $V_{0}$ and a strain jump from [math] to $\gamma^{-}$ across the shock. The momentum balance equation (8) and the integral form of (13b), when applied to the case of a propagating pure shock, give us the jump conditions

[TABLE]

where, $\dot{s}$ is the shockspeed and $\llbracket f\rrbracket=f^{+}-f^{-}$ denotes the jump in a field $f$ across a shock, from its value $f^{+}$ ahead of the shock to $f^{-}$ behind it. Eliminating $\dot{s}$ in (138) we obtain

[TABLE]

Using $\gamma^{+}=\tau(\gamma^{+})=v^{+}=0$, and $v^{-}=V_{0}$ in (139), and defining the auxiliary function $\phi(\gamma)=\tau(\gamma)\gamma$ we have

[TABLE]

Consider cases in which $t_{\textrm{shock}}>t_{0}$, that is the first shock does not form before the ramping is completed. Then the largest strain (in magnitude) realized before first shock formation, $\gamma_{1}$, would be the strain state at the trailing edge of the waveform where the velocity is equal to $V_{0}$, and is obtained from (using $\eqref{Qsol}$, $\eqref{Q_def}^{1}$, and definition of $c$ in $\eqref{waveeqn}$)

[TABLE]

Squaring (141) and defining the non-dimensional auxiliary function $\chi(\gamma)=(\int\displaylimits_{0}^{\gamma}\sqrt{\tau^{\prime}(\gamma)}\differential{\gamma})^{2}$ we obtain

[TABLE]

For a general strain stiffening solid the functions $\chi(\gamma)$ and $\phi(\gamma)$ are different, and hence, from (140) and (142) the strain state at the trailing edge of the evolving waveform, $\gamma_{1}$, is different from the strain state behind a steady state shock. As mentioned in Sec. 4.1, when applying the method of characteristics solution (16) behind the shock, $f(\beta)$ is not constant and the $\beta$ characteristics come into play, meaning that there is also information traveling back to the loading surface from the growing shock. This might facilitate the change in strain state from $\gamma_{1}$ to $\gamma^{-}$ for possible steady state shock propagation but this is subject for future research. See Ch. 4 in Hamilton et al. (1998) and the references therein for a discussion on reflection of waves by shocks in gas dynamics. There, it is demonstrated that for weak shock propagation the “reflected wave” is extremely small and can be neglected.