Abstract

Dynamic elastography attempts to reconstruct quantitative maps of the viscoelastic properties of biological tissue, properties altered by disease and injury, by noninvasively measuring mechanical wave motion in the tissue. Most reconstruction strategies that have been developed neglect boundary conditions, including quasi-static tensile or compressive loading resulting in a nonzero prestress. Significant prestress is inherent to the functional role of some biological tissues, such as skeletal and cardiac muscle, arterial walls, and the cornea. In the present article a novel configuration, inspired by corneal elastography but generalizable to other applications, is studied. A polymer phantom layer is statically elongated via an in-plane biaxial normal stress while the phantom's response to transverse vibratory excitation is measured. We examine the interplay between biaxial prestress and waveguide effects in this plate-like tissue phantom. Finite static deformations caused by prestressing coupled with waveguide effects lead to results that are predicted by a novel coordinate transformation approach previously used to simplify reconstruction of anisotropic properties. Here, the approach estimates material viscoelastic properties independent of the nonzero prestress conditions without requiring advanced knowledge of those stress conditions.

1 Introduction

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 initial studies have focused on larger organs, such as the liver or brain, where boundary effects were assumed negligible. But, as elastography expands 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 vessels, 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 and Lamb wave dispersion [24]. Limited studies on cardiac elastography have also acknowledged the frequency-dependent (i.e., wavelength-dependent) waveguide behavior of the heart wall [5,6].

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 has 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 uni-axially prestressed cylindrically-shaped structures, as well as biaxially prestressed plate-like structures [7]. In another article, the impact of neglecting the prestress effect in cornea elastography is quantified [8].

In the present study, focused on biaxially prestressed plate-like 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.

2 Theory

2.1 Transverse Wave Motion in a Prestressed Viscoelastic Material.

Dynamic elastography methods are based on the assumption that the measured transverse wave speed or wavelength for small amplitude (linear theory) motion is directly related to the material's elastic or viscoelastic properties [9]. Assuming isotropy, homogeneity, and neglecting boundary effects or variation in density ρ in a viscoelastic material, the frequency-dependent shear wave phase speed, c[ω], for harmonic excitation at circular frequency, ω, is
c=ω/Real[ksh[ω]]=1/ρReal[1/μ[ω]]
(1)

Here, ksh[ω] 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[ksh[ω]]. 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 are generalizable to the linear viscoelastic problem for harmonic motion by adding the imaginary shear loss modulus to form the complex shear modulus.

Consider the introduction of a biaxial static prestress, σ=σî+σĵ, that exists in the xy or rφ plane, perpendicular to the z axis, as illustrated in Fig. 1. Both Cartesian (x,y,z) and polar (r,φ,z) coordinate systems are used in this study. If the static deformation due to the prestress is assumed to be small or incremental, infinitesimal strain theory can be used to incorporate σ into the equations of motion leading to the following [10,11]
ρu,tt=(κ+4μ3)u,xx+μu,yy+(μσ2)u,zz+(κ+μ3)v,xy+(κ+μ3+σ2)w,xz
(2)
ρv,tt=(κ+4μ3)v,yy+μv,xx+(μσ2)v,zz+(κ+μ3)u,xy+(κ+μ3+σ2)w,yz
(3)
ρw,tt=(κ+4μ3)w,zz+(μ+σ2)w,yy+(μ+σ2)w,xx+(κ+μ3σ2)u,xz+(κ+μ3σ2)v,yz
(4)
Fig. 1
Plate of thickness h and radius R under biaxial tensile stress σ subject to harmonic z-direction displacement wzRejωt
Fig. 1
Plate of thickness h and radius R under biaxial tensile stress σ subject to harmonic z-direction displacement wzRejωt
Close modal

Here u, v, and w refer to the displacement component in the x, y, and z direction, respectively, and subscripted x,y,z,andt 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 more than a few orders of magnitude greater than μ such that the Poisson's ratio for the isotropic material approaches, but does not equal 0.5.

The prestress alters the otherwise isotropic, direction-invariant, nature of the medium. Taking the case of harmonic excitation at circular frequency, ω, solving the above equations leads to expressions for two shear waves where motion polarization is perpendicular to the direction of wave propagation. For one of the shear waves, referred to here as the slow shear wave with phase speed in the elastic case given as cs, the motion polarization, and direction of propagation, while perpendicular to each other, both lie within the x-y plane. The other shear wave, referred to here as the fast shear wave cf, will have its propagation direction and polarization direction forming a plane aligned with the z-axis and perpendicular to the x-y plane. Expressions for these phase speeds squared are as follows:
cs2=μρ
(5)
cf2[θ]=μρ(1σ2μcos[2θ])
(6)

Here, θ is the angle between the direction of propagation and the z-axis.

The above analysis based on infinitesimal strain theory loses accuracy as the strain caused by σ becomes significant, in other words beyond a few percent. As prestrain increases it may be necessary to use finite strain theory, also known as large deformation theory, to account for changed geometry, and with it use a hyperelastic model of the material properties that may introduce material nonlinearity into the static analysis. If the harmonic waves imposed upon the finitely deformed medium are, themselves, of small amplitude oscillatory motion, it may still be reasonable to use a linearized analysis to describe wave motion, with linearized parameter values dependent on the degree of static deformation. In the analysis of small amplitude wave motion imposed upon finitely deformed hyperelastic materials, it is common to employ a finite strain energy function ψ to describe the material's properties as static prestress or strain is applied. Many different material models and functions have been introduced and used to describe the behavior of nearly incompressible materials, like soft biological tissue [1217]. For the present study, limiting our focus to biaxial prestress of a nearly incompressible isotropic material and considering finite strains of no more than 20%, we will use a Gent hyperelastic model, which can be defined using the following finite strain energy function [18]
ψ=μJm2ln[1I¯13Jm]+κ2(J212ln[J])
(7)

