## 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 [2–4]. 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.

Here, $ksh[\omega ]$ is the complex-valued, frequency-dependent shear wave number and $\mu [\omega ]$ is the complex-valued, frequency-dependent shear modulus, comprised of the shear storage modulus, $\mu R[\omega ]$, and the shear loss modulus, $\mu I[\omega ]$, such that $\mu [\omega ]=\mu R[\omega ]+j\mu I[\omega ]$ where $j=\u22121$. The attenuation rate of the wave as it propagates is governed by the imaginary part of the wavenumber: $Imag[ksh[\omega ]]$. In a viscoelastic material, both the shear storage and loss moduli affect both the phase speed and attenuation rate. In a purely elastic material, $\mu I=0$, there is no attenuation and the phase speed is independent of frequency (nondispersive) and reduces to $c=\mu /\rho $. 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.

Here *u*, *v*, and *w* refer to the displacement component in the x, y, and z direction, respectively, and subscripted $x,\u2009y,\u2009z,\u2009and\u2009t$ after a comma refer to partial derivatives with respect to that spatial or time dimension. The term $\kappa $ 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, $\kappa $ is more than a few orders of magnitude greater than $\mu $ such that the Poisson's ratio for the isotropic material approaches, but does not equal 0.5.

Here, $\theta $ is the angle between the direction of propagation and the z-axis.

where $I\xaf1=J\u221223I1$, with $I1=tr[C]$, $C=FTF$, $J=Det[F]$, and $Jm$ is a limiting parameter for $I\xaf1\u22123$. 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\u2202\psi \u2202I1$, which in the incompressible limit is equal to $\mu $ 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.

*R*and thickness

*h*, with free boundary conditions on the top and bottom flat surfaces at $z=\xb1h/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:

*J*denotes a zeroth order Bessel function of the first kind. Also, we define

_{0}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.

*z*motion is governed by Eq. (9), evaluated at $z=0$, and is given by

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 [23–26] 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.

When the formulas are expressed in terms of $\u03f5R$, it is clear that one does not need to know $\sigma $ or $\mu $ 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.

Parameter | Nomenclature | Value(s) |
---|---|---|

Bulk modulus | $\kappa $ | 2 $GPa$ |

Shear storage modulus in plane of isotropy | $\mu R$ | 39.840$\u2009kPa$ |

Limiting parameter | $Jm$ | 50 |

Ratio of shear loss to storage moduli | $\eta =\mu I\u2009/\mu R$ | 0.1642 |

Undeformed plate thickness | $H$ | 8.5 mm |

Undeformed phantom radius | $R$ | 25.4 mm |

Biaxial tensile strain | $\u03f5R$ | 0, .025, .05, 0.1, 0.2 |

Density | $\rho $ | 1070 $kg/m3$ |

Frequency | $f$ | $300\u2009Hz$ |

Parameter | Nomenclature | Value(s) |
---|---|---|

Bulk modulus | $\kappa $ | 2 $GPa$ |

Shear storage modulus in plane of isotropy | $\mu R$ | 39.840$\u2009kPa$ |

Limiting parameter | $Jm$ | 50 |

Ratio of shear loss to storage moduli | $\eta =\mu I\u2009/\mu R$ | 0.1642 |

Undeformed plate thickness | $H$ | 8.5 mm |

Undeformed phantom radius | $R$ | 25.4 mm |

Biaxial tensile strain | $\u03f5R$ | 0, .025, .05, 0.1, 0.2 |

Density | $\rho $ | 1070 $kg/m3$ |

Frequency | $f$ | $300\u2009Hz$ |

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 $\mu R$ and bulk $\kappa $ moduli (specified in ansys as a compressibility factor $2/\kappa $), 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 $\rho $ and shear loss modulus $\mu I$, which is specified via beta damping in ansys, where the beta damping term will equal the ratio of shear loss to storage moduli, $\eta $, divided by the frequency in radians/second, $2\pi 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.

*H*, again using a least squares optimization algorithm. From $ksh/TAE$ an estimate of the material's complex shear modulus, $\mu R+j\mu I$, is found since

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, $\mu R+j\mu 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 $\mu R+j\mu I=39.840+j6.542\u2009kPa$, the prediction at 0% prestrain is $40.1+j6.78\u2009kPa$. 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.

Prestrain | 0% | 2.5% | 5% | 10% | 20% |
---|---|---|---|---|---|

$\lambda sh0\xb0\u2009(mm)$ | 20.408 | 20.15 | 19.89 | 19.36 | 18.35 |

$\lambda sh90\xb0\u2009(mm)$ | 20.408 | 20.66 | 20.91 | 21.40 | 22.36 |

