Nonlinear elasticity in rocks : A comprehensive three-dimensional description

We study theoretically and experimentally the mechanisms of nonlinear and nonequilibrium dynamics in geomaterials through dynamic acoustoelasticity testing. In the proposed theoretical formulation, the classical theory of nonlinear elasticity is extended to include the effects of conditioning. This formulation is adapted to the context of dynamic acoustoelasticity testing in which a low-frequency “pump” wave induces a strain field in the sample and modulates the propagation of a high-frequency “probe” wave. Experiments are conducted to validate the formulation in a long thin bar of Berea sandstone. Several configurations of the pump and probe are examined: the pump successively consists of the first longitudinal and first torsional mode of vibration of the sample while the probe is successively based on (pressure) P and (shear) S waves. The theoretical predictions reproduce many features of the elastic response observed experimentally, in particular, the coupling between nonlinear and nonequilibrium dynamics and the three-dimensional effects resulting from the tensorial nature of elasticity.


I. INTRODUCTION
Nonlinear mesoscopic elastic materials (NMEMs) refer to a class of materials that can be mechanically described as an assembly of mesoscopic-sized "hard" elements (e.g., grains with characteristic lengths ranging from tens to hundreds of microns) embedded in a "soft" bond system (e.g., cement between grains, pore space, fluid) [1].The microscopic-sized imperfections at the interfaces between the "hard" and "soft" subsystems are believed to be responsible for a number of exotic properties, beyond the classical theory of nonlinear elasticity [2,3], manifesting themselves in a number of situations.The identification, quantification, and theoretical representation of these properties have been active topics of research over the last 30 years.The fundamental mechanisms behind these effects are not completely understood and delineated yet though, probably due to the multiscale nature of the problem, e.g., microscopic and mesoscopic features of the material structure lead to noticeable effects at the macroscopic scale.NMEMs are not simply a scientific curiosity.Understanding their complex mechanical behavior is a growing need for many applications, including the monitoring of building integrity in civil engineering [4,5], nuclear safety [6,7], bone fragility (e.g., microfractures) [8], the mechanisms of earthquake triggering [9], or the design of novel materials [10].
Hysteresis with discrete memory is one of these exotic properties and has first been evidenced carefully in quasistatic experiments on rocks using nonperiodic stress protocols [11,12].Hysteresis was subsequently modeled using the Preisach-Mayergoyz formalism borrowed from the field of electromagnetics [13,14].Practically, when a sinusoidal elastic wave propagates in a medium exhibiting such a property, it progressively evolves into a triangular wave as demonstrated theoretically [15] and experimentally [16].In the frequency domain, this distortion leads to the generation of odd harmonics only [13].
Slow dynamics is another exotic property typically observed in many dynamic experiments.In brief, under a continuous harmonic excitation, the elastic moduli decrease (i.e., material softening) and eventually reach a new but unstable equilibrium state tuned to the dynamic strain amplitude induced by the excitation, a process known as "conditioning."Conversely, if elastic energy is no longer injected into the system, the material recovers slowly its original elastic properties, with a rate proportional to the logarithm of time [17,18], a process known as "relaxation." The effects of hysteretic nonlinearity and slow dynamics have first been evidenced with resonance experiments [19,20].Typically, a long thin bar with free boundary conditions, which is representative of a one-dimensional (1D) unconstrained system, is vibrated around one of its resonance frequencies.In NMEMs, the resonance frequency starts decreasing as a function of strain once the strain reaches a certain amplitude.Depending on the strain amplitude, there exist different regimes in which the variation of the resonance frequency with strain is dominated by classical nonlinearity, slow dynamics, or both [21,22].These observations have been described theoretically through a number of models derived from 1D theory of elasticity [13,15,[23][24][25].However, as soon as the geometry of the sample or testing conditions deviate from the 1D assumption, these models are no longer applicable [6,26].
Classical nonlinear elasticity is intrinsically three dimensional and should be treated using the appropriate mathematical tools, i.e., tensors.In that sense, NMEMs are different from fluids.One indisputable evidence is the ability to perform noncollinear three-wave mixing in nonlinear elastic solids but not in fluids.Noncollinear wave mixing refers to the experiments where two waves with two different frequencies, f 1 and f 2 , are launched at two different angles mix to generate waves with frequencies [27,28].This process is governed by the third-order elastic constants (TOECs) of the solid.TOECs have been first introduced by Murnaghan in 1937, as part of a general tensorial description of the nonlinear stress-strain relationship in elastic isotropic materials [2].Classical nonlinear elasticity in an isotropic solid is described by three TOECs (l, m, n in [2] or A, B, C in [3]) in addition to the two Lamé constants (λ, μ) used in the theory of linear elasticity.Measurement of these constants in solids has been made possible via the simplification of Murnaghan's theory by Hughes and Kelly [29].They consider the particular case of a high-frequency (HF) pulse of relatively small amplitude (probe) used to probe the local change of elastic wave speed induced by a large deformation (pump).This configuration has been widely adopted for the measurement of TOECs [30,31] or to investigate the effects of stress on wave speed [32][33][34].Murnaghan's theory does not predict hysteresis or slow dynamics in the stress-strain relation.These effects can be observed in a dynamic experiment of acoustoelasticity where the large deformation is no longer static but dynamic [35,36].However, in these experiments, the analysis of the dynamic acoustoelastic (DAE) data is often supported by the 1D theory of elasticity enhanced with nonlinear terms, thus neglecting the three-dimensional nature of the strain field.This may lead to wrong interpretations of the anisotropy induced by the softening effects.Additionally, hysteresis and slow dynamics are treated as two independent mechanisms whereas it has been shown that they are strongly related [37].A number of models have been proposed where such connection between hysteresis and slow dynamics is accounted for [24,25,38,39].However, these models are also all derived from 1D theory of elasticity.
More recently, Lott et al. [40,41] have proposed an extension of Hughes and Kelly's formulation [29] by taking into account the softening of geomaterials under dynamic vibration.This new formulation starts from a volume of material that is small enough to fit within the assumption of continuum mechanics but large enough to contain a relatively large number of defects.In this volume, the defects have random orientations, thus resulting in an isotropic conditioning effect quantified by the scalar parameter α.This parameter is then projected onto the principal strain axes, resulting in the three-dimensional (3D) manifestations of conditioning, as already observed experimentally [26].
In this paper, dynamic acoustoelasticity is studied in three dimensions via an extensive set of experiments where the large amplitude strain induced by the pump is no longer restricted to a longitudinal motion as in previous work but extended to a torsional motion.Likewise, the HF pulses used to probe the elastic state of the material is no longer restricted to (compressional) P waves but extended to (shear) S waves.The experimental findings are compared to theoretical predictions from the approach proposed by Lott et al. [40].

II. THEORETICAL BACKGROUND
The mechanisms responsible for slow dynamics can be described at multiple scales.At a microscopic scale, it can be explained by surface interactions (e.g., clapping and friction) that exist within the interface of a defect (e.g., grain boundaries Wave speed (i axis: direction, j axis: polarization) in most rocks).At larger scales, these interactions can be homogenized and described by thermodynamic laws [18]: a volume of material with a sufficiently large density of defects may be represented as a smooth spectrum of energy barriers that need to be hopped to reach a new elastic state.The existence of this spectrum is supported by the fact that the rate at which the relaxation process occurs is proportional to the logarithm of time.This thermodynamic description can be tied to a friction-based mechanism, as suggested in [18].The global stiffness is proportional to the surface area of microscopic contacts within the media through a static friction coefficient.Under dynamic stress, the total area of microscopic contacts diminishes and so the material appears "weaker."For clarity, the symbols used in this paper are listed and described in Table I.
In this paper, we assume that the microcontacts are randomly oriented and distributed, with a length scale much smaller than the typical acoustic wavelength used in this study.A random orientation and distribution of nonlinear sources implies an isotropic effect on the elastic response.Under this assumption, we extend the formulation of the stiffness tensor proposed by Hughes and Kelly [29] to introduce conditioning.The approach we follow is borrowed from theory of damage evolution in fracture mechanics, described in detail by Zubelewicz [42].The tensorial product between the strain and stress vector bases δ ij is the natural basis for the elastic tensor and should now include conditioning effects as where n i is the principal strain direction, * ij the strain amplitude, the star symbol * denotes the basis formed by the principal strain axes, and α is a scalar quantifying conditioning.The stiffness tensor is then expressed as where λ and μ are the Lamé constants and l, m, and n are the Murnaghan constants [2].Practically, conditioning is applied in the basis formed by the principal strain axes.This basis is obtained from an eigendecomposition of the strain field measured in the geometric basis.This decomposition, = P * P −1 , is always possible because the strain tensor is always real and symmetric.Once conditioning has been applied to the stiffness tensor in the basis formed by the principal strain axes, according to Eq. ( 2), the following transformation is used to express the conditioned stiffness tensor in the geometric basis: C ij kl = P ir P js P kt P lu C * rstu .This result is equivalent to classical acoustoelasticity theory when the term quantifying conditioning is equal to the identity matrix, = I 3 .

III. EXPERIMENTAL ARRANGEMENT
The experimental arrangement is depicted in Fig. 1.Experiments were conducted on a cylindrical sample of Berea sandstone (Cleveland Quarries, Amherst, OH) with a diameter of 25.8 mm, a length of 305.5 mm, a mass density of 2054 kg/m 3 , and a nominal permeability ranging between 500 and 1000 mD.The linear elastic properties of this sample were characterized in previous work using resonant ultrasound spectroscopy [43].It was found that the sample is well described by a homogeneous and isotropic material with a Young's modulus E = 9.9 GPa and a Poisson's ratio ν = 0.068.The sample was instrumented with a piezoelectric disk (PZT-5H ceramic with a diameter of 25.54 mm and a thickness of 6.35 mm) epoxied on one of its flat ends and five shear piezoelectric plates (PZT-5A ceramic with dimensions of 20 × 20 × 3 mm 3 ) epoxied on the round surface of the sample, near the opposite end.The piezoelectric disk induced a longitudinal motion whereas the shear plates, in this particular configuration and when driven in phase, induced a torsional motion in the sample.Longitudinal and torsional motions were recorded on the surface of the sample with an in-plane laser vibrometer (Polytec OFV-552).The two laser beams were shined near the end hosting the piezoelectric disk.Depending on the orientation, the laser beams can sense the particle velocity in the e 1 (axial) or e 3 (tangential) direction [26].Note that, at the particular sensing position shown in Fig. 1, v 3 in Cartesian coordinates is equivalent to v φ in cylindrical coordinates, which is the coordinate system chosen to treat the torsional motion.As part of our experiments of dynamic acoustoelasticity, it is necessary to know the position of maximum strain in the sample since it will be where the HF transducers will be placed to probe the instantaneous elastic state of the sample (i.e., where nonlinear and slow-dynamic effects will be the strongest).Because of the mass added by the piezoelectric disk and plates, the position of the maximum strain is not exactly at the middle of the sample.The position of maximum strain 11 for the first longitudinal mode of resonance and φx for the first torsional mode of resonance were measured accurately with a 3D laser vibrometer (Polytec PSV-3D-500) following the setup described in the Supplemental Material in [26].The relations between the particle velocities measured with the in-plane laser near the end of the sample and the maximum strains of interest for the longitudinal and torsional modes were also established as part of this measurement step.Resonance frequencies were found near 3 kHz for the first longitudinal mode and near 1.8 kHz for the first torsional mode.
Experiments of dynamic acoustoelasticity (DAE) were conducted according to the procedure described by Renaud et al. [35] and Rivière et al. [36].Relative changes in the propagation speed of the HF probe waves, which will be referred to as relative velocity modulation (RVM) in the rest of this paper, are monitored under dynamic loading.One transducer is used to send short pulses centered at 1 MHz that are generated by an arbitrary waveform generator (National Instrument PXI-5412).Once the pulses have been propagated across the diameter of the sample, they are received by the second transducer and digitized at 50 MHz (National Instrument PXIe-5122).The probe system is activated while the sample is excited at the frequency of the first resonance mode of longitudinal or torsional vibration (LF pump).The RVM is obtained by cross correlating the modulated transmitted pulses with the one transmitted before the pump activation.The pulses are sent with a repetition frequency corresponding to 1.01 times the period of the LF pump signal.Provided that the response of the sample has reached a steady state, the RVM induced by one cycle of the pump can be fully captured over 100 cycles of the pump, with a regular spacing equal to 0.01 times the period of the LF harmonic signal.Using this stroboscopic effect, the original time frame can be compressed by a factor of 100 to define a new time frame in which the dynamic response of the material will be described.A representative example of the measured DAE data is shown in Fig. 2(b).In the figure, V 21 /V 0 21 denotes the RVM of the P wave in the "e 2 " direction.The time history of the strain 11 induced by the LF pump at the point of maximum strain in the sample is depicted in Fig. 2(a).When the LF pump is activated, the speed of sound drops rapidly and oscillates in a near steady-state regime (i.e., conditioning process).This regime never actually reaches a steady state, as shown in Fig. 2(b), but for all practical purposes a steady state can be assumed for the 100 cycles during which the stroboscopic probe is analyzed.When the LF pump is stopped, the equilibrium value of V 21 /V 0 21 is recovered smoothly (i.e., relaxation process).

IV. P WAVES AS HF PROBE WAVES
Two series of DAE experiments were conducted using P waves generated via a pair of HF compressional-wave transducers (Olympus V303-SU) to probe the elastic state of the sample locally.In the first series of experiments, the dynamic LF pump consisted of the first longitudinal mode of vibration of the sample, which is the typical DAE configuration found in the literature.In the second series of experiments, the dynamic LF pump is changed to the first torsional mode of vibration.Experiments were conducted for eight amplitudes of the pump.For each configuration (i.e., pump type and amplitude level), experiments were conducted three times to ensure repeatability.

A. Longitudinal mode as LF pump
For a longitudinal mode, the strain field is homogeneous over a cross section of the sample in the plane e 2 − e 3 , which is also the propagation plane of the HF pulses.The strain field can be written in Cartesian coordinates as where A L is a constant equal to the peak amplitude of the strain 11 at the axial position where the strain is maximum for the first longitudinal mode of vibration (where the HF probe is positioned in the experiments).This tensor is already diagonal.Therefore, the transformation matrix P is the identity matrix I 3 and the expression of the stiffness tensor is the same in the geometric basis and in the basis formed by the principal strain axes.The propagation speed V 22 of the HF P wave in the strain field induced by a longitudinal mode is then written as The coupling between the LF pump strain and propagation speed of the HF probe waves through the third-order elastic constants (e.g., Murnaghan constants) is given by Eq. ( 4).The conditioning effects are contained in the term 22 , which is written as 22 = 1 − ανA L .Typical results obtained in this set of experiments are shown in Fig. 3.The RVM of the P wave is not symmetric with respect to the axis 11 = 0, as shown in Fig. 3(a).This asymmetry between the tension and compression phases is due to the thirdorder elastic constants.This data set also exhibits the effects of hysteresis (the same path is not followed for increasing and decreasing strain) and material softening (the average steadystate value of V /V 0 is not zero).These results are expected and consistent with previous work [35,36].

B. Torsional mode as LF pump
For a torsional mode, the strain field is not uniform over a cross section of the sample in the plane e 2 − e 3 .The strain field can be written in polar coordinates as where A T is a constant equal to the peak amplitude of φx /R at the axial position where the strain is maximum for the first torsional mode of vibration (where the HF probe is positioned).The strain field corresponding to the torsional mode is plotted in polar and Cartesian coordinates in Fig. 4. Across the diameter of the sample (r varying from −R to R), half of the propagation is affected by a positive strain amplitude whereas the other half is affected by a negative strain amplitude.The strain field can be expressed on its principal axes as Finally, the propagation speed V 22 of the HF P wave in the strain field induced by a torsional mode is written as In this theoretical description, 11 is equal to one and there is no dependence of the wave speed on the third-order elastic constants.In other words, the propagation speed of the P wave should not be affected by the strain field of a torsional mode.More generally, a probe wave is insensitive to the strain field if both its polarization and propagation direction are aligned with the principal strain axes with zero strain amplitudes.Here, the principal strain axes with nonzero amplitudes are linear combinations of the two basis vectors e 1 and e 3 , as suggested by the transform matrix P T .On the other hand, the polarization and direction of propagation of the HF probe wave are on the "e 2 " axis.This phenomenon has already been observed in experiments by Gallot et al. [44].
Results from the DAE experiments are shown in Fig. 5.There are obvious differences between the results obtained in this configuration and those obtained with a longitudinal mode.First, the the RVM of the P wave in the "e 2 " direction is now symmetric with respect to the zero-strain axis, which agrees with the theoretical prediction on the absence of coupling between the propagation speed of the P wave in the "e 2 " direction and the third-order elastic constants.This result is important as it highlights the necessity of treating dynamic acoustoelasticity with the 3D theory of elasticity, as opposed to a superposition of 1D descriptions.In this configuration, we can also observe that the magnitude of the hysteresis has been considerably reduced, despite the fact that the strain amplitude of the torsional mode is more than twice as large as the strain amplitude of the longitudinal mode.However, we still observe a large curvature of the RVM.This curvature most likely originates from the second-order term of classical nonlinearity δ in the 1D description of elasticity discussed earlier, which is not included in the 3D model of elasticity formulated in this paper.More importantly, we still observe elastic softening, which was not expected from the theoretical prediction (i.e., 11 = 1).A first possible reason for the discrepancies between theory and experiments could be the departure from the plane-wave assumption in the experiments.The theoretical description is based on the propagation of FIG. 5. DAE data for the configuration where the LF pump is the first torsional mode of vibration of the sample and the HF probe is a P wave.(a) RVM of the P wave in the "e 2 " direction as a function of strain for the largest drive amplitude of the pump.(b) Time history of the RVM of the P wave in the "e 2 " direction.(c) Time history of the strain induced by the LF pump at the axial position of the HF probe.Time histories are displayed in the compressed stroboscopic time frame: to obtain the three cycles depicted in (b) and (c), 300 cycles of real data were acquired.The blue and red colors indicate increasing and decreasing strain, respectively.plane waves, which ignores the scattering induced by the HF transducers and cylindrical shape of the sample.This scattering results in volume effects on the wave-speed measurements.A two-dimensional simulation of a P wave traveling across the sample illustrates these effects in Fig. 6.These results show that the direction of propagation is not purely collinear with the "e 2 " axis and is consequently affected by the nonuniform strain field away from this axis.Another possible reason could be the assumption of a scalar parameter α in the theoretical description whereas this parameter is probably tensorial, as suggested by Remillieux et al. [26] in a series of resonance experiments.The anisotropy of microcracking in Berea sandstones was also formally evidenced elsewhere [45].

C. Harmonic content
Data analysis can also be conducted in the frequency domain to extract additional information about the nonlinearity of the material (e.g., quantifying the parameters of classical nonlinearity and hysteresis).The use of a compressed stroboscopic time frame provides a constant time step (in the new compressed time frame) and thus allows such analysis by simply taking the fast Fourier transforms (FFTs) of the RVM waveforms.Figures 7(a) and 7(b) show the spectral content of the RVM of the P wave in the "e 2 " direction at the highest excitation levels of the longitudinal and torsional modes, respectively.The variations of the amplitudes of the fundamental frequency and harmonics of the RVM as a function of the excitation level is shown in Fig. 8 for the two mode types as well.To understand and analyze these results, it is beneficial to make a connection with the 1D theory of nonlinear elasticity used typically to describe DAE experiments.In this description, the relevant modulus of elasticity M is strain dependent and written as M = M 0 (1 + β + δ 2 + H.O.T.), where M 0 is the linear elastic modulus, β is the coefficient of first-order nonlinearity, which can be expressed in terms of in terms of second-and third-order elastic constants, and δ is the coefficient of second-order nonlinearity, which can be expressed in terms of second-, third-, and fourth-order elastic constants in the direction of propagation.In the context of an experiment where an elastic pulse propagates in a nonlinear FIG. 8. Amplitudes of the fundamental frequency (f 0 ), second (2f 0 ), fourth (4f 0 ), sixth (6f 0 ) harmonics as a function of the strain amplitude.The LF pump consists of the (a) longitudinal mode of vibration and (b) torsional mode of vibration.Amplitudes are shown in log scale, relative to their respective amplitude at the lowest strain value.
elastic material, the parameter of classical nonlinearity β can be quantified by monitoring the amplitude growth of the second harmonic of this pulse [46].Since the RVM is a perturbed quantity, the parameter β affects its fundamental frequency (and not its second harmonic).More generally, if a parameter of nonlinearity affects the amplitude of the nth harmonic of an elastic pulse, it will affect the amplitude of the (n − 1)th harmonic of the RVM.
Figure 7 shows that by switching the LF pump from a longitudinal mode to a torsional mode of vibration, the amplitude of the fundamental frequency f 0 vanishes.This means that the effect of β is inhibited.This result is in perfect agreement with the results shown in Figs.3(a) and 5(a) where the RVM of the P wave as a function of strain evolves from being asymmetric with respect to the zero-strain axis (effect of β) and having a curvature (effect of δ) to only having a curvature.In the case of a torsional mode, only even harmonics (2f 0 , 4f 0 , etc) can be observed on the spectrum [see Fig. 7(b)].It can be demonstrated theoretically that these harmonics are mainly generated either from the term of classical nonlinearity δ or from hysteresis.A brief look at Fig. 7(a) can give us a first hint about which of the two is dominating.In this figure, the amplitudes of the fundamental frequency and second harmonic are not sensibly different.This result cannot be explained by classical nonlinearity.Indeed, it was demonstrated that when the clapping mechanism of a contact is modeled with a strong stiffness asymmetry (i.e., different stiffnesses are involved when opening or closing the crack), the harmonic amplitude rapidly decreases with the harmonic order [47].More recently, Zhao et al. [48] proposed a more refined model to study simultaneously the effects of clapping and friction in solids with microcracking.They showed that the parameter of nonlinearity β is mostly caused by the clapping mechanism and stiffness asymmetry whereas it is rather insensitive to the frictional mechanism.Based on these facts, the generation of even harmonics observed in our DAE experiment is most likely dominated by the hysteresis and softening mechanisms.Figures 8(a) and 8(b) show the evolution of the first six harmonics of the RVM as a function of strain amplitude for the longitudinal and torsional modes, respectively.Quadratic and cubic dependencies associated to the terms of classical nonlinearity β and δ, respectively, are also represented for comparison purposes.For a longitudinal pump, the amplitude of the fundamental frequency f 0 grows with the square of the strain amplitude (second power law) and follows almost perfectly the theoretical prediction obtained using the term of classical nonlinearity β.The amplitude of the second harmonic, however, does not follow the cubic dependence expected from the term of classical nonlinearity δ.In other words, the term δ is not fully responsible for the curvature of the RVM observed in Fig. 3(a), as hinted previously.When the LF pump consists of a torsional mode, the amplitude of the fundamental frequency is almost independent of the strain amplitude.This is expected since there is no coupling with third-order elastic constants in this setup.The amplitudes of the second, third, and fourth harmonics grow twice as fast as when a longitudinal mode is used as the LF pump.The asymmetry of the strain field induced by a torsional mode [see Fig. 4(b)] is responsible for the faster growth.In a torsional mode, the principal strain axes change their direction across the diameter of the sample (in Cartesian coordinates, strain experiences a sign change).It may be seen as two independent sources of nonlinearity acting on the spectrum of the RVM, thus artificially increasing the harmonic generation.Those frequencies (2f 0 , 4f 0 , and 6f 0 ) may be identified as the signatures of nonclassical nonlinear elasticity.

V. S WAVES AS HF PROBE WAVES
In this section, the elastic state of the sample is probed via a pair of HF shear-wave transducers (Olympus V154-RM).As in the previous section, the HF probe is placed at the axial position of maximum strain for the LF pump of interest, namely, a longitudinal or a torsional mode of vibration.Strong coupling between the HF transducers and the sample is required for an efficient shear-wave transmission.This is achieved by using clamps, whose effect on the response of the sample is naturally mitigated by the fact that the HF transducers are mounted at a node of vibration of the sample (i.e., zero displacement).The elastic state of the sample is probed at different polarizations of the HF probes, for the highest level at which the LF pump can be operated.The convention used to describe the polarization of the probe is shown in Fig. 9.The angle θ parametrizes the polarization of the probe.In the case of a longitudinal mode, the probe is used from θ = 0 • to 90 • in steps of 10 • .In the case of a torsional mode, the probe is used from θ = 0 • to 180 • in steps of 10 • .In both cases, results at other angles can be inferred by symmetry.For each polarization angle, the stiffness tensor experienced by the HF S wave is expressed in the local basis of the HF probe, i.e., (e 1 ,e 2 ,e 3 ) in Fig. 9. Transformation from (e 1 ,e 2 ,e 3 ) to (e 1 ,e 2 ,e 3 ) is simply achieved by a rotation around the "e 2 " axis, which can be formulated by the rotation matrix.

A. Longitudinal motion
The propagation speed V 21 , of the HF S wave in the strain field induced by a longitudinal mode can be written in the local basis of the HF probe as where is expressed in the basis formed by the principal strain axes and has the following components: 11 = 1 − α , 22 = 1 − αν , and 33 = 1 − αν .Typical results obtained in this set of experiments are shown in Fig. 10.As previously, the RVM waveform shown in Fig. 10(a) can be plotted as a function of the LF pump strain when near steady-state conditions have been reached (e.g., around time t 1 ) to visualize the nonlinear effects experienced by HF This step is repeated for all polarizations of the HF probe, with results for three polarizations shown in Figs.10(b)-10(d).As in the case of the HF P -wave probe, strong hysteresis and asymmetry with respect to the zero-strain axis is observed.When the polarization angle increases from θ = 0 • to 90 • , the magnitude of the conditioning decreases and that of the curvature increases.These data sets can be further reduced to visualize the magnitude of the conditioning as a function of polarization angle, as shown by the polar plot in Fig. 10(e).Conditioning is extracted for all polarization angles at times t 1 and t 2 .Smallest (resp.largest) conditioning effects are observed at θ = 90 • (resp.θ = 0 • ) when the polarization direction of the probe is aligned with the principal strain axis with smallest strain amplitude νA L (resp.largest strain amplitude A L ).This effect from the Poisson's ratio was already observed by Lott et al. [40].During the slow relaxation (e.g., at time t 2 ), the induced anisotropy remains.Eventually, as the material recovers back to its original elastic properties, this nonlinear strain field signature collapses to one point.The theoretical formulation can be used to reproduce, at least partially, the experimental findings.The predicted conditioning as a function of polarization angle is shown in Fig. 10(f).For the predictions, the following parameters were used: λ = 0.73 GPa and μ = 4.64 GPa measured on this sample using resonant ultrasound spectroscopy [43], l = −1600 GPa, m = −3400 GPa, and n = −450 GPa measured by Winkler and Liu [31] on Berea sandstone using static acoustoelasticity, and α = 1600.The theoretical model captures well the features observed in experiments, except near the angle θ = 0 • .One possible reason for the discrepancy could be the change in the contact area between the shear-wave transducer and the sample as the polarization angle varies from 0 to 90 • .At θ = 0 • , the shear wave is launched from a line contact whereas at θ = 90 • it is launched from a point contact.Such issue could be easily resolved by using a sample with a rectangular cross section.However, in this case, the strain field for a torsional mode would be more complex than the one depicted in Fig. 4(b).

B. Torsional motion
The propagation speed V 21 of the HF S wave in the strain field induced by a torsional mode can be written in the local basis of the HF probe as  In this expression, 11 is still equal to one but 22 = 33 = 1 − αA T r.In order to model the asymmetry of the strain field, Eq. ( 9) is evaluated at the radial positions with the largest strain, i.e., at r = −R and R, where φx ≈ ±4 microstrain for this Data from this experiment can be analyzed in the same fashion as previously, for the case of the longitudinal mode.Raw and processed data are shown in Fig. 11.Similarly to the results obtained with a HF P -wave probe (see Sec. IV B), the variation of the RVM as a function of strain [see Figs.11(b)-11(d)] does not exhibit any slope or strong asymmetry with respect to the zero-strain axis, meaning that the effects of the third-order elastic constants on the propagation of the HF S wave are inhibited by the strain field of a torsional mode.In the present case, strong hysteresis and conditioning are observed at all polarizations of the probe.In the case of P waves, no coupling between the P -wave propagation and material softening was expected from the theoretical model.Such coupling is now enforced in the case of S waves by the terms 22 and 33 .Figure 11(e) clearly indicates that the probe wave is most affected by conditioning at θ = 45 • and −45 • , which is aligned with two of the principal strain axes.This is well captured in the theoretical predictions shown in Fig. 11(f).For these predictions, the same second-order and third-order elastic constants were used but the parameter α had to be reduced to 700 to match the magnitude of conditioning observed in experiments.The ratio of α predicted for a longitudinal mode (α = 1600) and that predicted for a torsional mode (α = 700) is approximately the same as the ratio observed experimentally using nonlinear resonant ultrasound spectroscopy by Remillieux et al. [26].Furthermore, it should be noted that the theoretical model does not include viscoelasticity.Therefore, in this model, the material is conditioned and relaxed instantaneously by the applied strain amplitude.For this reason, the shape observed experimentally in Fig. 11(e) can only be reproduced by superposing predictions of the instantaneous conditioning as a function of polarization angle at the radial positions r = −R and R. Including viscoelasticity in the theoretical model as well as other relevant physics is currently undertaken by the authors.

VI. CONCLUSION
A tensorial model of conditioning was proposed and partially validated in this paper.The original formulation of Hughes and Kelly [29] was extended to introduce conditioning effects, through a scalar parameter α that is projected onto the principal strain axes.This model was derived in a very general form and therefore can be applied to arbitrary wave polarizations and strain fields.It captures many features observed in experiments of dynamic acoustoelasticity including (1) inhibition of the effect of third-order elastic constants when a high-frequency P wave propagates in the strain field of a torsional mode (for the particular experimental configuration described in this paper), (2) polarization angles for which the propagation direction of a high-frequency S wave is most affected by the strain fields of longitudinal and torsional modes.By comparing experimental findings with theoretical predictions, it was also suggested that the model could be greatly improved through the use of a tensorial α (the parameter α would then be different for longitudinal and torsional modes) and the use of viscoelasticity theory.

FIG. 2 .
FIG.2.Time histories of the LF pump and HF probe DAE data in the compressed stroboscopic time frame.(a) Time history of the strain 11 when the sample is driven at the resonance frequency of its first longitudinal mode of vibration (LF pump).The LF pump is activated at time t = 0.5 μs.(b) Time history of the RVM of the P wave in the "e 2 " direction (HF probe).

FIG. 3 .
FIG.3.DAE data for the configuration where the LF pump is the first longitudinal mode of vibration of the sample and the HF probe is a P wave.(a) RVM of the P wave in the "e 2 " direction as a function of strain for the largest drive amplitude of the pump.(b) Time history of the RVM of the P wave in the "e 2 " direction.(c) Time history of the strain induced by the LF pump at the axial position of the HF probe.Time histories are displayed in the compressed stroboscopic time frame: to obtain the three cycles depicted in (b) and (c), 300 cycles of real data were acquired.The blue and red colors indicate increasing and decreasing strain, respectively.

FIG. 4 .
FIG. 4. Normalized strain field for a torsional motion in polar (left) and Cartesian (right) coordinates.

FIG. 6 .
FIG.6.Two-dimensional simulation of a P wave propagating from an immersed transducer through the cylindrical sample.Two instants are depicted: soon after the wave penetrates the sample (t = 5 μs) and soon before the wave reaches the receiving transducer (t = 15 μs).The slowness surface of the P wave in the strain field of a torsional mode and a schematic representation of the probed region are also depicted to understand the experimental observations.

FIG. 9 .
FIG. 9. Schematic representation of the polarization angle of the HF S-wave probe.

FIG. 10 .
FIG. 10.DAE data for the configuration where the LF pump is the first longitudinal mode of vibration of the sample and the HF probe is an S wave.(a) Time history of the RVM of the S wave in the "e 2 " direction.(b)-(d) RVM of the S wave in the "e 2 " direction as a function of strain for the largest drive amplitude of the pump at polarization angles of the probe of θ = 0 • , 40 • , and 90 • .(e) Polar plot (function of the polarization angle of the probe) of the measured average value of the RVM (conditioning) around times t 1 and t 2 .(f) Polar plot of the predicted instantaneous value of the RVM (conditioning) at the largest strain amplitude A L .

FIG. 11 .
FIG. 11.DAE data for the configuration where the LF pump is the first longitudinal mode of vibration of the sample and the HF probe is an S wave.(a) Time history of the RVM of the S wave in the "e 2 " direction.(b)-(d) RVM of the S wave in the "e 2 " direction as a function of strain for the largest drive amplitude of the pump at polarization angles of the probe of θ = 0 • , 40 • , and 90 • .(e) Polar plot (function of the polarization angle of the probe) of the measured average value of the RVM (conditioning) around times t 1 and t 2 .(f) Polar plot of the predicted instantaneous value of the RVM (conditioning) at the largest strain amplitudes A T , at the radial positions r = −R and R.

TABLE I .
List of the symbols.