Experimental characterization of high-strain-rate viscoelastic and damage behavior in anisotropic soft materials using laser-induced inertial cavitation

TL;DR

This study uses laser-induced cavitation to measure how anisotropic soft materials behave under ultra-high strain rates, revealing directional mechanical responses and damage mechanisms.

Contribution

A novel experimental method combining laser-induced inertial cavitation and constitutive modeling to quantify anisotropic ultra-high-rate mechanics and damage in soft materials.

Findings

LIC experiments showed anisotropic bubble elongation along fiber directions in PVA hydrogels and chicken breast tissue.

Model fitting extracted ultra-high-rate directional moduli and critical stretch thresholds for damage initiation.

Directional damage models were found essential for accurately capturing experimental observations in anisotropic soft materials.

Abstract

Characterizing soft materials at ultra-high strain rates (> 103 s−1) remains a significant challenge due to their nonlinear, large-deformation, and rate-dependent mechanical behavior. Despite these challenges, understanding material response in this regime is essential for a wide range of engineering and biomedical applications, including laser eye surgery, lithotripsy, and high-rate energy deposition in soft tissues. Laser-induced cavitation (LIC) has recently emerged as a powerful experimental approach for subjecting soft materials to extreme, localized loading rates on the microscale, regimes that are difficult to access using conventional mechanical testing methods. However, most biological tissues are anisotropic, optically opaque, and prone to damage at high strain rates, and their orientation-dependent ultra-high-rate mechanical behavior remains poorly understood.

The objective of this study is to develop and validate an experimental method that is based on laser-induced inertial cavitation and enables quantitative measurement of anisotropic soft materials and biological tissues at ultra-high strain rates, while explicitly accounting for fiber-directional mechanics and damage evolution.

LIC experiments were performed in synthetic anisotropic polyvinyl alcohol (PVA) hydrogels and fresh chicken breast tissue using nanosecond-duration pulsed lasers to generate ultra-high-rate deformation. Simultaneously, cavitation bubble dynamics were captured using ultra-high-speed videography at 1–2 million frames per second. Two nonlinear hyper-viscoelastic constitutive models (a Poynting-Thomson model and a generalized Maxwell model) were developed and integrated with the measured bubble dynamics to quantify fiber-directional material response and investigate rate-dependent damage mechanisms. In addition, complementary quasistatic and oscillatory shear rheometry tests were conducted to provide low-rate mechanical benchmarks and investigate the rate dependency of soft materials’ mechanical behavior.

LIC experiments revealed pronounced anisotropic bubble dynamics, with preferential elongation along fiber directions in both PVA hydrogels and chicken breast tissue. Model-based fitting of major-axis bubble radius-time histories enabled the extraction of effective ultra-high-rate directional moduli and critical stretch thresholds for damage initiation. Comparison between constitutive models demonstrated that while viscoelastic relaxation is necessary to capture high-rate damping, incorporation of directional damage is essential to accurately reproduce experimental observations in soft materials and biological tissues.

Figures