$RE\u2009(mm)$ | 25.4 | 25.1 | 24.9 | 24.5 | 24.1 |

$HE\u2009(mm)$ | 8.50 | 8.41 | 8.36 | 8.40 | 9.33 |

$kRL\u2009(m\u22121)$ | 368−j24.9 | 355−j24.3 | 343−j23.8 | 322−j22.9 | 285−j21.3 |

$\lambda RL\u2009(mm)$ | 17.1 | 17.7 | 18.3 | 19.5 | 22.0 |

$kRL/TAE\u2009(m\u22121)$ | 368−j24.9 | 368−j25.2 | 368−j25.5 | 367−j26.1 | 361−j26.9 |

$\lambda RL/TAE\u2009(mm)$ | 17.1 | 17.1 | 17.1 | 17.1 | 17.4 |

$ksh/TAE\u2009(m\u22121)$ | 305−j25.6 | 304−j25.9 | 304−j26.3 | 302−j26.9 | 304−j27.6 |

$\lambda sh90\xb0/TAE\u2009(mm)$ | 20.6 | 20.7 | 20.7 | 20.8 | 20.7 |

$\mu R+j\mu I/TAE\u2009(kPa)$ | 40.1+j6.78 | 40.2+j6.90 | 40.3+j7.03 | 40.6+j7.27 | 40.3+j7.38 |

$\mu $ Estimated error (%) | 0.78+j0.47 | 1.02+j0.74 | 1.39+j1.00 | 2.15+j1.48 | 1.34+j1.89 |

$\sigma TAE\u2009(kPa)$ | 0 | 6.06 | 12.2 | 24.7 | 49.3 |

$\sigma FEA\u2009(kPa)$ | 0 | 6.03 | 12.1 | 24.4 | 48.3 |

$\sigma $ Estimated error (%) | 0 | 0.49 | 0.82 | 1.21 | 2.03 |

Prestrain | 0% | 2.5% | 5% | 10% | 20% |
---|---|---|---|---|---|

$\lambda sh0\xb0\u2009(mm)$ | 20.408 | 20.15 | 19.89 | 19.36 | 18.35 |

$\lambda sh90\xb0\u2009(mm)$ | 20.408 | 20.66 | 20.91 | 21.40 | 22.36 |

$RE\u2009(mm)$ | 25.4 | 25.1 | 24.9 | 24.5 | 24.1 |

$HE\u2009(mm)$ | 8.50 | 8.41 | 8.36 | 8.40 | 9.33 |

$kRL\u2009(m\u22121)$ | 368−j24.9 | 355−j24.3 | 343−j23.8 | 322−j22.9 | 285−j21.3 |

$\lambda RL\u2009(mm)$ | 17.1 | 17.7 | 18.3 | 19.5 | 22.0 |

$kRL/TAE\u2009(m\u22121)$ | 368−j24.9 | 368−j25.2 | 368−j25.5 | 367−j26.1 | 361−j26.9 |

$\lambda RL/TAE\u2009(mm)$ | 17.1 | 17.1 | 17.1 | 17.1 | 17.4 |

$ksh/TAE\u2009(m\u22121)$ | 305−j25.6 | 304−j25.9 | 304−j26.3 | 302−j26.9 | 304−j27.6 |

$\lambda sh90\xb0/TAE\u2009(mm)$ | 20.6 | 20.7 | 20.7 | 20.8 | 20.7 |

$\mu R+j\mu I/TAE\u2009(kPa)$ | 40.1+j6.78 | 40.2+j6.90 | 40.3+j7.03 | 40.6+j7.27 | 40.3+j7.38 |

$\mu $ Estimated error (%) | 0.78+j0.47 | 1.02+j0.74 | 1.39+j1.00 | 2.15+j1.48 | 1.34+j1.89 |

$\sigma TAE\u2009(kPa)$ | 0 | 6.06 | 12.2 | 24.7 | 49.3 |

$\sigma FEA\u2009(kPa)$ | 0 | 6.03 | 12.1 | 24.4 | 48.3 |

$\sigma $ Estimated error (%) | 0 | 0.49 | 0.82 | 1.21 | 2.03 |

This is compared to the biaxial stress value $(\sigma 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].

