Resonant ultrasound spectroscopy for crystalline samples containing initial strain

Resonant ultrasound spectroscopy (RUS) describes an acoustic or ultrasonic technique for measuring elastic properties of anisotropic crystalline solids by analyzing a sample’s free vibrational frequencies of resonance.^1,2 For measuring the components of the anisotropic elastic stiffness tensor $cijkl$⁠, the experimental frequency response is combined with a predictive inverse model to determine unknown $cijkl$⁠. Prior to the creation of the inverse model, researchers developed the forward model based on variational methods applied to the elastodynamics of a freely vibrating anisotropic solid.^1–5 Visscher et al. enabled the possibility of the inverse model by introducing basis functions amenable to straightforward and computationally efficient integrations.¹ 

Inherent to the traditional RUS technique for determining $cijkl$ is the assumption that the crystalline solid obeys Hooke’s law, which states that the stress and strain are linearly related through the stiffness components $cijkl$⁠. Hooke’s law is fundamental to the vast majority of solid mechanic modeling. However, real crystalline solids exhibit anharmonic potential energy functions, and Hooke’s law is valid only for small deformations in which points in the material remain close to equilibrium. In cases where residual stress or a prestrain is present, the small dynamic deformation couples to the prestrain nonlinearly and requires the use of a nonlinear constitutive relationship. The first-order anharmonicity can be accounted for by introducing a quadratic stress/strain relationship, which is commonplace in several physical models (see Clayton’s book⁶ on nonlinear mechanics, for example).

The integration of anharmonic effects into RUS modeling is not common. A few researchers have investigated anharmonic effects on RUS by considering pressure loading on spheres to observe pressure-dependent elastic constants.^7–9 For this situation, the pressure exhibited on the surface of the sphere must be included as a boundary condition, which complicates modeling as pressure loading is not a traditional natural boundary condition as opposed to rigid or traction free conditions.⁴ More recently, Ji et al. measured third-order elastic constants by measuring the change in cross-sectional resonance frequencies of an aluminum sample under tensile loading and inversion of a semi-analytical finite element model (SAFEM).¹⁰ 

For the case when a sample contains spatially inhomogeneous residual stresses, the wave displacement is nonlinearly coupled with the prestrain. The constitutive behavior for coupling two or more strain fields is aptly described using incremental elasticity theories.^11,12 In the previous work, the authors developed a micromechanical model to help study the role of residual stress on RUS.¹³ A significant influence of residual stress on RUS is also supported by recent experiments on additively manufactured metallic samples.¹⁴ Unlike the pressure loading case, a sample containing residual stress or prestrain is in equilibrium. Thus, the usual traction free boundary condition may be applied, and the RUS model can closely follow the traditional formulation.^1,2,4,5

Our previous work¹³ was focused on the specific problem involving polycrystalline samples containing residual stress and included micromechanical modeling. Often, RUS is applied to non-polycrystalline or homogeneous samples (at least at the length scale of the wavelengths seen in RUS measurements). The previous work may be overly complicated for several RUS applications. Thus, our objective is to provide a simpler model for describing RUS of homogeneous, anisotropic samples with a spatially inhomogeneous prestrain. The choice of including the prestrain rather than initial (residual) stress simplifies the model equations and eliminates the need to include compliance constants to convert to stress. The simplified model equations enable more efficient computations,¹³ which is expected to lead to efficient inverse models in the future. Additionally, focusing on strain offers a natural connection to diffraction based residual stress measurements, such as x-ray and neutron diffraction. Such diffraction measurements are explicitly dependent on strain, whereas residual stress estimates arrive from converting the strain, usually from a linearized strain–stress relationship. Future work will seek to apply the model to simultaneous RUS and diffraction measurements as reported by Torres et al.¹⁵ 

The theory is presented in Sec. II that culminates in the RUS eigenvalue problem dependent on initial strain or prestrain. The eigenvalue problem is solved numerically through the use of the standard Rayleigh–Ritz procedure using polynomial basis functions.¹ The numerical evaluation results from needing to evaluate spatially inhomogeneous volume integrals, which results from the spatially dependent prestrain. In Sec. III, prestrain fields are generated using a simulation of a synthetic polycrystalline sample of a parallelepiped shape undergoing plastic deformation on two of its surfaces, representing the type of residual stresses/strains left in the material after a shot-peening process. Using the resulting simulated prestrain fields, the RUS model is evaluated in Sec. IV.