where I¯1=J23I1, with I1=tr[C], C=FTF, J=Det[F], and Jm is a limiting parameter for I¯13. 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 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 is left for future study. It is also possible to incorporate the prestress field directly into the strain energy function [19,20]. If it is assumed that the superimposed wave motion on top of the finite deformation is small, one may still recover Eqs. (5) and (6) for planar shear wave phase speed.

2.2 Two-Dimensional Waveguide With Polar Wave Front Without Prestress.

We adapt Meral et al. [21], with the geometry shown in Fig. 1. We have a cylindrical homogeneous, isotropic viscoelastic specimen of radius R and thickness h, with free boundary conditions on the top and bottom flat surfaces at z=±h/2. It is axisymmetrically driven in the axial “z” direction at its curved outer radial surface, r = R, with a harmonic displacement (nonhomogeneous boundary condition) of amplitude wzR such that:
wz(R,z,t)=wzRejωt
(8)
We have the following steady-state solution for the antisymmetric wave case depicting Rayleigh–Lamb wave behavior throughout the phantom:
wz(r,z,t)=Aα¯β¯2α¯2{cos(α¯z)2kRL2kRL2β¯2cos(α¯h/2)cos(β¯h/2)cos(β¯z)}J0(kRLr)ejωt
(9)
Here, A is an unknown constant to be determined based on boundary conditions, and J0 denotes a zeroth order Bessel function of the first kind. Also, we define
α¯2=kp2kRL2
(10)
β¯2=ksh2kRL2
(11)
kp=ω/(λ+2μ)/ρ
(12)
ksh=ω/μ/ρ
(13)
Here, kp, ksh, and kRL are the compression, shear and Rayleigh–Lamb wavenumbers. Possible solutions for kRL are obtained by satisfying the following equation [22]
tan(β¯h/2)tan(α¯h/2)=ksh44kRL2α¯β¯
(14)

Applying the nonhomogeneous boundary condition, Eq. (9), we expect only the lowest order solution from Eq. (14) to be mechanically driven, and the corresponding wavenumber kRL will be close in value to, but slightly larger than ksh, associated with a slightly slower wave propagation speed and shorter wavelength.

On the central plane of the plate normal z motion is governed by Eq. (9), evaluated at z=0, and is given by
wz(r,z=0,t)=Akp2kRL2ksh2kp2{ksh2ksh22kRL2}J0(kRLr)ejωt=wzaJ0(kRLr)J0(kRLR)ejωt
(15)

Note, in this equation the only “r” dependence is within J0(kRLr).

2.3 Accounting for in-Plane Biaxial Stress Using Transformation Acousto-Elastography.

A novel approach is proposed to account for biaxial prestress aligned in the plane of the plate. In previous studies on transversely isotropic materials not under prestress [2326] 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 transversely isotropic problem was called “Transformation Elastography”. Biaxial prestress causes a similar, though not identical, direction and polarization dependence of the planar shear wave phase speed, as exhibited in particular in Eq. (6). The same approach is adapted to the acoustoelastic problem here. Note, different geometric distortions are needed depending on whether the wave motion of interest is a slow or a fast shear wave, based on propagation direction and polarization. The out-of-plane axisymmetric wave motion detailed in the previous subjection matches the criteria for a fast shear wave.

Referring to Eq. (6), the phase speed in the direction of the prestress (radial direction) divided by the phase speed in the unstressed medium is cf[θ=90°]μ/ρ=1+σ2μ and for propagation perpendicular to the stress plane is cf[θ=0°]μ/ρ=1σ2μ. Approximating the material as incompressible, the prestress results in a static strain changing its radius to R(1+σ6μ)=R(1+ϵR) and its thickness to H/(1+ϵR)2. Here, ϵR is the radial strain due to the biaxial prestress acting in the radial direction. The equivalent isotropic unstressed plate of thickness He and radius Re needs to have its dimensions adjusted so that the planar shear wave phase speed is μ/ρ in all directions. This results in:
Re=R(1+σ6μ)/1+σ2μ=R(1+ϵR)/1+3ϵR
(16)
He=H/(1+ϵR)21σ2μ=H/(1+ϵR)213ϵR
(17)

When the formulas are expressed in terms of ϵR, it is clear that one does not need to know σ or μ a priori in order to calculate the equivalent dimensions. One only needs know the radial strain. If the compressibility is non-negligible, Eqs. ((16) and (17)) will need some adjustment.

3 Numerical Case Study—Transformation Acousto-Elastography Validation

A numerical case study was conducted to understand the interactions between biaxial 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. These material property and geometry values are chosen to match those for the material and fixture used in experimental studies described later in this article.

Table 1

Parameter values for case studies

ParameterNomenclatureValue(s)
Bulk modulusκ2 GPa
Shear storage modulus in plane of isotropyμR39.840kPa
Limiting parameterJm50
Ratio of shear loss to storage moduliη=μI/μR0.1642
Undeformed plate thicknessH8.5 mm
Undeformed phantom radiusR25.4 mm
Biaxial tensile strainϵR0, .025, .05, 0.1, 0.2
Densityρ1070 kg/m3
Frequencyf300Hz
ParameterNomenclatureValue(s)
Bulk modulusκ2 GPa
Shear storage modulus in plane of isotropyμR39.840kPa
Limiting parameterJm50
Ratio of shear loss to storage moduliη=μI/μR0.1642
Undeformed plate thicknessH8.5 mm
Undeformed phantom radiusR25.4 mm
Biaxial tensile strainϵR0, .025, .05, 0.1, 0.2
Densityρ1070 kg/m3
Frequencyf300Hz