^{TM}, 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 $\mu 0$ in parallel with a fractional order springpot that is defined by $\alpha $ and $\mu \alpha $ such that, at frequency, $\omega =2\pi f\u2009(rad/s)$, the shear storage $\mu R$ and loss $\mu I$ moduli are given by the following equations.

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: $\mu 0=13.3\u2009kPa$, $\mu \alpha =2\u2009kPa\xb7s\alpha $, and $\alpha =1/3$. For Ecoflex-30, measurements conducted over the range of 200 Hz–1 kHz yielded: $\mu 0=27\u2009kPa$, $\mu \alpha =1.5\u2009kPa\xb7s\alpha $, and $\alpha =0.3$. Thus, while Ecoflex-30 is “stiffer” (higher $\mu 0$) under static conditions, due to its lower viscosity the magnitude of its complex shear modulus increases at a slower rate with frequency ($\alpha =0.3$), as compared to that of Ecoflex-10 ($\alpha =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+\u03f5R)$. In the experimental case the radial strain, $\u03f5R$, 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%.)

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 $\lambda RL=2\pi /Real[kRL].$ Next, as in the numerical study both $\lambda 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 $\lambda 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, $\mu R+j\mu 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, $\u03f5R$.

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, $\mu R+j\mu 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 $\mu R+j\mu I=39.840+j6.542\u2009kPa).$ 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 H_{2}O, which induces a radial prestrain of $\u03f5R=7.03%$. Also, it is found that estimates of the radial biaxial prestress, $\sigma 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.

Air pressure (cm H_{2}0) | 0 | 25 | 50 | 75 | 100 |
---|---|---|---|---|---|

Radial strain $\u03f5R\u2009(%)$ | 0 | 0.75 | 2.52 | 4.71 | 7.03 |

$RE\u2009(mm)$ | 25.4 | 25.3 | 25.1 | 24.9 | 24.7 |

$HE\u2009(mm)$ | 8.50 | 8.47 | 8.41 | 8.37 | 8.35 |

$kRL\u2009(m\u22121)$ | 368−j24.9 | 366−j25.1 | 361−j25.6 | 353−j26.4 | 345−j25.8 |

$\lambda RL\u2009(mm)$ | 17.1 | 17.2 | 17.4 | 17.8 | 18.2 |

$kRL/TAE\u2009(m\u22121)$ | 368−j24.9 | 370−j25.4 | 374−j26.6 | 377−j28.2 | 380−j28.4 |

$\lambda RL/TAE\u2009(mm)$ | 17.1 | 17.0 | 16.8 | 16.7 | 16.6 |

$ksh/TAE\u2009(m\u22121)$ | 305−j25.6 | 306−j26.1 | 310−j27.3 | 312−j28.9 | 315−j29.2 |

$\lambda sh90\xb0/TAE\u2009(mm)$ | 20.6 | 20.5 | 20.3 | 20.1 | 19.9 |

$\mu R+j\mu I/TAE\u2009(kPa)$ | 40.1+j6.78 | 39.6+j6.79 | 38.7+j6.87 | 38.0+j7.11 | 37.3+j6.97 |

Estimated difference % | 0.78+j0.47 | 0.47+j0.70 | 2.76+j1.29 | 4.27+j2.12 | 5.99+j2.05 |

$\sigma TAE\u2009(kPa)$ | 0 | 1.78 | 5.85 | 10.7 | 15.7 |

$\sigma 1\u2009(kPa)$ | 0 | 2.52 | 5.70 | 9.16 | 12.8 |

Air pressure (cm H_{2}0) | 0 | 25 | 50 | 75 | 100 |
---|---|---|---|---|---|

Radial strain $\u03f5R\u2009(%)$ | 0 | 0.75 | 2.52 | 4.71 | 7.03 |

$RE\u2009(mm)$ | 25.4 | 25.3 | 25.1 | 24.9 | 24.7 |

$HE\u2009(mm)$ | 8.50 | 8.47 | 8.41 | 8.37 | 8.35 |

$kRL\u2009(m\u22121)$ | 368−j24.9 | 366−j25.1 | 361−j25.6 | 353−j26.4 | 345−j25.8 |

$\lambda RL\u2009(mm)$ | 17.1 | 17.2 | 17.4 | 17.8 | 18.2 |

$kRL/TAE\u2009(m\u22121)$ | 368−j24.9 | 370−j25.4 | 374−j26.6 | 377−j28.2 | 380−j28.4 |

$\lambda RL/TAE\u2009(mm)$ | 17.1 | 17.0 | 16.8 | 16.7 | 16.6 |

$ksh/TAE\u2009(m\u22121)$ | 305−j25.6 | 306−j26.1 | 310−j27.3 | 312−j28.9 | 315−j29.2 |

$\lambda sh90\xb0/TAE\u2009(mm)$ | 20.6 | 20.5 | 20.3 | 20.1 | 19.9 |

