Decoupling Uniaxial Tensile Prestress and Waveguide Effects From Estimates of the Complex Shear Modulus in a Cylindrical Structure Using Transverse-Polarized Dynamic Elastography

Dynamic elastography methods—based on optical, ultrasonic, and magnetic resonance imaging modalities—aim to quantitatively map the shear viscoelastic properties of biological tissue, which are often altered by disease and injury. Most studies to date have focused on larger organs, such as the liver or brain, where boundary effects were assumed negligible. But, as elastography is applied to other anatomical regions where dimensions in at least one direction are smaller or of comparable length to bulk shear wavelengths—such as in slender skeletal muscles, blood vessel walls, the heart wall, and the cornea—boundary effects become non-negligible and must be considered. Researchers using optical elastography to assess the viscoelastic properties of the cornea have long recognized this, adapting models to include waveguides by treating the cornea as a plate-like structure. Here, transverse wave motion on the cornea is modeled as Rayleigh–Lamb waves [1]. Blood vessels, as well, have been modeled using cylindrical shell equations considering fluid-structure interaction [2–4]. Limited studies on cardiac elastography have also acknowledged the frequency-dependent (i.e., wavelength-dependent) waveguide behavior of the heart wall [5].

Often, when elastography studies are done under varying nonzero quasi-static prestress conditions, observed changes in mechanical wave behavior are attributed solely to the nonlinear property of the tissue: it's been suggested that its shear and viscous constants are highly dependent on the tensile load and associated deformation. A recent article provides a summary of the literature relevant to this issue, in particular for uniaxially prestressed cylindrically-shaped structures, as well as biaxially prestressed plate-like structures [6].

In the present study, focused on uniaxially prestressed cylindrical structures, the confounding effects of finite dimensions and prestress are further explored, analytically and experimentally. Furthermore, we articulate and evaluate a strategy for decoupling prestress and waveguide effects from estimates of material shear viscoelastic properties.

Here, $k[ω]$ is the complex-valued, frequency-dependent shear wave number and $μ[ω]$ is the complex-valued, frequency-dependent shear modulus, comprised of the shear storage modulus, $μR[ω]$⁠, and the shear loss modulus, $μI[ω]$⁠, such that $μ[ω]=μR[ω]+jμI[ω]$ where $j=−1$⁠. The attenuation rate of the wave as it propagates is governed by the imaginary part of the wavenumber: $Imag[k]$⁠. In a viscoelastic material, both the shear storage and loss moduli affect both the phase speed and attenuation rate. In a purely elastic material, $μI=0$⁠, there is no attenuation and the phase speed is independent of frequency (nondispersive) and reduces to $c=μ/ρ$⁠. While some studies have assumed pure elasticity (no viscosity), often their analyses can be generalized to the linear viscoelastic problem for harmonic motion by simply adding the imaginary shear loss modulus to form the complex shear modulus.

Here u, v, and w refer to the displacement component in the x, y, and z direction, respectively, and subscripted $x, y, z, and t$ after a comma refer to partial derivatives with respect to that spatial or time dimension. The term $κ$ is the bulk modulus of the material, which will affect compression wave behavior, but not shear wave behavior. In biological soft tissue or other nearly incompressible materials, $κ$ is multiple orders of magnitude greater than $μ$ such that the Poisson's ratio for the material approaches, but does not equal 0.5.

Here, $θ$ is the angle between the direction of propagation and the axis of uniaxial stress.

where $I¯1=J−23I1$⁠, with $I1=tr[C]$⁠, $C=FTF$⁠, $J=Det[F]$⁠, and $Jm$ is a limiting parameter for $I¯1−3$⁠. Here, $F$ is the deformation gradient, $C$ is the right Cauchy-Green tensor, and $J=1$ if the material is incompressible. The first term of the strain energy function is based on the isochoric deformation of the isotropic material, and the second term only exists if the material is compressible. The shear modulus is given by $2∂ψ∂I1$⁠, which in the incompressible and small strain limit is equal to $μ$ and is consistent with infinitesimal strain theory. Consequently, Eqs. (5) and (6), governing planar shear wave phase speed should remain valid here for nearly incompressible materials. Utilization of other strain energy functions that work for larger deformations is left for future study. It is also possible to incorporate the prestress field directly into the strain energy function [17,18]. If it is assumed that the superimposed wave motion on top of the finite deformation is small, one will still recover Eqs. (5) and (6) for planar shear wave phase speed.