An axisymmetric finite element (FE) nonlinear static deformation study followed by a time-harmonic analysis based on linearization of the stiffness matrix after deformation was conducted using ansys Mechanical APDL 2020 revision 2 (ansys, Canonsburg, PA). A plate-shaped structure with radius R = 25.4 mm and thickness H = 8.5 mm was defined using Plane183 eight-node elements with individual element side dimensions of ∼0.1 mm. This resulted in 16,341 nodes and 5,334 elements. (Use of different element dimensions confirmed resolution was sufficient for these studies.) A mixed u-P (displacement-pressure) solution approach was used due to the near-incompressibility of the element material properties, given in Table 1, in order to avoid locking, which can occur when using a displacement formulation that is dependent on volumetric strain calculations that are difficult to accurately calculate when Poisson's ratio is near 0.5. In this formulation, as opposed to a more efficient formulation based only on displacements, the hydrostatic pressure or volume change rate is interpolated on the element level and solved on the global level independently in the same way as displacements. The material elasticity was defined using a Gent model, which only requires specification of the shear storage μR and bulk κ moduli (specified in ansys as a compressibility factor 2/κ), as well as the parameter Jm. With the nonlinear deformation analysis enabled, a radial displacement boundary condition was applied to the nodes on the outer edge of the model, with displacements in the r direction of 0, 0.635, 1.27, 2.54, and 5.08 mm in five separate simulations, resulting in percentage strains in the radial direction of 0%, 2.5%, 5%, 10%, and 20%, respectively. After application of each deformation, a time-harmonic perturbation analysis was performed by applying a displacement to the same nodes on the outer edge of the phantom, but in the out-of-plane z direction with amplitude of 10 microns and at a frequency f =300 Hz. For the harmonic (dynamic) analysis, in addition to shear storage and bulk moduli, it is also necessary to specify density ρ and shear loss modulus μI, which is specified via beta damping in ansys, where the beta damping term will equal the ratio of shear loss to storage moduli, η, divided by the frequency in radians/second, 2πf. The static displacement and stress conditions due to the initial radial deformation, and the dynamic displacement, due to the harmonic excitation, were output for every node into text files, which were imported into matlab for analysis.

In matlab, the “lsqcurvefit” command was used to optimize the fit of Eq. (15) to the FEA-generated out-of-plane harmonic displacement along a radial line. Specifically, it was necessary to identify best-fit values for complex-valued kRL and wza that result in Eq. (15) matching the in-phase (real part) and out-of-phase (imaginary part) simulated response along a radial line. Once kRL is determined the Rayleigh–Lamb wavelength λRL=2π/Real[kRL]. Next, both λRL and kRL can be scaled to estimate the undeformed case via TAE by a distortion of the dimensions as given in Eqs. ((16) and (17)). Specifically
λRL/TAE=λRLReR(1+ϵR)=λRL/1+3ϵR
(18)
kRL/TAE=kRLR(1+ϵR)Re=kRL1+3ϵR
(19)
The complex-valued shear wave number ksh/TAE and the shear wavelength λsh/TAE in the undeformed state then are estimated via solution of Eq. (14) with He from Eq. (17) in place of H, again using a least squares optimization algorithm. From ksh/TAE an estimate of the material's complex shear modulus, μR+jμI, is found since
ksh=2πfρ/(μR+jμI)
(20)

In summary, we have articulated here a strategy to estimate the material's shear storage and loss moduli based on elastography measurements acquired in the presence of a known static deformation, but potentially unknown prestress condition and unknown material moduli. Equations ((16) and (17)) are expressed in terms of the prestrain, independent of the prestress and material properties.

Predictions based on the TAE concept to estimate material shear viscoelastic property values at 300 Hz, μR+jμI, independent of prestress conditions are in Table 2. The case of no prestress (0% prestrain) serves as a check on the accuracy of the finite element model and matlab -based least squares curve fit in the process given above in this section. In other words, this level of error in estimate is independent of the TAE approach. Whereas the actual value of the complex shear modulus is μR+jμI=39.840+j6.542kPa, the prediction at 0% prestrain is 40.1+j6.78kPa. The “Estimate Error” shown in the table is the difference of the calculated and actual complex moduli, divided by the actual complex modulus. It is 0.78% and 0.47% for estimates of the shear storage and lost moduli, respectively. For 2.5%, 5%, 10%, and 20% radial prestrain cases, the TAE approach yields calculations of the complex shear modulus with the error in the shear storage and loss moduli predictions usually within 2%, except reaching 2.15% in the case of 10% prestrain.

Table 2

Shear storage and loss moduli estimates with unknown prestress conditions using TAE