Using the standard Rayleigh–Ritz procedure provides a strong connection to the previous foundational work that led to RUS becoming the standard method for inversion to measurements of elastic constants.^1,2,5 An inverse problem toward nondestructive characterization of the prestrain field is an exciting future extension to this work.

One of the key constitutive relations in the field of mechanics is the relation between stress and strain. Traditionally, RUS is formulated through a Hookean or linear relationship between stress $σ$ and strain $ϵ$⁠, i.e., $σij=cijklϵkl$⁠, where the components of $cijkl$ are the second-order elastic constants or stiffnesses of the sample. For the present problem, the sample is allowed to contain a spatially inhomogeneous strain field or prestrain, which will be accounted for by the Green or Lagrangian strain tensor $E˚$⁠. While $E˚$ is assumed to be present in the sample prior to the RUS measurement, the constitutive relationship must account for the coupling of $E˚$ and the dynamic strain developed in the sample during the measurement. This situation is handled by the theory of incremental elasticity¹¹ and has been applied previously to similar problems involving elastic wave propagation in samples containing initial (residual) stress.^12,16,17 An appropriate constitutive relationship, which has been derived previously,^12,13,17 relates the first Piola–Kirchhoff stress $P$ to the coupled prestrain and the dynamic displacement gradient $∇u$⁠,

where $cijklmn$ is a tensor with components of the third-order elastic constants or stiffnesses.

Cauchy’s law or the balance of linear momentum can be formed with the first Piola–Kirchhoff stress from Eq. (1),

where $u$ is the displacement, $ρ$ is the density of the material (in the strained configuration), $T˚$ is the initial Cauchy stress, and the dynamic displacement is assumed to have an $eiωt$ time dependence. The tensor $B$ can be viewed as an effective strain-dependent elastic stiffness,

However, $Bijkl$⁠, unlike $cijkl$⁠, does not enjoy the minor symmetries $Bijkl≠Bjikl≠Bijlk$⁠. Taking the dot product of Eq. (2) with $u$⁠, integrating over the volume, applying the product rule followed by the divergence theorem leads to the integral relationship

where the term involving $T˚$ vanishes as the volume averaged stress over the entirety of the sample’s volume must satisfy equilibrium.¹⁸ Furthermore, to arrive at Eq. (4), the traction free boundary conditions,

are applied to the surface integrals formed when applying the divergence theorem. The condition $P⋅n=0$ is the usual boundary condition assumed for the free vibration RUS problem, whereas the condition $T˚⋅N=0$ arrives from the material being in stress equilibrium.¹⁸ Note that the unit vectors $n$ and $N$ are directions normal to surface elements during and prior to the dynamic excitation and, in general, are not equal.

Rayleigh’s quotient, which is an explicit expression for the frequencies $ω$⁠, is commonplace in vibration problems and follows from Eq. (4),

which highlights that the resonances are governed by the volume integrals involving the effective stiffness components $B$ and mass density $ρ$⁠. The volume integrations in Eq. (6) permit both $B$ and $ρ$ to be spatially inhomogeneous, where $B$ is necessarily inhomogeneous because of the physical processes leading to $E˚$⁠.¹⁸ Frequency solutions satisfying Eq. (6) can be obtained through traditional eigenvalue methods.^1,3–5 For evaluation of the model, we employ a quasi-analytical approach by making use of trial functions in the form $ui=∑μCiμΨμ$⁠, where the index $μ$ is summed over the number of trial functions used. For present purposes, we consider the sample to be a rectangular parallelepiped. For this case, the basis functions can be polynomials of the form¹ 

where the sets $n1$⁠, $n2$⁠, and $n3$ are integers mapped to respective index values of $μ$⁠. All permutations of $n1$⁠, $n2$⁠, and $n3$ are included up to a chosen cutoff value of $μ$⁠. Typically, the expansion is truncated at a value $μ$ such that it satisfies the condition $n1+n2+n3≤nmax$⁠, where $nmax$ is a sufficiently large integer appropriately selected to allow for converged solutions. In the previous work,^1,2 a value of $nmax=10$ is appropriate for convergence.