Here, $J#$ denotes a Bessel function of the first kind of # = 0th, 1st, or 2nd order, and $α$ and $β$ are wavenumbers in the r (radial) direction. The wavenumber in the z (axial) direction, $ξ$⁠, is determined by satisfying a complex frequency equation found in the reference (Eqs. 8.2. and 24-25) [19]. Note, $ξ$ is in all three equations. In the present study our focus is transverse wave motion, which with respect to the polar coordinate system, will have motion polarization in the r and $θ$ direction. Our focus will be to measure the evolution of that motion along the z direction in order to estimate $ξ$⁠. In the present study, we are interested in incorporating uniaxial prestress $σ$ and, if possible, would like to identify simpler analytical expressions than the above analysis affords.

We have two propagating waves (α) in the + or −z direction, and 2 nonpropagating (near field or evanescent) waves (β) in the + or −z direction. For the propagating waves, the phase speed will be: $cph=ωReal[α]$⁠. Taking the limit that $σA≪2EI$ we see that $α=(ωa)1/2$ and thus for the elastic case $cph=ω1/2(EIρA)1/4$⁠, which is the classic thin (Euler–Bernoulli) beam transverse vibration solution. On the other hand, taking the limit of tension $σA≫2EI$ we then drop $EIu,zzzz$ from Eq. (11) and reformulate the solution to find that there are two propagating solutions with $cph=(σρ)1/2$⁠. This is the classic transverse thin-string vibration solution. A case between the extremes of either neglecting $σA$ or $EI$ still does not match the phase speed of bulk shear waves, which is given by Eq. (5).

The thin beam formulation presented in Sec. 2.3 is only a reasonable approximation when the wavelength of the propagating transverse waves is at least an order of magnitude greater than the beam's cross-sectional radius R, or equivalent radius for a noncircular cross-section. A formulation yielding a valid approximation for shorter transverse waves is based on the Timoshenko beam theory, which allows for shear deformation and accounts for rotational inertia.

In this case, we cannot separate the wavenumber solutions into propagating and nonpropagating pairs. Rather, there will be certain prestress-dependent frequency ranges where we have only one propagating pair and others where we have two propagating pairs.

If the dimensions of the cylinder are such that neither thin nor thick beam approximations are sufficiently accurate, one is left to somehow integrate the prestress condition introduced in Sec. 2.1 into Eqs. ((8)–10)) of Sec. 2.2. A novel approach is proposed to accomplish this with reduced mathematical complexity. In previous studies on transversely isotropic materials not under prestress [20–23] the last author has shown that, by distorting the geometry based on direction and polarization-dependent planar phase speeds, one can then solve an equivalent isotropic problem. This approach to the transverse isotropic problem was called transformation elastography. Uniaxial prestress causes a similar, though not identical, direction and polarization dependence of the planar shear wave phase speed, as exhibited in Eqs. (5) and (6). The same approach is adapted to the acoustoelastic problem here. Note, different geometric distortions may be needed depending on whether the wave motion of interest is a slow or a fast shear wave, based on propagation direction and polarization. The wave polarization of interest here is perpendicular to the uniaxial stress direction and this is governed by slow shear waves.

Here, $ϵL$ is the axial strain. When the formulas are expressed in terms of this, it's clear that one does not need to know $σ$ and $μ$ a priori in order to calculate the equivalent dimensions. One only needs know the axial strain. If the compressibility is not negligible Eq. (16) will need some adjustment.

An analytical and numerical case study was conducted to understand the interactions between uniaxial prestress and waveguide behavior, as well as to validate the proposed transformation acousto-elastography (TAE) approach proposed above. Case study geometric and material property values are provided in Table 1.