Prestrain0%2.5%5%10%20%
λsh0°(mm)20.40820.1519.8919.3618.35
λsh90°(mm)20.40820.6620.9121.4022.36
RE(mm)25.425.124.924.524.1
HE(mm)8.508.418.368.409.33
kRL(m1)368−j24.9355−j24.3343−j23.8322−j22.9285−j21.3
λRL(mm)17.117.718.319.522.0
kRL/TAE(m1)368−j24.9368−j25.2368−j25.5367−j26.1361−j26.9
λRL/TAE(mm)17.117.117.117.117.4
ksh/TAE(m1)305−j25.6304−j25.9304−j26.3302−j26.9304−j27.6
λsh90°/TAE(mm)20.620.720.720.820.7
μR+jμI/TAE(kPa)40.1+j6.7840.2+j6.9040.3+j7.0340.6+j7.2740.3+j7.38
μ Estimated error (%)0.78+j0.471.02+j0.741.39+j1.002.15+j1.481.34+j1.89
σTAE(kPa)06.0612.224.749.3
σFEA(kPa)06.0312.124.448.3
σ Estimated error (%)00.490.821.212.03
Prestrain0%2.5%5%10%20%
λsh0°(mm)20.40820.1519.8919.3618.35
λsh90°(mm)20.40820.6620.9121.4022.36
RE(mm)25.425.124.924.524.1
HE(mm)8.508.418.368.409.33
kRL(m1)368−j24.9355−j24.3343−j23.8322−j22.9285−j21.3
λRL(mm)17.117.718.319.522.0
kRL/TAE(m1)368−j24.9368−j25.2368−j25.5367−j26.1361−j26.9
λRL/TAE(mm)17.117.117.117.117.4
ksh/TAE(m1)305−j25.6304−j25.9304−j26.3302−j26.9304−j27.6
λsh90°/TAE(mm)20.620.720.720.820.7
μR+jμI/TAE(kPa)40.1+j6.7840.2+j6.9040.3+j7.0340.6+j7.2740.3+j7.38
μ Estimated error (%)0.78+j0.471.02+j0.741.39+j1.002.15+j1.481.34+j1.89
σTAE(kPa)06.0612.224.749.3
σFEA(kPa)06.0312.124.448.3
σ Estimated error (%)00.490.821.212.03
The subsequent estimate of the biaxial stress using the TAE approach is
σTAE=6ϵRμR
(21)

This is compared to the biaxial stress value (σFEA) determined from the finite element simulation, which is the averaged radial stress in all elements of the model. For the given loading condition, the stress was uniform throughout the model. Here, the percent error in prestress prediction stays within a few percent reaching a maximum of 2.03% when the radial prestrain reaches 20%.

4 Numerical Case II and Experiment: Transformation Acousto-Elastography Adaptation and Evaluation for “Curved Plate”

4.1 Test Configuration and Methods.

To experimentally evaluate TAE-based estimates of material viscoelastic properties, a test configuration similar to, though not the same as the numerical study of Sec. 3 was constructed. This configuration was inspired by earlier studies of our group that identified the confounding effect of intraocular pressure (IOP) on estimates of cornea shear elasticity properties based on optical coherence elastography measurements [8]. It was shown that in-plane biaxial tensile stress in the cornea due to nonzero IOP affected measured Rayleigh–Lamb wave motion induced on the cornea, independent of the cornea shear elastic properties.

The experimental configuration, shown in Fig. 2, includes a circular plate-like phantom 8.5 mm thick and with 25.4 mm radius made of a moldable polymer, Ecoflex-30™. The polymer plate is fixed around its perimeter to a cylindrical container that is vibrated axially to create a radially converging, ideally axisymmetric, flexural wave pattern on the phantom material, similar to the FEA study in Sec. 2. Simultaneously, the cylindrical container is statically pressurized, which both creates some curvature of the phantom and places it under a nonzero axisymmetric multi-axial and nonhomogeneous prestress field. A scanning laser Doppler vibrometer (SLDV; PSV-400, Polytec, Waldbronn, Germany), is used to measure out-of-plane flexural motion on the phantom, as described in previous studies [21,27].

Fig. 2
Experimental setup: (a) pressure-backed phantom plate vibrated by shaker, with SLDV measurement of flexural wave motion on surface; (b) no static pressure; and (c) 75 cm H2O static pressure (note surface bulging)
Fig. 2
Experimental setup: (a) pressure-backed phantom plate vibrated by shaker, with SLDV measurement of flexural wave motion on surface; (b) no static pressure; and (c) 75 cm H2O static pressure (note surface bulging)
Close modal
Our group has characterized the dynamic shear viscoelastic properties of Ecoflex materials under infinitesimal (0% prestrain) conditions over a wide frequency range [28]. Previous publications focused on Ecoflex-10TM, a similar polymer, with greater viscosity, as compared to Ecoflex-30. For harmonic excitation of both materials using small perturbations about the unstressed state, we've found that a “fractional Voigt” rheological model best describes the frequency-dependent shear storage and loss moduli of the material. The 3-parameter fractional Voigt model is comprised of a purely elastic element of strength μ0 in parallel with a fractional order springpot that is defined by α and μα such that, at frequency, ω=2πf(rad/s), the shear storage μR and loss μI moduli are given by the following equations.
μR=μ0+μαωαcos(πα/2)
(22)
μI=μαωαsin(πα/2)
(23)

From the Ref. [28], it was found that for Ecoflex-10, the following values accurately describe shear viscoelastic properties over the range of 200 Hz to 7.75 kHz: μ0=13.3kPa, μα=2kPa·sα, and α=1/3. For Ecoflex-30, measurements conducted over the range of 200 Hz–1 kHz yielded: μ0=27kPa, μα=1.5kPa·sα, and α=0.3. Thus, while Ecoflex-30 is “stiffer” (higher μ0) under static conditions, due to its lower viscosity the magnitude of its complex shear modulus increases at a slower rate with frequency (α=0.3), as compared to that of Ecoflex-10 (α=1/3).

The measured flexural wave motion on the phantom material using the SLDV system is exported and read into matlab for further analysis. Typical experimental measurements of the Rayleigh–Lamb wave motion are shown in Fig. 3. As in the numerical study of Sec. 3, the “lsqcurvefit” command was used to optimize the fit of Eq. (15) to the experimentally-measured out-of-plane harmonic displacement along a radial line of the phantom (averaged from measurements along 40 radial lines spaced 9 degrees apart). Unlike in the previous numerical study, it is recognized that in the experimental case, the curvature of the plate means that the SLDV is not measuring motion exactly normal to the plate surface at all locations. The SLDV “imposes” a planar geometry with a radius R that is constant, not changing with increasing pressure-induced prestrain. To account for this in the TAE analysis the initial least squares curve fit to determine kRL is multiplied by (1+ϵR). In the experimental case the radial strain, ϵR, is estimated based on its correlation to the out-of-plane deformation caused by the static pressure, which is measured using a micrometer. This relationship is established via finite element studies using the phantom geometry. (The experimental measurement of vertical displacement with the micrometer was judged to be accurate to the nearest 0.1 mm with confidence, though the micrometer precision was 0.01 mm. We believe this translates to estimates of radial deformation of the order of 0.1 mm, which in turn provides strain precision and accuracy to 0.4%.)