A(i)](#F1)). The vertical temperature gradient promoted ice nucleation and subsequent vertical ice growth in parallel alignment, driving PVA chains to aggregate and self-organize into honeycomb-like porous microstructures ([Fig. 1B(i)](#F1) and [Fig. 1F(i)](#F1)). Meanwhile, the abundant hydroxyl groups (-OH) along the PVA chains formed intermolecular hydrogen bonds upon concentration, and the combined effects of molecular entanglement and hydrogen bonding provided the initial physical cross-linking framework of the hydrogel [[51](#R51)]. After freezing, the frozen block was then immersed in a

a](#F6)–[c](#F6)) and chicken breast tissue ([Fig. 6d](#F6)–[f](#F6)). Each loading orientation (0°, 45°, and 90° relative to the fiber direction) was tested three times. For PVA hydrogels, the mean Young’s moduli are $[eqn]$ , and $[eqn]$ . Thus, the modulus measured along 0° is approximately three times larger than that measured along 90°, confirming strong anisotropy. The PVA response is approximately linear elastic in this regime, and loading–unloading hysteresis is negligible. The resulting averaged stiffness matrix $[eqn]$ of anisotropic PVA hydrogel is:

a](#F9)), the no-damage model over-predicts rebound amplitudes and bubble persistence after the first collapse, indicating excessive elastic energy storage. In contrast, the fully damaged model underestimates resistance during the expansion phase and fails to capture the observed rebound timing. Model 1 can capture the bubble radius versus time, except for near-bubble collapse time points. Among the tested formulations, Model 2, which allows damage to accumulate preferentially during tensile expansion but preserves compressive stiffness during collapse, best reproduces both the expansion peak

a](#F10)), the generalized Maxwell model captures the initial expansion peak and early rebound behavior reasonably well across all damage assumptions. However, similar to the Poynting–Thomson case, the no-damage model overpredicts rebound amplitudes, while the fully damaged model underestimates resistance during expansion. Model 1 improves agreement but still exhibits phase errors near collapse. Model 2 again yields the best overall agreement, reproducing both the amplitude and timing of successive rebounds. These results indicate that, even when viscoelastic dissipation is represented through

Full text

Introduction

1

Understanding and quantifying the mechanical properties of soft materials and biological tissues at high strain rates (> 10^3^ s^−1^) is important to tissue injury diagnostics and various medical and engineering applications. For example, soft polymers are widely used as protective materials to mitigate impact and can be applied to prevent blast-induced traumatic brain injury (bTBI) [1–5] and high-energy ballistic impacts to the torso [6]. In clinical settings, high-energy laser and ultrasound pulses are utilized in microscale and nanoscale surgical procedures, including lithotripsy [7] and histotripsy [8]. In all of these scenarios, soft materials and biological tissues are subjected to extreme, transient deformation at strain rates far exceeding those accessible by conventional mechanical tests. However, most currently available tissue and soft material experimental characterization methods are hard to achieve such high strain-rate regimes [9]. For instance, conventional uniaxial [10], biaxial [11], triaxial [12], and blister [13–19] mechanical testers usually load samples at quasistatic or intermediate strain rates on the order of 10^−4^ ~ 10^1^ s^−1^ [20–23]. The oscillatory shear plate rheometer [10] and the recently developed needle-induced-cavitation based methods have strain rate on the order of 10^−1^ ~ 10^1^ s^−1^ [24–28]. The split-Hopkinson pressure bar (SHPB) and impact experiments can achieve high-rate mechanical testing at strain rates of 10^2^ ~ 10^4^ s^−1^ or even higher. But conventional SHPB systems using relatively hard bars are ill-suited for soft, compliant biological tissues due to a severe impedance mismatch, low signal-to-noise ratio, and difficulty in measuring large strain regimes [20, 22, 29–32].

In addition to rate limitations, many biological tissues are inherently heterogeneous and anisotropic. Many of the aforementioned techniques lack sufficient spatial resolution to resolve microscale deformation and are unable to quantitatively capture directional mechanical behavior. Atomic force microscopy (AFM) can probe tissue properties at small length scales and across a range of loading rates, but it is restricted to surface measurements. Shear-wave elastography enables subsurface characterization but suffers from limited spatial resolution and cannot readily access ultra-high strain-rate regimes [33, 34]. As a result, reliable experimental data describing anisotropic tissue behavior at extreme loading rates remain scarce.

More recently, a laser-induced cavitation (LIC)-based method, called Inertial Microcavitation Rheometry (IMR), has been developed to characterize soft materials at ultra-high strain-rates (> 10^3^ s^−1^) [35–45]. In IMR, a focused laser pulse induces optical breakdown within a soft tissue phantom and nucleates a rapid inertial cavitation bubble that loads the surrounding material at ultra-high strain rates. By analyzing the associated bubble dynamics (i.e., temporal evolution of bubble radius), we are able to quantify the hyper-viscoelastic properties of the surrounding soft material under extremely high loading rates. However, extending LIC-based methods to biological tissues presents additional challenges. Most tissues are optically opaque and strongly scattering, making it difficult to precisely deliver focused laser energy and to capture in situ bubble dynamics. To mitigate this limitation, Surya et al. [46] reduced the light scattering within a tissue by sandwiching a thinner sample (thickness is about $[eqn]$ ) between two glass coverslips, which facilitated the successful recording of the cavitation bubble dynamics. However, the thickness of the required thin samples can approach the size of the microbubbles, and the resulting confinement of the soft material layer between two boundaries changes the fundamental character of the assumed nearly infinite domain of the IMR theoretical framework. This thin-layer method will introduce the additional effects of the interfacial adhesive properties between soft materials and the glass coverslips, which further influence measured responses and complicate the interpretation of bulk material properties.

Beyond experimental limitations, the development of computational frameworks capable of modeling ultra-high-rate soft tissue dynamics and damage remains under-explored [47]. While some models account for anisotropy [48, 49] and viscoelasticity [50] exist, they are typically formulated and validated using low-to-intermediate strain-rate data. Extrapolating these models to the ultra-high-rate regime is an unvalidated assumption that yields predictions with low confidence. Therefore, there is a critical need for experimental techniques that provide microscale spatial resolution, explicitly account for anisotropy, and enable quantitative characterization of hyper-viscoelastic behavior and damage at extreme strain rates.

To address these challenges, in this work, we integrate the recently developed hydrogel fabrication manufacturing, tissue optical clearing method, with an extended laser-induced inertial microcavitation rheometry framework to be able to measure anisotropic soft viscoelastic material behaviors at ultra-high strain rates. We demonstrated this new experimental capability using both synthesized anisotropic hydrogel tissue phantoms and fresh chicken breast tissue as representative biological samples. This paper is organized as follows. We begin by introducing our experimental technique in Section 2 “Methodology”, where we utilized quasi-static uniaxial tension, shear-plate rheometer test, and laser-induced cavitation to synthesize anisotropic polyvinyl alcohol (PVA) gel and fresh anisotropic biological materials. In Section 3, we present the theoretical framework and numerical implementation, including modifications to the Keller-Miksis equation, adoption of Ogden-type constitutive models, fiber-directional bubble dynamics, and incorporation of damage mechanisms. Section 4 reports experimental and modeling results, including reconstructed stiffness characteristics, low-rate viscoelastic behavior, and analysis of radius-time ( $[eqn]$ vs. $[eqn]$ ) response during cavitation. Lastly, we conclude this paper in Section 5, “Conclusion and Discussion.”

Materials and Methods

2

Synthesis of anisotropic tissue phantoms

2.1

Biological soft tissues are complex, hierarchical composites, typically consisting of stiff collagen fibers embedded within a soft, hydrated extracellular matrix. Before directly testing fresh biological tissue, we first employed a synthetic anisotropic tissue phantom based on polyvinyl alcohol (PVA) hydrogel. This system was selected because its microstructural features and degree of anisotropy can be systematically controlled and precisely tuned. As illustrated in Fig. 1A, in our experiments, aligned PVA fibrous hydrogels were fabricated through a process combining directional freezing, solvent exchange, and mechanical training.

First, a 10 wt% PVA solution was prepared by homogeneously dissolving PVA in deionized water. The solution was then transferred into a custom-designed mold placed in direct contact with a cold finger (steel plate; Fig. 1A(i)). The vertical temperature gradient promoted ice nucleation and subsequent vertical ice growth in parallel alignment, driving PVA chains to aggregate and self-organize into honeycomb-like porous microstructures (Fig. 1B(i) and Fig. 1F(i)). Meanwhile, the abundant hydroxyl groups (-OH) along the PVA chains formed intermolecular hydrogen bonds upon concentration, and the combined effects of molecular entanglement and hydrogen bonding provided the initial physical cross-linking framework of the hydrogel [51]. After freezing, the frozen block was then immersed in a pre-cooled glycerol/H_2_O mixture (1:1, v/v) and stored at 0 °C for three days to ensure complete solvent exchange (Fig. 1A(ii)). During this process, the preformed honeycomb-like microstructures became more aggregated and stabilized, and the initially white PVA bulk underwent a macroscopic transition to a semi-transparent state.

To further enhance structural alignment, a programmable mechanical training protocol was applied following previously reported methods. Owing to the high mechanical stability of the PVA-glycerol/H_2_O gels, this process was conducted under ambient conditions. The training consisted of two sequential cyclic stretching stages. In the first stage, samples were stretched for 50 cycles at a loading rate of 1 mm/s to a maximum strain of 200% to pre-orient the polymer network. In the second stage, samples were subjected to 200 cycles at a slower loading rate of 0.1 mm/s to a maximum strain of 300%, promoting nanofibril alignment and crystalline reorganization.

We used scanning electron microscopy to confirm a marked enhancement in structural alignment after mechanical training (see Fig. 1B(ii)). The trained gels were then re-soaked in deionized water for an extended period to completely remove glycerol. The resulting hydrogels exhibited high optical transparency (see Fig. 1F(ii)) and a well-ordered hierarchical structure, spanning from microfibrous architectures observed by electron microscopy to nanocrystalline features revealed by small-angle X-ray scattering (SAXS) in the fully hydrated state (see Fig. 1B(iii) and Fig. 1E). The aligned microfibrous architecture was further verified using confocal fluorescence microscopy (see Fig. 1D).

For mechanical testing, including uniaxial tensile test, shear plate rheometry, and laser-induced inertial cavitation (LIC) experiments, the PVA hydrogel was cut into sheets with a uniform thickness of 2 mm and immersed in deionized water for at least two hours prior to each experiment to allow swelling equilibrium to be reached. This swelling process did not alter the mechanical properties of the hydrogel, as confirmed by both quasistatic uniaxial tensile tests and laser-induced inertial cavitation experiments. The average fiber width in the aligned PVA hydrogel was approximately $[eqn]$ , and the fibers were predominantly oriented along a single direction. For uniaxial tensile testing, samples were prepared as rectangular strips with a width of 5 mm and a length of 15 mm. For laser-induced cavitation experiments, the swollen samples were immersed in deionized water and affixed to the bottom of a glass-bottom Petri dish using a thick layer of grease lubricant to prevent delamination during testing. A glass coverslip was placed on top of the sample to ensure a flat and well-defined top surface for optical measurements.

Biological tissue sample preparation

2.2

Laser-induced inertial cavitation (LIC) experiments require sufficient optical transparency to enable visualization of internal bubble dynamics. Optical tissue clearing techniques generally seek to reduce the scattering between solid and liquid tissue components and phases [52]. At present, two major classes of optical clearing methods are commonly used. The first type of optical clearing method–for example, the CLARITY method [53]–removes lipids from the original tissue. In this approach, the matrix is stripped of optically turbid fat proteins and perfused by a more compliant hydrogel [53]. However, this method will change the mechanical behavior of the matrix of tissue. The second type method is to add a liquid with a closer refractive index to that of collagen and fat (refractive index: 1.43–1.53) than that of water (refractive index: 1.33), such as 50% glycerol in water and dimethylsulfoxide (DMSO) [54, 55] and more recently, tartrazine, a common FDA-approved water-soluble food coloring [56]. We recently demonstrated that using a 0.6 M tartrazine solution will enhance the optical transparency of a 2 mm-thick fresh chicken breast sample, while the material’s averaged viscoelastic material properties remain almost the same, as verified using a shear-plate rheometer and LIC using tissue phantoms (see Section 2.6).

In our experiments, unprocessed chicken breast tissue was first frozen at −20 °C and then sectioned into uniform sheets with a thickness of 2 mm using a meat slicer. The sliced tissue sheets were kept in sealed food storage bags and allowed to equilibrate at room temperature for at least one hour prior to further mechanical testing. The tartrazine clearing solution was prepared by dissolving tartrazine at a concentration of 0.6 mol per liter together with 0.2 wt% polyvinyl alcohol ( $[eqn]$ , 98% hydrolysis, Aldrich) in deionized water in a water bath of 37 °C [56]. The tartrazine dye solution was then applied to the sliced chicken breast sheet with a cotton swab. It should be noticed that the opaqueness of the tissue sample will decrease and then increase as the solution is being continuously brushed. Once opacity began to increase, the dye solution was replaced with deionized water and applied in the same manner until maximum optical transparency was achieved (see Fig. 2).

For oscillatory shear plate rheometer tests, tissue samples were punched into an 8 mm diameter circular shape. For LIC experiments, tissue samples were affixed to the bottom of a glass-bottom petri dish using a thick layer of grease lubricant to prevent delamination during cavitation. A glass coverslip was placed on top of the tissue to ensure a flat and well-defined upper surface. For uniaxial tensile testing, tartrazine dye was not applied. The 2 mm-thick tissue sheets were cut into rectangular strips approximately 5 mm wide and 15 mm long at prescribed fiber orientations. These samples were stored in sealed plastic bags at room temperature to defrost prior to testing and were exposed to ambient air during tensile experiments.

Quasistatic mechanical characterization

2.3

Quasistatic uniaxial tensile experiments were performed using a universal testing machine (FM6800; Instron, Norwood, MA, USA) to characterize the anisotropic mechanical response of synthesized polyvinyl alcohol (PVA) hydrogels and fresh chicken breast tissue samples.

To probe material anisotropy, specimens were cut such that the dominant fiber direction formed angles of 0°, 45°, and 90° relative to the loading direction. Prior to testing, samples were stored in deionized water to maintain hydration. Immediately before each test, samples were gently blotted dry using a paper towel, and a thin layer of tissue paper was wrapped around both ends of the specimen to improve grip and prevent slippage within the mechanical clamps.

Load cells with capacities of 1 kN and 100 N were used for testing PVA hydrogels and chicken breast tissue samples, respectively. All tests were conducted at a constant crosshead displacement rate of 0.2 mm/s, corresponding to a nominal strain rate on the order of 10^−2^ s^−1^. Specimens were loaded to a maximum engineering strain of 20%. To minimize the influence of grip-related artifacts, all samples were pre-stretched prior to data acquisition. For each fiber orientation, three loading-unloading cycles were performed before the fourth loading cycle, which was used for analysis in order to mitigate Mullins-type effects [58].

In addition to force and displacement data recorded by the testing machine, specimen deformation was monitored using a digital camera (FLIR BFS-U3–200S6M-C) equipped with a 75 mm lens (Kowa LM75FC24), operating at a frame rate of approximately 1/3 frames per second (fps). The recorded images were processed using a custom-developed image analysis code to track changes in specimen length and width during loading. These measurements were used to independently quantify axial strain and Poisson’s ratio, supplementing crosshead-based displacement measurements. The resulting stress-strain responses for all tested orientations are presented in Fig. 6 and discussed in Section 4.1.

Oscillatory shear plate rheometry

2.4

To characterize the dynamic mechanical properties of biological tissue at slow to intermediate loading rates and to quantify the influence of tartrazine-based optical clearing, oscillatory shear rheometry was performed using a strain-controlled shear plate rheometer (ARES-G2, TA Instruments). Both PVA and chicken breast tissue (with and without tartrazine treatment) samples were prepared to have a thickness of 2 mm, were punched into circular specimens with a diameter of 8 mm. Prior to testing, samples were stored in sealed plastic bags at room temperature to maintain hydration.

All rheological measurements were conducted using an 8 mm diameter flat circular shear plate geometry. Strain sweep tests were performed over a strain amplitude range from 0.1% to 30% at a fixed angular frequency of 1 radian per second. Frequency sweep tests were conducted over an angular frequency range from 1 to 100 radians per second at a constant strain amplitude of 1%. All tests were performed at room temperature.

Each strain sweep or frequency sweep test lasted approximately 10 minutes. Over this timescale, changes in tissue water content were assumed to be negligible. The resulting storage modulus, loss modulus, and complex viscosity data are presented and discussed in the Supplementary Materials Section S1.

Laser-induced inertial microcavitation rheometry (IMR)

2.5

We employed our developed inertial microcavitation rheometry (IMR) to measure soft materials’ ultra-high strain rate dynamic behavior with strain rate > 10^3^ s^−1^ [35, 36, 39]. In IMR, a single cavitation bubble was nucleated through a single 3–5 ns pulse from an adjustable 1–25 mJ Q-switched Nd:YAG laser (Minilite II, Continuum, Milpitas, CA) platform frequency-doubled to 532 nm [36]. The laser pulses were spatially expanded to fill the back aperture of a Nikon Plan Fluor 10×0.3 NA or a 20×0.5 NA imaging objective and were aligned through the back camera port of an inverted Nikon Ti:2E microscope (Nikon Instruments, Japan). The resulting time-dependent deformation was captured using an ultra-high-speed camera up to 10 million fps (HPV-X2, Shimadzu, Japan) for 256 frames per LIC event, with an effective spatial resolution of 3.2 or $[eqn]$ per pixel, depending on the selected objective.

Synthesized anisotropic PVA tissue phantoms and fresh chicken breast tissue samples were placed in a glass-bottom petri dish mounted on the microscope stage and covered with a glass coverslip to ensure a flat and well-defined top surface. Cavitation events were initiated at varying depths within the samples and using different laser energies in order to generate bubbles spanning a range of maximum sizes. A nondimensional depth parameter, $[eqn]$ , was defined as the ratio between the distance from the laser focal point to the bottom surface of the petri dish, $[eqn]$ , and the maximum major-axis bubble radius, $[eqn]$ . A schematic figure regarding the definition of major-axis bubble radius $[eqn]$ is provided in Fig. 7(b).

Quantifying the effects of tartrazine on soft materials’ mechanical properties

2.6

To isolate and quantify the purely mechanical effect of tartrazine on soft materials’ viscoelastic mechanical properties, we conducted oscillatory shear plate rheometer tests to verify their effects are negligible at low to intermediate strain rates (see Supplementary Materials Section S1). For ultra-high-rate mechanical behavior, we employed a cuvette geometry as in our prior work [45, 59]. Several polyacrylamide (PAAm) specimens were prepared in 4.5 mL polystyrene cuvettes by mixing 5/0.3 v/v acrylamide/bisacrylamide according to previously developed protocols [35, 59] and curing for 45 minutes. After curing, $[eqn]$ of 0.6 M tartrazine solution in water was added to the top of each cuvette (see Fig. 3(a–b)). Images were taken over the span of approximately 24 hours, and the diffusion front was monitored optically as a function of time. Concurrently, a linear 1D diffusion model was constructed to quantify tartrazine transport in PAAm numerically in MATLAB. The model assumed a well-mixed 0.6 M tartrazine solution to initially sit above the quiescent, tartrazine-free PAAm gel, after which tartrazine is permitted to flow into the gel. The gel was assumed to have a height of 31.75 mm, the bottom of which has a Neumann (zero-flux) boundary condition enforced. The 1D diffusion process was simulated using a finite-difference approach, in which concentrations were updated at each time step using a forward Euler scheme. The diffusion front images were used to calibrate the diffusivity of tartrazine in 5/0.3% PAAm for the model, which was found to be 8.5 × 10^−10^ m^2^*/*s (see Fig. 3(c)). Five identically prepared PAAm specimens with tartrazine gradients and two control PAAm specimens were included in the mechanical characterization analysis; all were tested at 1,250 minutes after curing.

Samples with tartrazine gradients were tested using LIC at 1 mm increments at locations corresponding to a range of concentrations approximated by the model. At tartrazine concentrations above approximately 0.07 M, the combination of scattering and 5 mm focusing depth into the sample led to poorly defined LIC events ^1^. Concentrations were binned into groups of 0.011 ± 0.002 M, 0.024 ± 0.005 M, and 0.048 ± 0.008 M, corresponding to approximately three LIC tests separated by 1 mm (important to ensure test independence) in each sample. The corresponding batch fits of material behavior for the three respective tartrazine groups to a Neo-Hookean Kelvin-Voigt viscoelastic material model [60] with shear modulus and viscosity $[eqn]$ were {10.2, 0.121} $[eqn]$ , {10.0, 0.131} $[eqn]$ , and {9.2, 0.121} $[eqn]$ , with individual control specimen batch fits of {8.9, 0.117} $[eqn]$ and {9.7, 0.121} $[eqn]$ .

Together, these experiments demonstrate that tartrazine-induced optical clearing does not introduce significant bias in biological tissues and soft materials’ ultra-high strain-rate mechanical properties within the tested concentration range. This validates the use of tartrazine as an effective optical clearing agent for laser-induced cavitation studies of soft materials.

Theoretical Models

3

Isotropic baseline laser-induced cavitation dynamics in soft biological materials

3.1

In this work, we extend our previously developed theoretical framework to model the laser-induced cavitation (LIC) in nonlinear viscoelastic soft biological materials [42, 60, 61]. For modeling soft materials’ nonlinear hyperelasticity, in which the material response is defined through a strain energy density function, a wide range of models has been proposed, spanning phenomenological to physically and statistically motivated approaches [62–65]. Among the most widely used are the generalized Mooney-Rivlin [66, 67], the Gent-Thomas [68], the Arruda-Boyce [69], and the generalized Ogden models [70]. By reviewing these models and considering their accuracy, computational performance, and implementation cost, we will adopt a generalized Ogden hyperelastic model [70, 71] and integrate it into the LIC modeling framework. The generalized Ogden model can describe strain stiffening, softening, and more complex nonlinear material behavior within a relatively compact formulation, and has been shown to exhibit comparable or superior performance relative to recently developed machine-learning and data-driven constitutive models [72].

To capture rate-dependent and energy-dissipating behavior of soft materials at ultra-high strain rates, we also have to move beyond simple linear viscoelastic formulations [73–75]. Classic models such as Kelvin-Voigt and Maxwell often fail to reproduce the full dynamic response observed in laser-material interactions. To overcome these limitations, we develop two nonlinear viscoelastic frameworks: a modified Poynting-Thomson [76, 77] and a Zener (or called generalized Maxwell) model [78]. As shown in Fig. 4(a), in the Poynting-Thomson model, a nonlinear hyperelastic spring A is connected in series with a Kelvin-Voigt element consisting of a second hyperelastic spring $[eqn]$ in parallel with a viscous dashpot $[eqn]$ . In the Zener (generalized Maxwell) model (Fig. 4(b)), a nonlinear hyperelastic spring $[eqn]$ acts in parallel with a Maxwell branch composed of a spring B and a dashpot $[eqn]$ in series.

The hyperelastic response of each spring is described using an Ogden strain energy density function [70]. The model parameters include shear moduli and nondimensional nonlinearity indices, while deformation is expressed in terms of principal stretch ratios in spherical coordinates. Specifically, we define spring $[eqn]$ and $[eqn]$ ’s strain energy densities as $[eqn]$ :

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ are spring $[eqn]$ and $[eqn]$ ’s shear moduli; $[eqn]$ are nondimensional parameters; $[eqn]$ , and $[eqn]$ are principal deformation stretches in a spherical coordinate system. $[eqn]$ , and $[eqn]$ represent radial, polar angular, and zximuthal angular coordinates, respectively. For small deformations, the conventional principal strain equals principal stretch minus one, i.e., infinitesimal strain $[eqn]$ , and $[eqn]$ . Here, we aim to describe large nonlinear deformations where using principal stretch ratios is more convenient. Cauchy stress tensor $[eqn]$ components $[eqn]$ can be calculated as

$$[eqn]$$

where $[eqn]$ is the hydrostatic pressure. If we assume the surrounding material is spherically symmetric and assume the surrounding material is nearly incompressible, we can present total stretch ratios $[eqn]$ , and $[eqn]$ using circumferential stretch ratio $[eqn]$ , with $[eqn]$ and $[eqn]$ are the current and undeformed distances from the bubble center:

$$[eqn]$$

Then, the Cauchy stress components of spring * ( $[eqn]$ for spring $[eqn]$ and $[eqn]$ , respectively) have the following forms:

$$[eqn]$$

$[eqn]$ and $[eqn]$ are “ $[eqn]$ ” and “ $[eqn]$ ” components of the Cauchy stress tensor $[eqn]$ . Other Cauchy stress components are zeros due to spherical symmetry.

Poynting-Thomson model formulation

3.1.1

According to the viscoelastic constitutive relationship in Fig. 4(a), we can write out the Cauchy stress and hydrostatic pressure of the dashpot as:

$$[eqn]$$

The deviatoric stress of viscous dashpot, $[eqn]$ , can be calculated as

$$[eqn]$$

We model the viscous dash-pot as a linear Newtonian whose deformation gradient tensor $[eqn]$ will obey the following evolution rule

$$[eqn]$$

where overdot means differentiation with respect to time $[eqn]$ is the viscous stretching tensor that can be linked with the mechanical properties of the dash-pot, with viscosity coefficient $[eqn]$ ,

$$[eqn]$$

The circumferential stretch rate of the dashpot will obey the following evolution rule:

$$[eqn]$$

The total deformation gradient tensor $[eqn]$ of the Fig. 4(a) model can be calculated as

$$[eqn]$$

Their circumferential stretch ratios follow

$$[eqn]$$

Combining Eqs (8–12), we can derive the temporal evolution of $[eqn]$ as

$$[eqn]$$

The momentum equation in the radial direction is

$$[eqn]$$

where $[eqn]$ is surrounding material mass density; $[eqn]$ is the radial velocity; $[eqn]$ and $[eqn]$ are “ $[eqn]$ ” and “ $[eqn]$ ” components of the deviatoric Cauchy stress tensor $[eqn]$ . By assuming the surrounding material is nearly incompressible, the radial velocity, $[eqn]$ , is

$$[eqn]$$

where $[eqn]$ is the bubble radius;. The radial acceleration can be further calculated as

$$[eqn]$$

Integrating Eqs (14–16) over the body from the bubble wall ( $[eqn]$ ) to infinite far field $[eqn]$ , and applying boundary conditions [79], we will derive the modified bubble dynamics in a nonlinear hyperelastic solid material. Taking into account acoustic radiation losses due to the compressibility of the medium, The bubble dynamics are represented by the modified Keller-Miksis equation [36, 60, 61]:

$$[eqn]$$

where $[eqn]$ is the speed of sound in the surrounding material in the far field. The pressure at the bubble wall is given by

$$[eqn]$$

where $[eqn]$ is the surface tension coefficient; $[eqn]$ the far-field driving pressure; variable $[eqn]$ the internal bubble pressure. $[eqn]$ is an integral to account for the mechanical interaction due to the surrounding viscoelastic material:

$$[eqn]$$

where $[eqn]$ is the circumferential stretch ratio at the bubble wall with $[eqn]$ the equilibrium stress free bubble radius.

Zener (generalized Maxwell) model formulation

3.1.2

In the second proposed Zener (genearlized Maxwell) model as shown in Fig. 4(b), where a nonlinear hyperelastic spring $[eqn]$ is in parallel with Maxwell units $[eqn]$ , spring $[eqn]$ has the same circumferential stretch ratio $[eqn]$ as all other branches:

$$[eqn]$$

The nonequilibrium branch $[eqn]$ shares the same deformation gradient tensor $[eqn]$ , which can be written as the multiplication of a hyperelastic deformation gradient tensor $[eqn]$ and a viscous deformation gradient tensor $[eqn]$ :

$$[eqn]$$

If we still consider dashpot $[eqn]$ ’s as linear Newtonian, their deformation evolution will follow:

$$[eqn]$$

where $[eqn]$ is the viscous stretching tensor of dashpot $[eqn]$ is dashpot $[eqn]$ ’s viscosity:

$$[eqn]$$

Tensor $[eqn]$ is the deviatoric component of the Cauchy stress of viscous dashpot $[eqn]$ whose “ $[eqn]$ ” component is denoted as $[eqn]$ . The circumferential stretch rate of $[eqn]$ is

$$[eqn]$$

For the non-equilibrium branch $[eqn]$ , due to Eq (21), its principal stretches obey the following relationship:

$$[eqn]$$

Because each spring is assumed to be nearly incompressible, we have the following constraint:

$$[eqn]$$

Spherical symmetry of deformation implies that

$$[eqn]$$ $$[eqn]$$

In this case, the stress integral $[eqn]$ has the following formula

$$[eqn]$$

Explicitly, using the Ogden material model, for spring $[eqn]$ ,

$$[eqn]$$

For the nonequilibrium branch $[eqn]$ ,

$$[eqn]$$

Plugging Eqs (30–31) into Eq (29) will provide the stress integral $[eqn]$ term, which can be further integrated into bubble dynamics (see Eq (17)).

Bubble internal pressure and boundary conditions

3.1.3

The internal bubble pressure $[eqn]$ is modeled as the sum of the partial pressures of the non-condensible gas $[eqn]$ due to optical breakdown and water vapor $[eqn]$ :

$$[eqn]$$

We assume that the partial water vapor pressure, $[eqn]$ , is equal to its saturation pressure at the bubble wall [80].

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are empirical constants. We additionally assume that the surrounding material remains isothermal everywhere with constant temperature, and the thermal boundary layer at the bubble wall is thin and negligible [80–82]. The noncondensible gas is modeled as an ideal gas, with its pressure evolving according to volume changes during bubble oscillation

$$[eqn]$$

The equilibrium non-condensible gas pressure can be obtained using the force balance at long times after cavitation,

$$[eqn]$$

where $[eqn]$ is the equilibrium bubble radius right after each LIC event (about $[eqn]$ after the laser pulse) but before it enters a much slower diffusion process with the time scale of $[eqn]$ .

Initial conditions

3.1.4

From the experimentally measured bubble radius-time curve $[eqn]$ , i.e., major radius $[eqn]$ in this paper for anisotropic PVA and chicken breast samples, we apply the initial conditions of cavitation at the moment of first time reaching $[eqn]$ . The bubble wall velocity at this time point is obtained by interpolation of experimentally measured bubble radius-time data. In our experiments, the bubble wall velocity at this time point is on the order of 10–100 m/s, which is much smaller than the wave longitudinal sound speed; therefore, the Keller-Miksis equation can be used since the bubble wall Mach number $[eqn]$ is smaller than 6.7%. Since our initial condition is applied after the beginning of the plasma optical breakdown period [83], we assume the surrounding material is still almost isothermal outside the bubble wall.

Numerical implementation

3.1.5

Our physical cavitation model is implemented using MATLAB “ode23tb” function. MATLAB’s “ode23tb” is an implementation of TR-BDF2 [84, 85], an implicit Runge-Kutta formula with a trapezoidal rule step as its first stage and a backward differentiation formula of order two as its second stage. By construction, the same iteration matrix is used in evaluating both stages. This solver provides stability and accuracy for stiff problems and has been validated in our prior cavitation numerical simulations [36, 39, 61].

Modeling surrounding material damage

3.2

Material damage model

3.2.1

When cavitation-induced hoop stress exceeds a threshold, the surrounding material fails. Previous studies by Glinsky et al. [86] and Luo et al. [37] employed simplified yield-stress models to describe cavitation damage in hydrogels. Movahed et al. [87] and Milner & Hutchens [88, 89] developed an energy-based threshold derived from Griffith’s fracture theory such that fracture initiates and propagates when the elastic strain energy released by crack growth exceeds the intrinsic surface energy required to form new fracture surfaces.

However, for LIC experiments, it is hard to quantify the newly created surfaces in the surrounding materials [90–92]. Here, we introduce a new material failure criterion that, when the local circumferential stretch exceeds the ultimate extensibility of the polymer chains, which is denoted by a threshold $[eqn]$ , the material locally loses its load-bearing capacity. Mathematically, this behavior is captured by scaling the strain energy density function with a scalar damage variable $[eqn]$ that ranges from zero (fully damaged) to one (undamaged). The damage variable $[eqn]$ depends on both spatial position and deformation history and records the maximum level of degradation experienced during tensile loading. In this work, we model spring $[eqn]$ as a damageable spring in both Poynting-Thomson and generalized Maxwell models (see Fig. 5):

$$[eqn]$$

where $[eqn]$ is a function with respect to $[eqn]$ , and is defined as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are radial coordinates of material points in deformed and reference configurations, respectively. The damage coefficient $[eqn]$ is treated as a history-dependent variable that records the maximum level of degradation experienced by the polymer network during the tensile expansion phase. If we assume that, during bubble expansion, the local material state transitions abruptly from fully intact to complete failure, the response corresponds to an ideally plastic behavior when $[eqn]$ . This assumption is consistent with the simplified rupture models reported in Refs. [37, 86, 93]. However, different from previous studies in Refs. [37, 86, 93] that assumed damage initiation during bubble expansion, our experimental observations indicate that significant damage manifests primarily during the collapse phase and subsequent oscillations.

To account for asymmetric behavior between expansion and compression, two damage models are considered. In “Model 1”, damage is symmetric with respect to radial expansion and compression. In “Model 2”, damage accumulates only during radially expanding loading, while materials become undamaged, i.e., $[eqn]$ when the circumferential stretch ratio is smaller than one.

Numerical implementation & parameter calibration

3.2.2

Because the damage variable evolves both spatially and temporally, we discretized the reference radial distance from the bubble wall outward to a sufficiently large cutoff radial distance, i.e., we used $[eqn]$ in our numerical implementations. At each time point and each reference radial distance, we calculated its local circumferential stretch ratio to determine the value of $[eqn]$ using the damage criterion in Eq (37).

In the parameter calibration procedure, we follow our previous fitting process [39] by optimizing the least squares (LSQ) error as the $[eqn]$ norm of the difference between the experimental data points and the numerical simulations for $[eqn]$ . A global optimization scheme, “pattern search” [94], was applied to search for a global minimum in the LSQ error, and the relative optimization tolerance parameter was set to be 0.01. Since damage is assumed to be obvious during collapse, elastic and viscous parameters are first calibrated using the expansion phase, followed by identification of the critical stretch during collapse.

Anisotropy: experimental manifestations and modeling limitations

3.3

Both anisotropic PVA and chicken breast materials are anisotropic that can be modeled as unidirectional fiber-reinforced composites, where the overall response is a combination of the ground matrix behavior and contributions from embedded fiber families. This modeling framework has successfully captured the nonlinear anisotropic behavior of tissues such as brain matter [95], skin [96–99], and arterial walls [100, 101].

Experimentally, anisotropy manifests as nonspherical cavitation bubble shapes, characterized by distinct major and minor axes (see Fig. 7(a–b) and Fig. 8, where the major and minor axes are labeled as $[eqn]$ and $[eqn]$ , respectively). Accurately modeling such behavior would require fully three-dimensional finite element simulations incorporating anisotropic hyper-viscoelastic constitutive laws. Because the present work focuses primarily on experimental characterization, we restrict our analysis to the major-axis bubble dynamics and assume that deformation along the minor axis does not significantly influence the dominant response. A fully coupled three-dimensional anisotropic modeling framework will be the subject of future work.

Results and Discussion

4

Quasistatic material characterization

4.1

The tested specimens exhibit a single dominant fiber direction; therefore, we model them as transversely isotropic materials. In the small-strain regime, the PVA hydrogel response is well approximated by linear elasticity and can be described by the generalized Hooke’s law,

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the Cauchy stress and infinitesimal strain vectors in Voigt notation, and $[eqn]$ is the fourth-order stiffness tensor represented as a symmetric 6 × 6 matrix. A general linear elastic solid has 21 independent elastic constants. Under transverse isotropy (with the symmetry axis aligned with $[eqn]$ ), the number of independent constants reduces to five, yielding

$$[eqn]$$

The shear modulus in the plane of isotropy ( $[eqn]$ plane) is related to the in-plane normal moduli

$$[eqn]$$

To quantify anisotropy experimentally, we performed uniaxial tensile tests with loading directions oriented at 0°, 45°, and 90° relative to the dominant fiber direction. In addition to force and crosshead displacement measurements, optical image post-processing was used to quantify transverse strain and extract Poisson’s ratios. Young’s modulus was fitted in the small-strain regime (0.05%−0.5%). The directional Young’s moduli are denoted by $[eqn]$ , and $[eqn]$ . The in-plane Poisson’s ratios measured from the 0°, 45°, and 90° tests are denoted by $[eqn]$ , and $[eqn]$ , respectively.

We first construct the compliance matrix $[eqn]$ using the measured Young’s moduli and Poisson’s ratios:

$$[eqn]$$

with components given by

$$[eqn]$$

Poisson’s ratio $[eqn]$ does not appear explicitly in Eq (42) because $[eqn]$ and $[eqn]$ have the following relation:

$$[eqn]$$

The transverse Poisson’s ratio $[eqn]$ was not directly measured in experiments. Because the PVA hydrogel is nearly incompressible under loading in the $[eqn]$ plane, we approximate $[eqn]$ and use $[eqn]$ in the analysis.

Figure 6 summarizes the measured cyclic stress–strain responses for PVA hydrogels (Fig. 6a–c) and chicken breast tissue (Fig. 6d–f). Each loading orientation (0°, 45°, and 90° relative to the fiber direction) was tested three times. For PVA hydrogels, the mean Young’s moduli are $[eqn]$ , and $[eqn]$ . Thus, the modulus measured along 0° is approximately three times larger than that measured along 90°, confirming strong anisotropy. The PVA response is approximately linear elastic in this regime, and loading–unloading hysteresis is negligible. The resulting averaged stiffness matrix $[eqn]$ of anisotropic PVA hydrogel is:

$$[eqn]$$

Compared to PVA, chicken breast tissues present larger hysteresis [58] where the unloading stress-strain curve does not overlap with the loading curve. For small linear loading regimes, the mean Young’s moduli of chicken breast tissues are $[eqn]$ , and $[eqn]$ . Based on our experimental observations, we measured the apparent Poisson’s ratios $[eqn]$ and $[eqn]$ from the captured image sequences using Fiji and found that they are close to zero in quasistatic tensile tests. We also assume that the chicken breast is nearly incompressible under loading in the $[eqn]$ plane; thus, we approximate $[eqn]$ and use $[eqn]$ in the analysis. The resulting averaged stiffness matrix $[eqn]$ of fresh chicken breast tissue is:

$$[eqn]$$

Laser-induced inertial cavitation results

4.2

For each recorded LIC event, we performed image post-processing to extract the bubble contour by thresholding the grayscale intensity. Figure 7(a) presents a series of selected images within a recorded LIC event in a PVA hydrogel. The extracted bubble contour at each time point was fitted with an ellipse, from which the two principal radii were obtained and denoted as the major axis $[eqn]$ and minor axis $[eqn]$ (see Fig. 7(b) for definitions). The fitted $[eqn]$ vs $[eqn]$ curves were displayed in Fig. 7(c).

To facilitate comparison across experiments, each $[eqn]$ curve was normalized by its maximum value $[eqn]$ . The nondimensional radii are defined as $[eqn]$ and $[eqn]$ . Nondimensional time is defined as $[eqn]$ , where $[eqn]$ is a characteristic speed, $[eqn]$ is the far-field atmospheric pressure, and $[eqn]$ is the material mass density. The normalized $[eqn]$ vs $[eqn]$ curves are shown in Fig. 7(d). The nondimensional ratio $[eqn]$ vs $[eqn]$ curves are shown in Fig. 7(e). The near overlap of all normalized curves in Fig. 7(d,e) indicates that the bubble dynamics are governed by the surrounding viscoelastic material response and are largely insensitive to bubble size over the range of $[eqn]$ .

We applied the same image processing and normalization techniques to the LIC events recorded in the fresh chicken breast tissue. Figure 8(a) presents a series of selected images within a recorded LIC event in a chicken breast tissue sample. Figures 8(d,e,f) display the fitted major radius $[eqn]$ vs time $[eqn]$ curves, nondimensional $[eqn]$ vs $[eqn]$ and $[eqn]$ vs $[eqn]$ curves, respectively. The time-lapse sequence in Fig. 8(a) and the radius-time plot in Fig. 8(d) reveal that the bubble undergoes a single expansion phase followed by a slow, heavily damped collapse. Unlike the PVA hydrogel (Fig. 7(c)), there are no strong second and third bubble rebounds. This pronounced asymmetry between expansion and collapse, with a significantly slower collapse phase, is direct evidence of substantial and irreversible energy dissipation. While the PVA hydrogel responds in a primarily viscoelastic manner (storing and returning elastic energy to drive rebounds), the chicken tissue’s response is more affected by inelastic processes.

The structural anisotropy of both the PVA tissue phantom and chicken breast tissue strongly influences bubble morphology. As shown in Fig. 7(b) and Fig. 8(b), the bubble expands non-spherically, with the major axis $[eqn]$ consistently aligned with the local fiber direction. This anisotropy is further quantified by the aspect-ratio evolution plots (Fig. 7(e) and Fig. 8(f)). For the PVA hydrogel, $[eqn]$ remains above 0.9 until the time of maximum expansion, indicating an almost spherical initial expansion. The ratio decreases during each collapse and recovers during subsequent rebounds, with its value strongly correlated with $[eqn]$ . After approximately four to five rebound cycles, $[eqn]$ approaches a steady value of about 0.75 as the bubble relaxes toward equilibrium.

In contrast, the chicken breast tissue exhibits markedly different behavior. The $[eqn]$ ratio decreases during the initial expansion and reaches a local minimum of approximately 0.6–0.8 at the time of maximum $[eqn]$ . The ratio then increases briefly during the first collapse but subsequently decreases monotonically in the absence of significant rebounds, reaching values of 0.4–0.6 as the bubble approaches its equilibrium size. These contrasting behaviors can be attributed to fundamental differences in material response. The PVA hydrogel exhibits nearly linear elastic behavior with negligible Mullins effects and minimal hysteresis under tensile loading, enabling strong elastic energy storage and multiple bubble rebounds. In contrast, chicken breast tissue displays pronounced hysteresis and damping, which suppress rebound dynamics. Additionally, LIC experiments in chicken breast reveal a characteristic lemon-shaped bubble morphology associated with fiber-aligned delamination (Fig. 8(c)). This observation suggests that the tissue fracture toughness is lower parallel to the muscle fibers than perpendicular to them, promoting inter-fiber delamination. The associated damage provides an additional dissipation pathway and explains both the rapid decay of bubble oscillations and the pronounced reduction in $[eqn]$ at the maximum expansion stage.

Directional moduli along fiber directions ultra-high rate characterization

4.3

Poynting-Thomson model results

4.3.1

Figure 9 compares the experimentally measured major-axis bubble dynamics with predictions from four different constitutive descriptions based on the Poynting-Thomson framework: a purely elastic model without damage (no-damage model; $[eqn]$ ), a fully damaged model $[eqn]$ , and two progressive damage models (Model 1 and Model 2). These comparisons provide direct insight into the effective directional mechanical response along the dominant fiber direction under ultra-high-rate loading.

For the anisotropic PVA hydrogel major radius vs time curves (Fig. 9a), the no-damage model over-predicts rebound amplitudes and bubble persistence after the first collapse, indicating excessive elastic energy storage. In contrast, the fully damaged model underestimates resistance during the expansion phase and fails to capture the observed rebound timing. Model 1 can capture the bubble radius versus time, except for near-bubble collapse time points. Among the tested formulations, Model 2, which allows damage to accumulate preferentially during tensile expansion but preserves compressive stiffness during collapse, best reproduces both the expansion peak and subsequent rebound amplitudes. This behavior is consistent with the quasi-elastic, low-hysteresis response observed in the quasistatic tensile tests, where PVA exhibits minimal Mullins effects and negligible permanent damage.

The fitted parameters in Table 1 further support this interpretation. The PVA hydrogel exhibits a large effective shear modulus in spring A $[eqn]$ and a relatively high critical circumferential stretch ratio $[eqn]$ , indicating that the polymer network can sustain substantial elastic stretching before losing load-bearing capacity. The resulting damage zone remains localized near the bubble wall, consistent with the limited irreversible deformation observed experimentally.

The fresh chicken breast tissue (Fig. 9b) displays different behavior. The no-damage model fails to reproduce the strong damping and absence of multiple rebounds, while the fully damaged model collapses too rapidly and underpredicts peak bubble size. Once again, Model 2 provides the closest agreement, accurately capturing the asymmetric expansion-collapse response and the rapid decay of oscillations. This result indicates that tissue damage initiates primarily during tensile loading along the fiber direction and persists through subsequent collapse phases, leading to irreversible energy dissipation.

Quantitatively, the fitted directional modulus for chicken tissue at ultra-high strain rates is over two orders of magnitude smaller than that of PVA $[eqn]$ , reflecting the relatively compliant nature of chicken tissue compared to PVA. The critical stretch for damage initiation, $[eqn]$ , is much smaller than the PVA and is close to one, $[eqn]$ , implying that chicken breast tissue exhibits early-onset damage and pronounced anisotropic softening along fiber directions. The inferred damage-zone radius is significantly larger than in PVA, consistent with the experimentally observed fiber-aligned delamination and “lemon-shaped” bubble morphology shown in Fig. 8(b).

Generalized Maxwell model results

4.3.2

Figure 10 presents fitting results for the same experimental data using the generalized Maxwell (Zener-type) viscoelastic framework, again considering four constitutive descriptions: no damage $[eqn]$ , fully damaged $[eqn]$ , and two progressive damage formulations (Model 1 and Model 2). Compared to the Poynting–Thomson model, the generalized Maxwell formulation places greater emphasis on stress relaxation through the non-equilibrium Maxwell branch, providing an alternative representation of rate-dependent dissipation at ultra-high strain rates.

For anisotropic PVA hydrogels (see Fig. 10a), the generalized Maxwell model captures the initial expansion peak and early rebound behavior reasonably well across all damage assumptions. However, similar to the Poynting–Thomson case, the no-damage model overpredicts rebound amplitudes, while the fully damaged model underestimates resistance during expansion. Model 1 improves agreement but still exhibits phase errors near collapse. Model 2 again yields the best overall agreement, reproducing both the amplitude and timing of successive rebounds. These results indicate that, even when viscoelastic dissipation is represented through a Maxwell-type relaxation mechanism, damage evolution during tensile expansion remains essential for accurately capturing bubble dynamics in PVA.

For chicken breast tissues (see Fig. 10b), the generalized Maxwell model reproduces the strong damping and absence of higher-order rebounds more naturally than the Poynting–Thomson model, owing to the enhanced stress relaxation inherent in the Maxwell branch. Nevertheless, the no-damage model still fails to capture the rapid decay of oscillations, while the fully damaged model collapses too quickly. Consistent with previous observations, Model 2 provides the closest match to experimental data, accurately reproducing the asymmetric expansion–collapse response and the monotonic relaxation toward equilibrium.

Overall, the generalized Maxwell results reinforce the conclusion that direction-dependent damage, rather than viscoelastic dissipation alone, governs cavitation-induced energy loss in anisotropic biological tissue. The Maxwell framework improves the representation of rate-dependent relaxation for both PVA and chicken breast tissue. These findings further highlight the necessity of coupling anisotropic damage evolution with high-rate viscoelastic models when interpreting laser-induced cavitation in soft biological materials.

Comparison between Poynting-Thomson and generalized Maxwell models

4.3.3

Figure 11 and Table 1 compare the best-performing damage formulation (Model 2) within the Poynting-Thomson and generalized Maxwell frameworks against experimentally measured major-axis bubble dynamics for both PVA hydrogel and chicken breast tissue. For both materials, the gray-shaded regions indicate the portions of the radius-time history used for parameter identification.^2^ Overall, both models successfully reproduce the dominant features of the experimental response within the calibration window. Nevertheless, systematic differences emerge in their treatment of rate-dependent dissipation, rebound dynamics, and inferred material parameters. Overall, both models successfully reproduce the dominant features of the experimental response within the calibration window; however, systematic differences emerge in their representation of rate-dependent dissipation, rebound dynamics, and inferred material parameters.

For the PVA hydrogel, both rheological frameworks capture the initial expansion, collapse, and subsequent rebound sequence with good fidelity (Fig. 11, left column). The fitted equilibrium and maximum bubble radii are identical between the two models (Table 1), confirming that the inferred geometric and kinematic descriptors are robust and model independent. Differences arise primarily in the partitioning of elastic and viscoelastic contributions. The Poynting-Thomson model attributes the high-rate response to a stiff primary elastic spring $[eqn]$ with a relatively small viscous contribution, whereas the generalized Maxwell model assigns a significantly lower elastic stiffness $[eqn]$ and relies more strongly on stress relaxation through the Maxwell branch. Notably, the stiffness calibrated using the generalized Maxwell model is of the same order as that obtained from low-loading-rate tensile tests. Despite these differences, both models converge on a similar critical circumferential stretch for damage initiation ( $[eqn]$ ; see Table 1 critical circumferential stretch “ $[eqn]$ ”) and an identical damage-zone radius, indicating that the onset and spatial extent of damage are largely insensitive to the specific viscoelastic architecture for this nearly elastic material.

In contrast, more pronounced differences between the two models are observed for chicken breast tissue (Fig. 11, right column). While both frameworks reproduce the strong damping and absence of higher-order rebounds when damage is included, the generalized Maxwell model provides a more natural representation of the smooth, monotonic relaxation toward equilibrium, even outside the fitting window. Quantitatively, the generalized Maxwell model yields a larger viscosity compared to the Poynting-Thomson model ( $[eqn]$ versus 0.23 Pa·s in the Poynting-Thomson model) and a higher stiffness contribution from the secondary elastic branch ( $[eqn]$ versus 0.30 MPa), reflecting the greater importance of stress relaxation mechanisms in biological tissue. Nevertheless, both models independently identify a critical stretch close to unity $[eqn]$ , consistent with early-onset damage and pronounced anisotropic softening along muscle fiber directions.

Despite their differing rheological constructions, both models predict comparable damage-zone radii for chicken breast tissue, on the order of several hundred micrometers (see Table 1), consistent with experimentally observed fiber-aligned delamination. This agreement underscores that damage initiation and spatial extent are governed primarily by stretch-based failure criteria rather than the specific form of viscoelastic dissipation.

Since the generalized Maxwell model provides better fitting accuracy than the Poynting–Thomson model, we fitted four samples for each material type and summarized the extracted material properties in Table 2. Taken together, these results indicate that although the Poynting–Thomson and generalized Maxwell models distribute elastic and viscous contributions differently, they converge on consistent damage thresholds and damage-zone sizes for both synthetic and biological materials. For relatively elastic systems such as PVA hydrogels, either framework provides comparable predictive capability; however, the generalized Maxwell model yields fitted spring stiffness values on the same order as those measured at lower strain rates, whereas the Poynting–Thomson model produces fitted spring stiffness values that are approximately two orders of magnitude higher. For highly dissipative biological tissues, the generalized Maxwell model offers an improved representation of relaxation behavior. For both models, accurate interpretation of laser-induced cavitation dynamics ultimately requires explicit incorporation of direction-dependent damage and, more fundamentally, a fully three-dimensional finite-element formulation. These extensions will be pursued in future work.

Conclusion

5

In this work, we investigated the rate-dependent and anisotropic viscoelastic behavior of polyvinyl alcohol (PVA) hydrogels and fresh chicken breast tissue by integrating quasi-static mechanical testing, oscillatory shear rheometry, and ultra-high-strain-rate inertial microcavitation rheometry. The combination of these complementary techniques allowed us to quantify soft materials and biological tissue mechanics over more than five orders of magnitude in strain rate, and to directly link their ultra-high-rate deformation to underlying damage mechanisms.

First, quasi-static uniaxial tensile tests performed along 0°, 45°, and 90° fiber orientations enabled the reconstruction of the full stiffness matrices for PVA hydrogel under small-strain conditions. PVA hydrogels exhibited strong transverse isotropy with orientation-dependent stiffness (with $[eqn]$ ). Cyclic loading along a single nominal orientation (0°) exhibits that the PVA hydrogel has minimal hysteresis and negligible Mullins effect, consistent with its nearly ideal elastic behavior in this regime. In contrast, cyclic loading experiments revealed a markedly different mechanical response of chicken breast tissue: the loading and unloading curves showed substantial hysteresis and pronounced Mullins softening even under small strains. These observations demonstrate that, unlike PVA, the native muscle tissue dissipates significantly more energy during repeated loading.

Second, oscillatory shear rheometry showed that applying tartrazine-based optical clearing, which is a critical step for enabling IMR in opaque biological tissues, induces less than a 5% difference in storage modulus, loss modulus, and complex viscosity across frequencies from 1–100 rad/s under a strain of 1% (see Supplementary Material Section S1). Using polyacrylamide (PAAm) tissue-mimicking phantoms and IMR, we further validated that the effect of tartrazine dye on ultra-high-strain-rate viscoelastic properties is also within 10%. These results confirm that the clearing protocol enhances optical transparency while largely preserving the intrinsic rheological behavior of the material.

Third, we employed laser-induced cavitation (LIC) experiments to probe the ultra-high-rate behavior of both PVA hydrogels and chicken breast tissue. Time-resolved imaging revealed preferential bubble elongation along the fiber direction, reflecting the anisotropy of the surrounding material. In both materials, the measured bubble radius-time curves exhibited pronounced asymmetry between the first expansion and subsequent collapse phases. In chicken breast tissue, we additionally observed fiber-parallel delamination once the circumferential stretch exceeded a critical threshold. These observations indicate that irreversible, inelastic dissipation plays a central role in governing ultra-high-strain-rate mechanical response in biological tissue.

Fourth, we developed and implemented two hyper-viscoelastic models, based on the Poynting-Thomson and generalized Maxwell frameworks, to interpret the measured major-axis radius-time curves. We further introduced a fiber-stretch-based material damage model. Analysis of the bubble temporally evolving morphology enabled the extraction of effective high-rate viscoelastic properties and critical stretches associated with the onset of material failure. It is important to note, however, that the present numerical framework fits only the major-radius trajectory, $[eqn]$ , and does not yet explicitly incorporate anisotropic constitutive behavior. As a result, while the model accurately reproduces the temporal dynamics of major-axis bubble expansion and collapse, it does not attempt to predict the observed bubble-shape anisotropy or preferential delamination along muscle fibers. A fully three-dimensional finite element formulation will be pursued in our future work.

In summary, this study introduces a new experimental methodology for characterizing anisotropic soft materials at ultra-high strain rates. The combination of optical clearing, IMR, low-rate to ultra-high-strain-rate mechanical characterization and modeling frameworks advances the current understanding and quantification capabilities of rate-dependent soft-tissue biomechanics. The experimental data and analysis approaches obtained from this introduction of experimental and computational methods will benefit future studies on high-strain-rate material damage and tissue injury, laser and ultrasound-related medical procedures, and anisotropic tissue constitutive modeling.

Supplementary Material

Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

LICchickenbreastsupp20260210.pdf