Substituting the trial functions into Eq. (6) gives

where the volume integrals are the tensors

Equation (9) is a restatement of the Lagrangian.¹ The Lagrangian is stationary, and thus, Eq. (9) may be differentiated with respect to the coefficients $Cij$ and set to zero¹ to yield a generalized matrix eigenvalue problem,

where pairs of indices $im$ and $kn$ are mapped to single indices $μ$ and $ν$⁠, respectively. In our evaluation, the volume integration required to form $Γ$ is completed using numerical integration using Simpson’s method. This results from $E˚$ and thus $B$ being spatially inhomogeneous, which excludes the possibility of pre-calculating the volume integrals, which is done when the integrands are homogeneous.¹ Once $M$ and $Γ$ are obtained, the eigenvalues $ω2$ and eigenvectors $C$ can be obtained using routine eigenvalue solvers and analyzed for different prestrains.

The model described in Sec. II explicitly depends on the inhomogeneous prestrain $E˚$ through the nonlinear constitutive relationship. It is important to note that the model is independent of the process leading to $E˚$⁠. Thus, $E˚$ is not restricted to arrive from strictly elastic deformations. This allows RUS to be potentially sensitive to a variety of processing mechanisms. As an example, the present analysis will consider $E˚$ being a consequence of a shot-peening process, used to improve the mechanical properties of a metallic part by inducing a compressive state in its surface. While the simulated shot-peening process is used here for providing input of prestrain into the RUS model, this integration demonstrates the possibility of extending RUS as a characterization tool to a variety of other complex material processes that lead to prestrain or residual stress formation.

The simulations are performed on a synthetic polycrystalline Cu sample of rectangular parallelepiped geometry having dimensions of $3×2×1mm3$ as seen in Fig. 1.

The prestrain fields to be used in the subsequent RUS analysis are simulated using an Elasto-ViscoPlastic Fast Fourier Transform-based (EVPFFT) formulation, whose full description can be found in Ref. 19.

Using Euler implicit time discretization and Hooke’s law, the local elasto-viscoplastic constitutive equation at $t+Δt$ adopted in our EVPFFT calculations is given by

where $σ$ is the Cauchy stress tensor, $cijkl$ is the elastic stiffness tensor, $ε$ is the total strain tensor, and $εp,t$ is the plastic strain-rate tensor at time t. The polycrystalline unit cell is partitioned into three regions, two shot-peening affected regions, corresponding to the material near both $3×2mm2$ surfaces, up to a depth of 0.03 mm, and the third region comprising the rest of the sample (bulk region).

In the shot-penning affected zones, a plastic strain-rate tensor is imposed of the form

which represents the average cumulative effect of individual shot-peening impacts, which produce plastic (isochoric) deformation of the surface. The coefficient $α$ scales the plastic strain rate, which corresponds to compression of the material along the surface’s normal direction, while expanding it along the in-plane directions, simulating the shot-peening conditions similar to those described previously.²⁰ This plastic strain rate is imposed for ten time-increments of $Δt=0.1μs$ on both surfaces, mimicking shot-peening with 100% coverage on the surface. By varying the coefficient $α$⁠, the total amount of plastic deformation that is deposited on the surfaces can be controlled. A higher amount of imposed plastic strain simulates a longer or stronger (higher impact velocity) shot-peening process that locally increases the plastic deformation at the surfaces of the material. In this way, we can adjust the magnitude of the residual state fields to test the sensitivity of RUS to different prestrain levels.

In the bulk region, the plastic strain-rate tensor $ε˙p$ is given by

where $τok$ and $mk$ are, respectively, the slip resistance and the Schmid tensor associated with the slip system $k$⁠, $γ˙o$ is a normalization factor, and $n$ is a stress exponent.

The copper single crystals deform plastically by dislocation glide on 12 ${111}⟨110⟩$ slip systems with an initial slip resistance ranging from 50 to 150 MPa and a linear hardening rate of $h=30$ MPa, where $h$ is the proportionality factor relating the increase of the slip resistance and the accumulated plastic strain.