Fig. 3
SLDV measurement showing the real (in-phase) Rayleigh–Lamb wave on the phantom at 300 Hz with (a)0 and (b) 25 cm H2O static pressure
Fig. 3
SLDV measurement showing the real (in-phase) Rayleigh–Lamb wave on the phantom at 300 Hz with (a)0 and (b) 25 cm H2O static pressure
Close modal

Next, dynamic small amplitude axisymmetric excitation is applied to the phantom along its outer boundary and dynamic displacement profiles along 20 equally angularly spaced diameters were averaged, and then used to identify best-fit values for kRL and wza.

Once kRL is determined, the Rayleigh–Lamb wavelength λRL=2π/Real[kRL]. Next, as in the numerical study both λRL and kRL are scaled to estimate the undeformed case via TAE using Eqs. ((18) and (19)). The complex-valued shear wave number ksh/TAE and the shear wavelength λsh/TAE in the undeformed state then are estimated via solution of Eq. (14) with He from Eq. (17) in place of H, again using a least squares optimization algorithm. From ksh/TAE an estimate of the material's complex shear modulus, μR+jμI, is found using Eq. (20).

In summary, this is a strategy to estimate the material's shear storage and loss moduli based on elastography measurements acquired in the presence of a known static deformation, but unknown prestress conditions and unknown material moduli. Once the material moduli are estimated, then this information can be used to estimate the averaged biaxial stress in the phantom, based on the measured estimate of the radial strain, ϵR.

The approach described in this section to estimate the complex shear modulus and biaxial stress is first tested using a finite element simulation that mimics the experiment by applying a pressure load to one side of the plate that is statically fixed at its perimeter. It is then applied to experimental data.

4.2 Results and Evaluation of Transformation Acousto-Elastography.

Predictions based on the TAE concept applied to the “curved surface” numerical study to estimate material shear viscoelastic property values at 300 Hz, μR+jμI, independent of prestress conditions are in Table 3. (As in the “flat surface” studies summarized in Table 2, the “ground truth” value for the complex shear modulus is μR+jμI=39.840+j6.542kPa). As compared to the flat surface numerical study of Sec. 3, the TAE approach does not do quite as well in enabling estimates of the complex shear modulus independent of the preload condition, with the error reaching 6% when the preload pressure reaches 100 cm H2O, which induces a radial prestrain of ϵR=7.03%. Also, it is found that estimates of the radial biaxial prestress, σTAE, based on Eq. (21), are not directly relatable to the actual prestress condition, as the actual condition is not simply biaxial and is not homogenous; it varies with both radial and axial location. In addition, the prestress field is not simply in the radial r direction or in the local tangential plane of the phantom; there is a significant through-plane predominantly z component, as well as a significant x-z shear stress component. Figure 4 is a plot of the principle stress obtained from the FE static simulation that is predominantly aligned with the radial direction but follows the local tangential plane of the phantom.

Fig. 4
Finite element prediction of the first principle stress under 100 cm H2O static pressure applied on its bottom surface. This is a cross section of the axisymmetrically stressed phantom from its center (left) to its fixed perimeter (right). Color scale is in kPa.
Fig. 4
Finite element prediction of the first principle stress under 100 cm H2O static pressure applied on its bottom surface. This is a cross section of the axisymmetrically stressed phantom from its center (left) to its fixed perimeter (right). Color scale is in kPa.
Close modal
Table 3

TAE estimates for numerical curved surface study

Air pressure (cm H20)0255075100
Radial strain ϵR(%)00.752.524.717.03
RE(mm)25.425.325.124.924.7
HE(mm)8.508.478.418.378.35
kRL(m1)368−j24.9366−j25.1361−j25.6353−j26.4345−j25.8
λRL(mm)17.117.217.417.818.2
kRL/TAE(m1)368−j24.9370−j25.4374−j26.6377−j28.2380−j28.4
λRL/TAE(mm)17.117.016.816.716.6
ksh/TAE(m1)305−j25.6306−j26.1310−j27.3312−j28.9315−j29.2
λsh90°/TAE(mm)20.620.520.320.119.9
μR+jμI/TAE(kPa)40.1+j6.7839.6+j6.7938.7+j6.8738.0+j7.1137.3+j6.97
Estimated difference %0.78+j0.470.47+j0.702.76+j1.294.27+j2.125.99+j2.05
σTAE(kPa)01.785.8510.715.7
σ1(kPa)02.525.709.1612.8
Air pressure (cm H20)0255075100
Radial strain ϵR(%)00.752.524.717.03
RE(mm)25.425.325.124.924.7
HE(mm)8.508.478.418.378.35
kRL(m1)368−j24.9366−j25.1361−j25.6353−j26.4345−j25.8
λRL(mm)17.117.217.417.818.2
kRL/TAE(m1)368−j24.9370−j25.4374−j26.6377−j28.2380−j28.4
λRL/TAE(mm)17.117.016.816.716.6
ksh/TAE(m1)305−j25.6306−j26.1310−j27.3312−j28.9315−j29.2
λsh90°/TAE(mm)20.620.520.320.119.9
μR+jμI/TAE(kPa)40.1+j6.7839.6+j6.7938.7+j6.8738.0+j7.1137.3+j6.97
Estimated difference %0.78+j0.470.47+j0.702.76+j1.294.27+j2.125.99+j2.05
σTAE(kPa)01.785.8510.715.7
σ1(kPa)02.525.709.1612.8