Based on these values, equations in Secs. 2.3 and 2.4 can be used to calculate the wavenumber $γ$ and wavelength $λ$ for transverse motion based on thin and thick beam theory, respectively. In order to numerically simulate the approach outlined in Sec. 2.5, finite element analysis (FEA) was performed using ansys Mechanical APDL Version 2022 R2. An axisymmetric mixed u-P formulation was used with Solid273 8-node by three circumferential plane generalized axisymmetric elements with individual element side dimensions in the axial and radial direction of 0.5 mm. Nonlinear attributes of the element were enabled and the material was defined by a hyperelastic Gent model, requiring specification of the bulk modulus, the real part of complex shear modulus, a limiting parameter, $Jm$⁠, as well as the density. One end of the cylinder was fixed and the other end was incrementally displaced in the z-direction, solving the problem in steps until the desired end displacement was reached that resulted in a uniform uniaxial strain throughout the model.

Once the static analysis was done, the solution routine was exited and then reentered this time specifying a harmonic solution routine using the modified stiffness matrix acquired at the end of the static solution routine and the frequency provided in Table 1 (600 Hz). Now, harmonic x-direction displacements were applied to the nodes either on the entire end of the cylinder or just on the nodes at the outer radial edge of the end of the cylinder representing a rigid transverse x-polarized oscillating end or ring in contact with the end of the phantom (see Fig. 1). In addition to the material properties specified for the static analysis, it is also necessary to specify viscous (beta) damping, which was done by taking the ratio of the imaginary to the real part of the shear modulus in Table 1 and dividing by $2πf$ where $f$ is the harmonic driving frequency.

The in-phase steady-state response over a sagittal slice is shown in Fig. 2 for selected cases in Tables 1 and 2. For the case that the entire end is transversely oscillated (upper row in Fig. 2), the response pattern is uniform and a line profile can be taken along the axial direction at any radial position to estimate the complex wavenumber $γ$ by fitting the profile to $Ae−jγz$⁠. For the case that the outer radial edge at the end has oscillated (lower row in Fig. 2), the response is more complex near the excitation source, but far enough away from the source begins to resemble the simpler case. For this case, the complex wavenumber $γ$ is estimated by fitting $Ae−jγz$ to axial line profiles near the outer edge $(r=R)$⁠. Once $γ$ is estimated, the wavelength $λ=2π/Real[γ]$⁠. The TAE-adjusted values for these, $γN/TAE$ and $λN/TAE$⁠, are determined by first distorting the axial z dimension according to Eq. (15). Specifically, the axial location of the nodes (which have been extended already by $1+ϵL$⁠) is then divided by $1+32ϵL$⁠. This alters both the wavenumber and wavelength.

Table 2 summarizes results obtained for 0, 5, 10, and 20% axial prestrain cases. The first four rows of the table provide the complex shear modulus value used in the thin and thick beam theoretical calculations and in the finite element simulations, the associated shear wavenumber $ksh$⁠, and the calculated shear wavelengths along the axis $λs0°$ and perpendicular to the axis $λs90°$⁠, based on Eq. (5). Next, estimates of propagating wavenumber and wavelength based on prestressed thin beam theory (Sec. 2.3), $γThin$ and $λThin$⁠, are provided, followed by estimates based on prestressed thick beam theory (Sec. 2.4), $γThick$ and $λThick$⁠. Thin beam theory wavelengths and wavenumbers are in error by more than a factor of 2, as compared to the true values in the first two rows of Table 2. For thick beam theory, there are two pairs of propagating wavenumber and wavelength values at this frequency. For the pairs that are closer to truth, marked with an asterisk*, the corresponding complex shear modulus is calculated based on Eq. (1) and provided in the line below $λThick$⁠. Its percent difference from the true value of 42.8+j8.05 $kPa$ is also provided, showing a consistent bias of about 9–10% on the high side for the shear storage modulus while being more accurate in estimating the relative value of the shear loss to shear storage modulus. It is expected that the thin or thick beam approximate theories will yield estimates with higher shear moduli, or lower shear wavenumbers, since such simplifying approximations provide greater constraints on how the material will deform, as compared to the full three-dimensional theory.