In the EVPFFT prestrain simulations, the boundary condition imposed to the unit cell corresponds to zero macroscopic strain, conforming with the general experimental observation that the shot-peening process does not produce significant macroscopic deformation of the material. The complementary boundary condition, i.e., zero macroscopic stress, was also investigated, resulting in very minor quantitative differences in the simulated prestrain fields. A total of four simulations were conducted resulting in various levels of prestrain. The case leading to the most significant prestrain is shown in Figs. 2 and 3, showing the residual stress and prestrain profiles. The prestrain components $E˚11≈E˚22$ and $E˚33$ from the midplanes of the sample are shown. The strain component $E˚11$ demonstrates large extensional strains within the interior of the sample volume, whereas compressional strains are present close to the shot-peened faces. The $E˚33$ component displays the opposite behavior. Each of the four cases simulated resulted in an approximately linear increase in the resulting prestrain fields. For example, the average $E˚11$ component over the $z=0$ plane was found to be 476.4, 652.3, 829.9, and 1377.0 (in units of $μ$m/m), respectively.

Section III described the simulation procedure to generate a realistic distribution of prestrain in a synthetic polycrystalline sample of Cu. Now, the RUS model provided in Sec. III is evaluated for the four cases of prestrains generated during the shot-peening simulations relative to the strain-free case. The analysis is performed through first inputting the bulk second- and third-order elastic constants (⁠$cijkl$ and $cijklkmn$⁠) to form the effective stiffness tensor $B$⁠. Here, the input (isotropic) elastic constants were $c12=123$ GPa, $c44=92$ GPa, $c123=−295.14$ GPa, $c144=−248.34$ GPa, and $c456=−197.94$ GPa, which were experimentally measured values reported previously.²¹ It is important to note these elastic constants are for isotropic polycrystalline Cu and different from the single-crystal elastic constants of the individual microstructural grains used in Sec. III.

The four cases of prestrain obtained from the shot-peening simulation led to a nearly linear dependence to the change in the resonance frequency. Thus, in what follows, the prestrain case will refer only to the most extreme case, which provided the largest average strain $E˚11$ at the z=0 midplane. For both the strain-free and prestrain cases, the tensor $B$ was discretized using the regular grid of 192 x 128 x 64 voxels. $B$ for the prestrain case is spatially inhomogeneous as it depends on $cijkl$⁠, $cijklmn$ and the prestrain $E˚$⁠. The spatial inhomogeneity forces the integrals in $Γ$ to be evaluated numerically, which is computationally expensive relative to the strain-free case where $B$ is homogeneous. For the strain-free case, the integrations in $Γ$ have analytical solutions as seen in Visscher et al.,¹ whereas the integrals in $Γ$ were calculated using Simpson’s 1/3 rule²² for the prestrain cases. The analytical solutions provide a useful check to the accuracy of Simpson’s rule calculations. The first two columns of Table I compare the resonant frequencies of the first 12 modes when calculated using the analytical and numerical evaluation.

Deviation of the numerical integration from the analytical solution is typically on the order of $0.0001$ kHz, which is sufficient to resolve changes due to the prestrains. The mode shapes of the first 12 resonant modes in the strain-free configuration are seen in Fig. 4.

The mode shapes are included to investigate whether particular mode types are more or less sensitive to the prestrain. The corresponding resonant frequencies of the first 12 modes for the most significant prestrain case are seen in column 3 of Table I, while the corresponding percent change in the resonant frequencies from the strain-free case is seen in the fifth column. In addition, a graphical representation showing the shift in the resonant frequency due to the prestrain is shown in Fig. 5.