Predictions based on the TAE concept applied to experimental data to estimate material shear viscoelastic property values at 300 Hz, μR+jμI, independent of prestress condition are in Table 4. For pressure-induced radial prestrain cases, the TAE approach yields calculations of the complex shear modulus with the error in the shear storage and loss moduli predictions of no more than 6% and 16%, respectively, over a range of pressure-induced prestrains up to 11.3%. And, as in the numerical studies, the biaxial radial stress (σTAE) is estimated using Eq. (21), though in light of the curved surface numerical study, the unknown actual prestress field for this case is expected to be more complex and not simply biaxial.

Table 4

TAE estimates for experimental study

Air pressure (cm H20)0255075100
Radial strain ϵR(%)06.428.359.9711.3
RE(mm)25.424.824.724.624.5
HE(mm)8.58.358.368.408.44
kRL(m1)348−j30.4333−j39.2324−j40.3319−j37.8312−j30.0
λRL(mm)18.118.819.419.720.1
kRL/TAE(m1)348−j30.4364−j42.8362−j45.1363−j43.0361−j34.7
λRL/TAE(mm)18.117.317.417.317.4
ksh/TAE(m1)283−j31.3299−j44.0297−j46.4299−j44.3298−j35.7
λsh90°/TAE(mm)22.221.021.121.021.2
μR+jμI/TAE(kPa)45.6+j10.239.8+j12.040.0+j12.839.9+j12.141.3+j10.1
Estimated difference %15.6+j6.632.09+j13.32.88+j15.22.48+j13.65.09+j8.07
σTAE(kPa)015.320.023.928.0
Air pressure (cm H20)0255075100
Radial strain ϵR(%)06.428.359.9711.3
RE(mm)25.424.824.724.624.5
HE(mm)8.58.358.368.408.44
kRL(m1)348−j30.4333−j39.2324−j40.3319−j37.8312−j30.0
λRL(mm)18.118.819.419.720.1
kRL/TAE(m1)348−j30.4364−j42.8362−j45.1363−j43.0361−j34.7
λRL/TAE(mm)18.117.317.417.317.4
ksh/TAE(m1)283−j31.3299−j44.0297−j46.4299−j44.3298−j35.7
λsh90°/TAE(mm)22.221.021.121.021.2
μR+jμI/TAE(kPa)45.6+j10.239.8+j12.040.0+j12.839.9+j12.141.3+j10.1
Estimated difference %15.6+j6.632.09+j13.32.88+j15.22.48+j13.65.09+j8.07
σTAE(kPa)015.320.023.928.0

5 Discussion and Conclusions

The numerical finite element and experimental studies of the previous two sections have quantified the confounding effects of finite dimensions, surface curvature and nonzero prestress on the elastography approach to estimating material viscoelastic properties in an isotropic plate-like structures under biaxial normal stress within the plane of the plate and under more complex but still axisymmetric stress conditions. Additionally, a novel coordinate transformation approach, TAE, has been introduced that, in theory, enables one to estimate material viscoelastic properties independent of the nonhomogeneous prestress conditions without requiring a priori knowledge of those stress conditions. Rather, only the amount of deformation, or strain, from the unstressed condition is required. If the prestress field is uniform, the method can then estimate the prestress.

The numerical and experimental studies show the promise of the TAE approach. In the flat surface numerical study, the material shear storage μR and loss μI moduli were determined with accuracy from less than 0.5% up to 2.15% error as prestrain levels increased to 20%. Subsequent estimates of the uniform radial stress field σ were comparably accurate over this prestrain range, reaching 2.03% error at 20% prestrain. The TAE approach does not account for geometric nonlinear effects that increase with deformation and thus, this may be a significant source of the error in estimates of both the shear moduli and the prestress. Another source of error or discrepancy is an incompressibility assumption used to formulate the approach. Because the bulk modulus is about five orders of magnitude greater than the shear moduli, we believe the assumption of incompressibility is justified. Though, the approach could be reformulated to account for compressibility. Accounting for geometric and material nonlinearity, which increase as prestrain increases, is left for future development of the TAE approach. Also left for future numerical simulation studies, is the robustness of the TAE approach when confronted with simulated noise in the measurements.

In the curved surface numerical study, which more closely approximates the experimental study, estimates of the complex shear modulus independent of prestress condition reach 6% error at the highest pressure considered—100 cm H2O, which generated a 7% radial strain. The estimate of a biaxial prestress value is less meaningful in this case, as the finite element study shows that the stress field is not simply biaxial and is not homogeneous. Rather, there is variation as a function of both radial r and axial z position (Fig. 4).

Finally, in the experimental study, one must also consider that the static shear modulus governing the initial deformation due to the prestress, is typically different from the frequency-dependent shear storage modulus that will govern response to harmonic excitation. This, and the fact that there may be some residual stress even in the unstressed condition, is why the deformation induced by the pressure loading is significantly greater than that in the curved surface numerical study. Despite this, estimates of the complex shear modulus under nonzero prestress conditions up to 11% prestrain, generally stayed close to each other and within 10% of the expected value. The value predicted in the unstressed case had the greatest difference and it is thought this may be due to unmodeled residual stresses caused by the assembly of the experimental fixture. Based on finite element simulations, we do not believe that the fact that the SLDV measured motion of the surface of the plate as opposed to the centerline was a major source of discrepancy in estimating wave numbers.

Logical next steps for the strategy introduced here for decoupling prestress from material shear stiffness estimates include consideration of more complex geometry and stress conditions, as well as anisotropic and nonuniform material properties. Conducting studies at multiple frequencies would enable identification of a rheological model that could help discern static (constant) versus frequency-dependent components of the shear storage and loss moduli. Use of torsional excitation of waves polarized in the plane of the plate, would excite the slow shear waves (Eq. (5)), which would provide an additional means of estimating the shear modulus [29]. Finite element models based on medical images that can provide detailed geometry [30], deformation 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.