$\mu R+j\mu I/TAE\u2009(kPa)$ | 40.1+j6.78 | 39.6+j6.79 | 38.7+j6.87 | 38.0+j7.11 | 37.3+j6.97 |

Estimated difference % | 0.78+j0.47 | 0.47+j0.70 | 2.76+j1.29 | 4.27+j2.12 | 5.99+j2.05 |

$\sigma TAE\u2009(kPa)$ | 0 | 1.78 | 5.85 | 10.7 | 15.7 |

$\sigma 1\u2009(kPa)$ | 0 | 2.52 | 5.70 | 9.16 | 12.8 |

Predictions based on the TAE concept applied to experimental data to estimate material shear viscoelastic property values at 300 Hz, $\mu R+j\mu 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 $(\sigma 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.

Air pressure (cm H_{2}0) | 0 | 25 | 50 | 75 | 100 |
---|---|---|---|---|---|

Radial strain $\u03f5R\u2009(%)$ | 0 | 6.42 | 8.35 | 9.97 | 11.3 |

$RE\u2009(mm)$ | 25.4 | 24.8 | 24.7 | 24.6 | 24.5 |

$HE\u2009(mm)$ | 8.5 | 8.35 | 8.36 | 8.40 | 8.44 |

$kRL\u2009(m\u22121)$ | 348−j30.4 | 333−j39.2 | 324−j40.3 | 319−j37.8 | 312−j30.0 |

$\lambda RL\u2009(mm)$ | 18.1 | 18.8 | 19.4 | 19.7 | 20.1 |

$kRL/TAE\u2009(m\u22121)$ | 348−j30.4 | 364−j42.8 | 362−j45.1 | 363−j43.0 | 361−j34.7 |

$\lambda RL/TAE\u2009(mm)$ | 18.1 | 17.3 | 17.4 | 17.3 | 17.4 |

$ksh/TAE\u2009(m\u22121)$ | 283−j31.3 | 299−j44.0 | 297−j46.4 | 299−j44.3 | 298−j35.7 |

$\lambda sh90\xb0/TAE\u2009(mm)$ | 22.2 | 21.0 | 21.1 | 21.0 | 21.2 |

$\mu R+j\mu I/TAE\u2009(kPa)$ | 45.6+j10.2 | 39.8+j12.0 | 40.0+j12.8 | 39.9+j12.1 | 41.3+j10.1 |

Estimated difference % | 15.6+j6.63 | 2.09+j13.3 | 2.88+j15.2 | 2.48+j13.6 | 5.09+j8.07 |

$\sigma TAE\u2009(kPa)$ | 0 | 15.3 | 20.0 | 23.9 | 28.0 |

Air pressure (cm H_{2}0) | 0 | 25 | 50 | 75 | 100 |
---|---|---|---|---|---|

Radial strain $\u03f5R\u2009(%)$ | 0 | 6.42 | 8.35 | 9.97 | 11.3 |

$RE\u2009(mm)$ | 25.4 | 24.8 | 24.7 | 24.6 | 24.5 |

$HE\u2009(mm)$ | 8.5 | 8.35 | 8.36 | 8.40 | 8.44 |

$kRL\u2009(m\u22121)$ | 348−j30.4 | 333−j39.2 | 324−j40.3 | 319−j37.8 | 312−j30.0 |

$\lambda RL\u2009(mm)$ | 18.1 | 18.8 | 19.4 | 19.7 | 20.1 |

$kRL/TAE\u2009(m\u22121)$ | 348−j30.4 | 364−j42.8 | 362−j45.1 | 363−j43.0 | 361−j34.7 |

$\lambda RL/TAE\u2009(mm)$ | 18.1 | 17.3 | 17.4 | 17.3 | 17.4 |

$ksh/TAE\u2009(m\u22121)$ | 283−j31.3 | 299−j44.0 | 297−j46.4 | 299−j44.3 | 298−j35.7 |

$\lambda sh90\xb0/TAE\u2009(mm)$ | 22.2 | 21.0 | 21.1 | 21.0 | 21.2 |

$\mu R+j\mu I/TAE\u2009(kPa)$ | 45.6+j10.2 | 39.8+j12.0 | 40.0+j12.8 | 39.9+j12.1 | 41.3+j10.1 |

Estimated difference % | 15.6+j6.63 | 2.09+j13.3 | 2.88+j15.2 | 2.48+j13.6 | 5.09+j8.07 |

$\sigma TAE\u2009(kPa)$ | 0 | 15.3 | 20.0 | 23.9 | 28.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 $\mu R$ and loss $\mu 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 $\sigma $ 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 H_{2}O, 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.