The percent change due to the prestrain is on the same order as the acoustoelastic effect, which describes the change in the wave speed relative to stress or strain, for example, a longitudinal elastic wave propagating in a Cu sample with the propagation direction parallel to uniaxial tensile loading. For a homogeneous uniaxial tensile stress of 100 MPa, the percent change in the phase velocity is around 0.3%.²³ Unlike the case of a tensile stress, the residual stress field must be in equilibrium, which causes locations of positive and negative stress components. Thus, the coupling between the dynamic displacement and the prestrain in the RUS problem competes with regions having different directions of stress components. This comparison leads to the belief that the percent change in the resonant frequency is consistent with acoustoelasticity or the wavespeed dependence on stress, which also utilizes a constitutive relationship like Eq. (1). Furthermore, acoustoelasticity depends on the mode type (longitudinal or shear) and their propagation and displacement directions relative to the direction of the stress. Similarly, the modes in Fig. 4 have different sensitivities to the prestrain. In particular, modes 3, 5, 10, and 11 have the lowest sensitivity to the prestrain, and all have a small displacement component normal to the shot-peened or $3×2mm2$ face. Interestingly, the categories of these modes are mode 3—bending, mode 5—shear, and mode 10—dilation, which indicates that the lack of sensitivity is not correlated to the particular mode type. Furthermore, modes 5, 6, and 11 are the only modes predicted to have their resonant frequency shift downward. The other modes are more sensitive to the prestrain, and all have appreciable displacement components in all three directions. In particular, modes 1, 4, and 12 are torsional modes in which mode 1 is the most sensitive to the prestrain and mode 4 and 12 are two of the most sensitive ones. This finding is similar to acoustoelasticity in that bulk wave modes with displacements in the direction of uniaxial stress have the highest sensitivity to the stress. The small changes in the resonant frequencies mitigate the possibility of mode switching, which complicates modal identification and is beyond the current scope of this paper. Furthermore, the current analysis has only considered Cu. Materials, such as aluminum, display greater nonlinearity and display greater sensitivity to residual stress as seen in the previous work.¹³ It is also expected that samples with larger aspect ratios will have enhanced sensitivity.

Up to this point, the results have been analyzed relative to the forward RUS model of predicting resonances and mode shapes from inputs of elastic constants. However, RUS is most often used in conjunction with the inverse model to predict the elastic constants from inputs of experimentally measured resonances. The present model offers a technique to estimate the error in determining the values of the elastic constants when residual stress or prestrain is neglected in a measurement. An inverse model was constructed using least squares minimization to obtain the elastic constants $c12$ and $c44$⁠. Using the resonant frequencies from the prestrain case, the elastic constants found using the inverse model were $c12=119.6$ GPa and $c44=92.46$ GPa. These values are significantly different from the strain-free case of $c12=123$ GPa and $c44=92$  GPa. Thus, various levels of prestrain can have an appreciable influence on the measured elastic constants. While these differences may seem small, the situation described is highly idealized. In actual experiments, the prestrain could lead to misidentifying modes, which will compound the error. Furthermore, modern RUS setups are extremely sensitive to changes in material properties. For example, Ghosh et al.²⁴ obtained elastic constants to within $≈0.01$ GPa. Thus, neglecting prestrain or residual stress has the potential to lead to spurious results if left unaccounted for.

In this work, resonant ultrasound spectroscopy (RUS) is formulated to explicitly include prestrains stemming from residual stress. By including prestrains, the application of RUS as a characterization tool is broadened as many modern material processing procedures and manufacturing processes cause prestrain and residual stress formation in materials. The influence of prestrain generated from simulating a shot-peening procedure was shown to significantly impact the resonance frequencies of a $3×2×1mm3$ sample. The lowest order bending mode was the most sensitive as its resonance frequency shifted from 386.22 to 388.94 kHz as a result of prestrain. While the shift might seem small, modern experimental RUS systems are extremely precise and will be influenced by the presence of prestrain. For the cases considered here, neglecting the prestrain would result in underpredicting the elastic constant $c12$ by 3.4 GPa.

The present model is related to a previous derivation developed strictly for polycrystalline materials as it included the influence of the polycrystalline microstructure.¹³ The current formulation can be applied more generally to samples that are sufficiently homogeneous, which includes most traditional applications of RUS. Additionally, this work has integrated a plasticity simulation of a commonly applied shot-peening procedure with the RUS framework to investigate how material processing influences resonances. Future extensions of this research will include other modern processing procedures, which will broaden the scope of applicability of RUS.

This research was supported by the Air Force Office of Scientific Research (AFOSR) under the Summer Faculty Fellowship Program (SFFP), and their financial support is graciously acknowledged. Z.F. and R.A.L. acknowledge support from Los Alamos National Laboratory's LDRD Project No. 20220488MFR. C.M.K. would like to thank Julian Maynard for providing references and valuable discussion. C.M.K. and M.C. also thank Anubhav Roy for help with computations in addition to Rasheed Adebisi, Ian Barrett, Tyler Lesthaeghe, John Aldrin, and Michael Uchic for their invaluable discussions and recommendations during the fellowship period.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.