The next section of the table summarizes results of the numerical FEA simulation coupled with the TAE coordinate transformation that was introduced in Sec. 2.5. The wavenumber and wavelength prior to the coordinate transformation are first provided, $γN$ and $λN$⁠, followed by their values after the TAE coordinate transformation, $γN/TAE$ and $λN/TAE$⁠. The complex shear modulus is subsequently calculated based on Eq. (1) and provided in the line below $λN/TAE$⁠, as well as its percent difference from the true value. Consistently, this approach has outperformed thick beam theory, with percentage differences in shear storage modulus ranging from 1 to 7.5% and excellent estimates of the ratio of shear loss to shear storage modulus. Note the importance of undergoing the TAE coordinate transformation in that the wavenumber before and after transformation significantly diverges with increasing prestrain.

The setup was designed and tested on an Agilent 9.4 Tesla horizontal, 30 cm bore, animal MR system. All fixture parts for the setup were designed in Solidworks (Solidworks 2021) and printed using a fused filament extrusion printer (Prusa Mk3, PRUSA REF) on Polyethylene terephthalate glycol (PETG) to limit interference with signal and keep parts light enough to be moved efficiently by the piezostack. Figure 3 shows the overall design, which prioritizes versatility in sample type, from phantoms, to muscle tissue or large structures that can be excised to fit the dimensions of the gradient coil opening.

Samples are gripped by two clamp types: a fixed clamp at the proximal end, and the tensioner clamp at the distal end. Both clamp types are slotted to allow wooden skewers to penetrate the sample, providing reasonably distributed tension, rather than surface tension that would be only available with typical clamps. The distal clamp is then attached to a reinforced nylon wire that is fed out the back of the magnet bore and attached to a pulley system where adjustable weight applies appropriate tension.

Simultaneous harmonic loading is achieved through a piezostack (P842.10, PI-USA LP) located proximally to the entry of the bore to limit interference from wiring, and attached to a Delrin counterweight to force transference of motion toward the sample. The sample is placed in a cylindrical tube, the actuator, which is attached directly to the piezo-and slotted for adjustable rings that uniformly and circumferentially grip the sample, providing a fixed location of the wave source that is 15 mm thick. While, schematically it appears that this setup would provide harmonic vibration that is polarized along the z-axis, the axis of the cylindrical phantom and of the tensile load, because of any small nonaxisymmetries, for example, due to gravity, the resulting excitation has both z-axis and orthogonal x-axis components.

The phantom was made with Ecoflex™ 00-30 (Smooth-On, Inc., Macungie, PA) platinum-catalyzed silicone. Ecoflex silicone bases were combined in a 1 A:1B ratio and cured at room temperature with minimum shrinkage. Air bubbles produced during the mixing process were removed using a vacuum chamber (5305-1212, Thermo Scientific-Nalgene, Rochester, NY) for 5 to 10 min to minimize porosity and nonhomogeneous bulk material properties.

SLIM MRE [24] provided motion encoding in three orthogonal directions, simultaneously. Sequence parameters were TR/TE = 1600/25 ms with 8-time steps and corresponding 180 deg offsets. A 250 mT/m MEG was used with ten cycles. The data matrix size was 64 × 64 × 40 with a FOV of 48 mm × 48 mm × 30 mm for an isotropic voxel size of 0.75 mm. The piezostack provided an input axial harmonic motion of 10 μm peak amplitude.

Transverse-polarized motion in the sagittal plane is shown in Fig. 4 for the cases of 0, 5, and 10% axial prestrain. The measurements reflect the fact that it is impossible to achieve a perfectly symmetric setup, as in the numerical study. Additionally, some slack in the system or uneven residual stresses may disproportionally affect the system response under 0% prestrain, and lead to the nonmonotonic trend in shear modulus estimates with respect to prestrain, given in Table 2. Images were analyzed the same way that numerical simulation data, shown in the lower half of Fig. 2, was processed, as described in Sec. 2.5, to obtain estimates of the shear wavenumber. Results are included in Table 2 and can be compared directly to the numerical simulation results, except it was not possible with our setup to achieve 20% prestrain.