Funding Data

  • 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).

Conflict of Interests

There is no conflict of interest.

References

1.
Larin
,
K. V.
, and
Sampson
,
D. D.
,
2017
, “
Optical Coherence elastography - OCT at Work in Tissue Biomechanics
,”
Biomed. Opt. Express
,
8
(
2
), pp.
1172
1202
.10.1364/BOE.8.001172
2.
Bernal
,
M.
,
Nenadic
,
I.
,
Urban
,
M. W.
, and
Greenleaf
,
J. F.
,
2011
, “
Material Property Estimation for Tubes and Arteries Using Ultrasound Radiation Force and Analysis of Propagating Modes
,”
J. Acoust. Soc. Am.
,
129
(
3
), pp.
1344
1354
.10.1121/1.3533735
3.
Roy
,
T.
, and
Guddati
,
M. N.
,
2021
, “
Shear Wave Dispersion Analysis of Incompressible Waveguides
,”
J. Acous. Soc. Am.
,
149
(
2
), pp.
972
982
.10.1121/10.0003430
4.
Roy
,
T.
,
Urban
,
M.
,
Xu
,
Y.
,
Greenleaf
,
J.
, and
Guddati
,
M. N.
,
2021
, “
Multimodal Guided Wave Inversion for Arterial Stiffness: Methodology and Validation in Phantoms
,”
Phys. Med. Biol.
,
66
(
11
), p.
115020
.10.1088/1361-6560/ac01b7
5.
Nenadic
,
I. Z.
,
Urban
,
M. W.
,
Mitchell
,
S. A.
, and
Greenleaf
,
J. F.
,
2011
, “
Lamb Wave Dispersion Ultrasound Vibrometry (LDUV) Method for Quantifying Mechanical Properties of Viscoelastic Solids
,”
Phys. Med. Biol.
,
56
(
7
), pp.
2245
2264
.10.1088/0031-9155/56/7/021
6.
Urban
,
M. W.
,
Qiang
,
B.
,
Song
,
P.
,
Nenadic
,
I. Z.
,
Chen
,
S.
, and
Greenleaf
,
J. F.
,
2016
, “
Investigation of the Effects of Myocardial Anisotropy for Shear Wave Elastography Using Impulsive Force and Harmonic Vibration
,”
Phys. Med. Biol.
,
61
(
1
), pp.
365
382
.10.1088/0031-9155/61/1/365
7.
Crutison
,
J.
,
Sun
,
M.
, and
Royston
,
T. J.
,
2022
, “
The Combined Importance of Finite Dimensions, Anisotropy, and Pre-Stress in Acoustoelastography
,”
J. Acoust. Soc. Am.
,
151
(
4
), pp.
2403
2413
.10.1121/10.0010110
8.
Sun
,
M. G.
,
Son
,
T.
,
Crutison
,
J.
,
Guaiquil
,
V.
,
Lin
,
S.
,
Nammari
,
L.
,
Klatt
,
D.
,
Yao
,
X.
,
Rosenblatt
,
M. I.
, and
Royston
,
T. J.
,
2022
, “
Optical Coherence Elastography for Assessing the Influence of Intraocular Pressure on Elastic Wave Dispersion in the Cornea
,”
J. Mech. Behav. Biomed Mater
,
128
, p.
105100
.10.1016/j.jmbbm.2022.105100
9.
Manduca
,
A.
,
Bayly
,
P. J.
,
Ehman
,
R. L.
,
Kolipaka
,
A.
,
Royston
,
T. J.
,
Sack
,
I.
,
Sinkus
,
R.
, and
Van Beers
,
B. E.
,
2021
, “
MR Elastography: Principles, Guidelines, and Terminology
,”
Magn. Reson. Med.
,
85
(
5
), pp.
2377
2390
.10.1002/mrm.28627
10.
Biot
,
M. A.
,
1965
,
Mechanics of Incremental Deformation
,
Wiley
,
New York
.
11.
Singh
,
I.
,
Madan
,
D. K.
, and
Gupta
,
M.
,
2010
, “
Propagation of Elastic Waves in Prestressed Media
,”
J. Appl. Math.
,
2010
, p.
e817680
.10.1155/2010/817680
12.
Gennisson
,
J. L.
,
Renier
,
M.
,
Catheline
,
S.
,
Barriere
,
C.
,
Berco
,
J.
,
Tanter
,
M.
, and
Fink
,
M.
,
2007
, “
Acoustoelasticity in Soft Solids: Assessment of the Nonlinear Shear Modulus With the Acoustic Radiation Force
,”
J. Acoust. Soc. Amer.
,
122
(
6
), pp.
3211
3219
.10.1121/1.2793605
13.
Destrade
,
M.
,
Gilchrist
,
M. D.
, and
Saccomandi
,
G.
,
2010
, “
Third- and Fourth-Order Constants of Incompressible Soft Solids and the Acousto-Elastic Effect
,”
J. Acoust. Soc. Amer.
,
127
(
5
), pp.
2759
2763
.10.1121/1.3372624
14.
Destrade
,
M.
, and
Ogden
,
R. W.
,
2010
, “
On the Third- and Fourth-Order Constants of Incompressible Isotropic Elasticity
,”
J. Acoust. Soc. Am.
,
128
(
6
), pp.
3334
3343
.10.1121/1.3505102
15.
Feng
,
Y.
,
Okamoto
,
R. J.
,
Genin
,
G. M.
, and
Bayly
,
P. V.
,
2016
, “
On the Accuracy and Fitting of Transversely Isotropic Material Models
,”
J. Mech. Behav. Biomed. Mater.
,
61
, pp.
554
566
.10.1016/j.jmbbm.2016.04.024
16.
Caenen
,
A.
,
Knight
,
A. E.
,
Rouze
,
N. C.
,
Bottenus
,
N. B.
,
Segers
,
P.
, and
Nightingale
,
K. R.
,
2020
, “
Analysis of Multiple Shear Wave Modes in a Nonlinear Soft Solid: Experiments and Finite Element Simulations With a Tilted Acoustic Radiation Force
,”
J. Mech. Behav. Biomed. Mater
,
107
, p.
103754
.10.1016/j.jmbbm.2020.103754
17.
Remenieras
,
J. P.
,
Bulot
,
M.
,
Gennisson
,
J. L.
,
Patat
,
F.
,
Destrade
,
M.
, and
Bacle
,
G.
,
2021
, “
Acousto-Elasticity of Transversely Isotropic Incompressible Soft Tissues: Characterization of Skeletal Striated Muscle
,”
Phys. Med. Biol.
,
66
(
14
), p.
145009
.10.1088/1361-6560/ac0f9b
18.
Steck
,
D.
,
Qu
,
J.
,
Kordmahale
,
S. B.
,
Tscharnuter
,
D.
,
Muliana
,
A.
, and
Kameoka
,
J.
,
2019
, “
Mechanical Responses of Ecoflex Silicone Rubber: Compressible and Incompressible Behaviors
,”
J. Appl. Polym. Sci.
,
136
(
5
), p.
47025
.10.1002/app.47025
19.
Ogden
,
R. W.
, and
Singh
,
B.
,
2011
, “
Propagation of Waves in an Incompressible Transversely Isotropic Elastic Solid With Initial Stress: Biot Revisited
,”
J. Mech. Mater. Struct.
,
6
(
1–4
), pp.
453
477
.10.2140/jomms.2011.6.453
20.
Destrade
,
M.
, and
Ogden
,
R. W.
,
2013
, “
On Stress-Dependent Elastic Moduli and Wave Speeds
,”
IMA J. Appl. Math.
,
78
(
5
), pp.
965
997
.10.1093/imamat/hxs003
21.
Meral
,
F. C.
,
Royston
,
T. J.
, and
Magin
,
R. L.
,
2011
, “
Rayleigh-Lamb Wave Propagation on a Fractional Order Viscoelastic Plate
,”
J. Acoust. Soc. Am.
,
129
(
2
), pp.
1036
1045
.10.1121/1.3531936
22.
Graff
,
K. F.
,
1975
,
Wave Motion in Elastic Solids
,
Ohio State University Press, Dover Publications Inc.
,
New York
.
23.
Guidetti
,
M.
, and
Royston
,
T. J.
,
2018
, “
Analytical Solution for Converging Elliptic Shear Wave in a Bounded Transverse Isotropic Viscoelastic Material With Nonhomogeneous Outer Boundary
,”
J. Acoust. Soc. Am.
,
144
(
4
), pp.
2312
2323
.10.1121/1.5064372
24.
Guidetti
,
M.
, and
Royston
,
T. J.
,
2019
, “
Analytical Solution for Diverging Elliptic Shear Wave in Bounded and Unbounded Transverse Isotropic Viscoelastic Material With Nonhomogeneous Inner Boundary
,”
J. Acoust. Soc. Am.
,
145
(
1
), pp.
EL59
EL65
.10.1121/1.5088028
25.
Guidetti
,
M.
,
Caratelli
,
D.
, and
Royston
,
T. J.
,
2019
, “
Converging Super-Elliptic Torsional Shear Waves in a Bounded Transverse Isotropic Viscoelastic Material With Nonhomogeneous Outer Boundary
,”
J. Acoust. Soc. Am.
,
146
(
5
), pp.
EL451
EL457
.10.1121/1.5134657
26.
Royston
,
T. J.
,
2021
, “
Analytical Solution Based on Spatial Distortion for a Time-Harmonic Green's Function in a Transverse Isotropic Viscoelastic Solid
,”
J. Acoust. Soc. Am.
,
149
(
4
), pp.
2283
2291
.10.1121/10.0004133
27.
Brinker
,
S.
,
Kearney
,
S. P.
,
Royston
,
T. J.
, and
Klatt
,
D.
,
2018
, “
Simultaneous Magnetic Resonance and Optical Elastography Acquisitions: Comparison of Displacement Images and Shear Modulus Estimations Using a Single Vibration Source
,”
J. Mech. Behav. Biomed Mater
,
84
, pp.
135
144
.10.1016/j.jmbbm.2018.05.010
28.
Yasar
,
T. K.
,
Royston
,
T. J.
, and
Magin
,
R. L.
,
2013
, “
Wideband MR Elastography for Viscoelasticity Model Identification
,”
Mag. Res. Med.
,
70
(
2
), pp.
479
489
.10.1002/mrm.24495
29.
Torres
,
J.
,
Faris
,
I. H.
,
Callejas
,
A.
,
Reyes-Ortega
,
F.
,
Melchor
,
J.
,
Gonzalez-Andrades
,
M.
, and
Rus
,
G.
,
2022
, “
Torsional Wave Elastography to Assess the Mechanical Properties of the Cornea
,”
Sci. Rep.
,
12
(
1
), p.
8354
.10.1038/s41598-022-12151-2
30.
Crutison
,
J.
, and
Royston
,
T. J.
,
2022
, “
The Design and Application of a Diffusion Tensor Informed Finite Element Model for Exploration of Uniaxially Prestressed Muscle Architecture in Magnetic Resonance Imaging
,”
Eng. Comput.
, epub.10.1007/s00366-022-01690-x