The wavenumber and wavelength prior to the coordinate transformation are first provided, $γE$ and $λE$⁠, followed by their values after the TAE coordinate transformation, $γE/TAE$ and $λE/TAE$⁠. The complex shear modulus is subsequently calculated based on Eq. (1) and provided in the line below $λE/TAE$⁠. What can be observed is that the coordinate transformation results in wavenumber values that deviate less with prestrain, since the TAE transformation is designed to remove the effect of the prestrain. In that sense, the TAE approach has succeeded. However, a direct comparison to the true value used for the theoretical and numerical studies is not as appropriate since the actual value for the complex shear modulus of the material used in the experiment is not exactly known. The theoretical value used is based on experiments conducted on samples of Ecoflex-30 harmonically vibrated at 600 Hz, but there can be some variation between samples.

The numerical finite element and experimental studies of the previous sections have identified and quantified some of the confounding effects of finite dimensions and nonzero prestress on the elastography approach to estimating material viscoelastic properties in an isotropic cylindrical structure under uniaxial normal stress aligned with the cylinder axis. Additionally, a coordinate transformation approach, transformation acousto-elastography (TAE), has been introduced that enables one to estimate material viscoelastic properties independent of the prestress condition without requiring a priori knowledge of either the viscoelastic properties or stress conditions. Rather, only the amount of deformation, or strain, from the unstressed condition is required. A priori knowledge of prestrain is more easily available than prestress, as medical images of geometry via MRI when doing MR Elastography, or via ultrasound if doing ultrasound elastography, can be acquired under varying prestrain conditions, with precise tracking of anatomical landmarks that will establish prestrain levels during a given measurement. These anatomical images provide no direct information about the corresponding prestress. Once viscoelastic properties are estimated under the known prestrain conditions, then the prestress can also be estimated.

The numerical and experimental studies show the promise of the TAE approach. In the numerical study of Sec. 3, the material shear storage $μR$ and loss $μI$ moduli, as well as the uniaxial normal stress $σ$ were determined within 5% for prestrains up to 10%, and within 7.5% for a prestrain of 20%. Without the TAE adjustment to the complex wavenumber, the error in the estimate of shear modulus would reach 40% at 20% prestrain. Potential sources of error or discrepancies that affect the estimates of both shear moduli and stress are: (1) an incompressibility assumption used to formulate the approach; (2) the assumption of linearity; and (3) ignoring more complex wave patterns that exist in the 3-dimensional structure. Since the bulk modulus is about five orders of magnitude greater than the shear moduli, we believe the assumption of incompressibility is justified. The approach could be reformulated to account for compressibility. With regard to nonlinearity, the numerical study was repeated assuming simple linear elastic material properties and not enabling the geometric nonlinearity feature in ansys. Results were the same, suggesting that, at least in the numerical study, nonlinearity is not a factor for the strain rates considered. The images in Fig. 2 show that, even with the uniform transverse harmonic displacement across the entire end of the cylinder, the steady-state wave pattern shows some variation as a function of radial position. This variation could result in the small variation seen in predictions based on the TAE approach.

In the experimental study, another source of error is that one must consider that the static shear modulus governing the initial deformation due to the prestress is different from the frequency-dependent shear storage modulus that will govern response to harmonic excitation. To overcome this, one possible solution is to conduct studies at multiple frequencies and then fit a rheological model that predicts the static shear storage modulus, as detailed in Sec. 4.

Logical next steps for advancing the strategy introduced here to detangle prestress from material shear stiffness estimates include consideration of more complex geometry and stress conditions, as well as anisotropic and nonuniform material properties. Finite element models based on medical images that can provide detailed geometry and deformation information under varying loading conditions and, if needed, measures of anisotropy and inhomogeneity, may provide a way to advance the TAE technique beyond simple geometries and assumptions of isotropy and homogeneity [25].

• National Science Foundation (NSF) (Grant No. 1852691; Funder ID: 10.13039/100000084).

• National Institutes of Health (NIH) (Grant No. AR071162; Funder ID: 10.13039/100000069).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.