Elastic cavitation and fracture via injection

Abstract

The cavitation rheology technique extracts soft materials mechanical properties through pressure-monitored fluid injection. Properties are calculated from the system's response at a critical pressure that is governed by either elasticity or fracture (or both); however previous elementary analysis has not been capable of accurately determining which mechanism is dominant. We combine analyses of both mechanisms in order to determine how the full system thermodynamics, including far-field compliance, dictate whether a bubble in an elastomeric solid will grow through either reversible or irreversible deformations. Applying these analyses to experimental data, we demonstrate the sensitivity of cavitation rheology to microstructural variation via a co-dependence between modulus and fracture energy.

---

Rapid, 3-D mechanical characterization of soft solids provides information necessary for the elucidation of health-related mechanisms, e.g., mechanotransduction in cells and tissue property evolution during onset of disease. A recently explored technique referred to as cavitation rheology (CR)^1–10 or the maximum bubble pressure (MBP)¹¹ method provides a promising avenue for quickly and inexpensively obtaining this data, especially in soft-solids whose properties can be otherwise difficult to measure using traditional techniques. Past CR measurements in a variety of synthetic, polymeric materials at a range of length scales have been found to correlate with elastic modulus.^1,2,4,7,11 These measurements arise from an instability governed mechanism that causes the sudden, rapid, and temporarily unstable expansion of a pressurized void at the tip of an embedded needle (Fig. 1). This expansion results in a pressure drop within the system and thus a characteristic peak pressure, P_c, which can, notably, be obtained without monitoring deformation visually.

Fig. 1 Characteristic cavitation data in a PMMA–PnBA–PMMA triblock co-polymer (7 wt%, 20 °C). Images show the bubble just prior to and after the cavitation event at critical pressure P_c (collected at 40 fps). Scale bars are 200 μm.

To date, only elementary analyses have been applied to these observations, with a limited quantitative understanding as to when and how the observed deformation mechanisms relate to the measured data.² Although the morphology of the expanded voids clearly indicate that fracture can be a dominant growth mechanism for certain materials, application of the simple ‘cavitation equation’² based only on elastic deformations can typically yield modulus values having less than 15% error. In contrast, correlating P_c to the fracture toughness, G_c, results in much greater error. This result points to the need for further understanding into how materials properties, including nonlinear behavior, control unstable void growth in an elastomeric solid. This insight is critical for enhancing the implementation of CR in biological tissues, where nonlinear mechanics play a dominant role. More generally, resolving apparent conflict between the observed deformation path and the critical pressure is important for understanding how damage is incurred in soft materials, a topic that has been studied extensively, most notably by Alan Gent, though many questions remain.

In their 1959 paper, Gent and Lindley¹² observed a discontinuity in the load-displacement response of rubber disks under uniaxial tension. Noting that the predominant stress state in this testing geometry is hydrostatic, they attributed the discontinuity to the rapid expansion of pre-existing spherical defects or voids within the sample. These observations ignited a discussion of cavitation (and fracture) within an elastomeric solid (e.g.ref. 13–15) that continues even now.^16,17 Our extension of this discussion in the context of CR provides new understanding of how elasticity, fracture toughness, and the full system thermodynamics, control void growth mechanics in soft solids.

In the following work, we utilize a simplified, spherical geometry under quasi-static loading in order to quantitatively determine the interplay between elastic and fracture mechanisms and the material parameters and behaviors that affect that interplay, in particular strain hardening. Knowledge gained from understanding these mechanisms individually is then unified into a mechanism map of anticipated behavior from which predictions for P_c are made for three distinct experimental systems. Specifically, we follow previous derivations for a spherical void in an incompressible neo-Hookean solid.^18–20 We incorporate energy-based fracture criteria for the void in order to determine the bubble deformation at which fracture occurs. From these models, we construct mechanism maps for two limiting-case loading scenarios. The resulting maps predict the critical pressure P_c as a function of the non-dimensionalized elastocapillary length γ/μA and relevant geometric and material parameters. The elastocapillary length γ/μ is the ratio of the surface tension, γ, to the elastic shear modulus, μ. A is the initial radius of the void, which we assume is equivalent to the needle radius used in CR. The relationship between these length scales is thus essential for interpreting experimental P_c data as an intrinsic material property.

1 Cavitation rheology: boundary and loading conditions

The details of CR experiments have been given previously.^1–9,11,21 Essentially, a needle is embedded within a soft solid. This is followed by pressurization of the fluid within the needle that leads to the growth and eventual rapid, unstable expansion of a bubble at the tip of the needle (Fig. 1). For the analyses developed here, we use expressions for an incompressible neo-Hookean solid in consideration of space, but leave the derivation open to the incorporation of any isotropic solid defined by a free energy function. In addition to an incompressible, neo-Hookean solid, we show results for an incompressible, strain hardening solid described by a modified Yeoh model (Section 7 and Fig. S2, ESI†). We consider CR loading conditions that arise from two limiting cases for the choice of pressurizing fluid.

In experiments, the pressurizing fluid has been either air or water^1–8,11 and in both cases unstable expansion at a critical pressure P_c is observed. Because the fluid reservoir has an inherent stiffness, k, associated with it (due to rubber o-rings and pressure sensors), we approximate the system as shown in Fig. 2c and d. This approximation leads to two limiting cases for loading conditions based on simplifications when the fluid is either a gas (ideal gas) or a liquid (incompressible fluid). In the case of the former, Zhu et al.²⁰ previously showed that sudden elastic cavity growth (elastic cavitation) requires both a sufficiently large bubble surface energy, γ, and an initial fluid reservoir volume, V_c0, capable of achieving a pressure that can overcoming the surface energy contribution the bubble. Their system was comprised of a reservoir of ideal gas molecules and assumed no system compliance. Here we use this same framework, setting system stiffness equal to infinity (compliance equal to zero) in the case of a gas pressurizing fluid. This approximation is valid even for a system with finite stiffness due to the fact that the ratio of the system stiffness, k = dP/dV|_spring, (k ≈ 1 kPa μL⁻¹, Fig. S7, ESI†) to the ‘stiffness’ of an ideal gas, dP/dV|_gas, at the onset of cavitation is much greater than one in all cases except that of a high modulus material (∼1 MPa) combined with a small reservoir (<1 mL). Materials with such high moduli are readily characterized using traditional techniques and therefore outside the main scope of CR. This stiffness ratio is determined for a chamber volume V_c (initially having volume V_c0 at atmospheric pressure P₀) evaluated at a critical pressure P_c that scales roughly with the material shear modulus so that kV_c0P₀/P_c² ≈ kV_c0P₀/μ² ≫ 1. Thus, for air as the pressurizing fluid, we will simply follow Zhu et al.²⁰ derivation. However, we include a small redefinition of the initial state of the system that allows for experimental comparison as will be discussed in the next section.

Fig. 2 Schematic drawing of a bubble within an infinite elastic solid illustrating the process by which bubble deformation, Λ, is defined via (a) sphere removal, (b) addition of surface tension, and (c) bubble expansion through pressurization. (d) Illustration of the functional components of the pressurizing syringe including the initial volume, V_c0, and system stiffness, k.

As shown by Zhu et al.,²⁰ the gas provides the only negative contribution to the free energy of the system just described (i.e., a system having an infinite stiffness spring) and thus enables its ability to cavitate. Application of this same model to an incompressible fluid⁴ therefore, cannot predict elastic cavitation. Because addressing both gas and system compliance contributions simultaneously becomes unwieldy, our second limiting case is simply that of a spring of finite stiffness only (no gas reservoir).¹⁹ We demonstrate that this system compliance provides the necessary negative free energy contribution to allow for elastic cavitation with an incompressible pressurizing fluid. In summary, our two limiting cases for loading conditions become: a reservoir of ideal gas molecules (no spring) and a linear spring of stiffness k = dP/dV. Note that in neither case are these load or displacement-controlled conditions. (Displacement-controlled, or more accurately volume-controlled, testing applies only to the case of a spring of infinite stiffness combined with an incompressible fluid.) Similar to a cantilever-type, ‘displacement-controlled’ compression test, the compliance of either the gas or the system (like the centilever) mediates any deformation of the test sample. In CR, that effect is an essential part of the observed, catastrophic response that leads to P_c.

2 Elastic model

We assume a spherical void geometry within an infinite, incompressible hyperelastic solid. (Previous calculations found that the needle-tip geometry plays a minimal role in the elastic pressure–stretch response of the void.²¹) The void is connected to a chamber representing the syringe via a hypothetical fluid transporter²⁰ (Fig. 2c and d); a linear spring accounts for loading system compliance (Fig. 2d and Fig. S7, ESI†). As mentioned above, this geometry and loading (without the spring) has been considered previously for the elastic deformation regime; however, the initial configuration and pressure condition used was not compatible with CR, making quantitative comparisons cumbersome.²⁰ To address this issue, we define the bubble deformation, Λ, as the bubble radius a normalized by the initial void radius A,where A₀ is the radius of a void in an unstrained elastomer in the absence of surface tension (Fig. 2a) such that Λ₀ = a/A₀ is the bubble stretch ratio that corresponds directly to the hoop stretch in the elastomer.²⁰ The initial, equilibrium deformation, Λ_eq, that sets A is determined for A₀ subject to surface tension at the void–solid interface (Fig. 2b). One can obtain Λ_eq, as Gent and Tompkins²² did by solving for the stretch that satisfies thermodynamic stability criteria for a total free energy, F_tot, comprised of only the surface energy and elastic strain energy, i.e., ∂²F_tot/∂Λ₀² = 0 (or equivalently ∂P/∂Λ₀ = 0) (see ESI†). Defining the bubble deformation, Λ, in such a way effectively sets the initial pressure within the void/chamber to atmospheric, P₀, (Fig. 2) consistent with an experimental CR setup which starts at atmospheric conditions.

The pressure–stretch (P vs. Λ₀) and free energy–stretch (F_totvs. Λ₀) responses²⁰ previously derived without this correction are then simply modified through the substitution of Λ₀ (eqn (1)) such that they become functions of Λ given that Λ_eq = Λ_eq(γ/μA, constitutive model parameter). For an incompressible neo-Hookean solid with strain energy function W = μ(I − 3)/2 (where I = 2Λ₀² + Λ₀⁻⁴: the first strain invariant of the right Cauchy–Green tensor) and assuming an ideal gas pressurizing fluid, these expressions are:

Subscripts ‘surf’, ‘elas’, ‘gas’, and ‘env’ correspond to contributions from the surface/interface of the void, elastic strain energy in the elastomer, the gas within the void/syringe, and the work done on the environment by the expanding, infinitely large elastomeric sphere. F_gas and F_env are determined from ideal gas and P–V work expressions, respectively. CF corresponds to the compression fraction 1 − V_c/V_c0 where V_c is the chamber volume, which decreases throughout an experiment as illustrated in Fig. 3c.

Fig. 3 (a) Pressure–deformation response for a neo-Hookean void having surface energy within an infinite elastic solid with A = 50 μm, μ = 250 Pa, and γ = 27.7 mN m⁻¹ (black dashed line) and for which elastic cavitation is impossible, γ/μA = 1.4953 (gray dashed line). Overlays of decreasing saturation correspond to decreasing gas chamber size V_c (i.e., depressing the syringe) for compression values that result in single (pre-cavitation), absolute stability, metastability, and single (post-cavitation) equilibria. (V_c0 = 1 mL, P₀ = 1 atm). (b) Free energy–stretch response for the combined void/syringe system plotted in (a) illustrating absolute and meta-stability. Saturation values correspond to those plotted in (a). (ΔF_tot = F_tot(Λ) − F_tot(1).) (c) Schematic showing chamber volume V_c in relation to the initial chamber volume V_c0 which defines the compression fraction, 1 − V_c/V_c0. In both (a and b), vertical lines exemplify three scenarios in which the deformation at the onset of fracture, Λ_F might relate to the deformation at the onset of fracture, Λ_c, and/or the deformation after elastic cavitation has occurred, Λ_cf. The latter locations are marked with gray arrows.

For the incompressible fluid approximation accounting for system stiffness, k, the expression for F_gas is replaced with one for the spring, F_spr.

The squared difference is the total displaced volume, i.e., the total applied volumetric displacement V_app = V_c0 − V_c less the increased volume of the bubble. k is the spring compliance, dP/dV, assumed constant here (Fig. S7, ESI†). By varying V_app, a set of curves similar to those for varying CF in Fig. 3b are obtained (Fig. S3, ESI†).

We obtain stability limits, and thus P_c and deformation at cavitation, Λ_c, from the expression for the free energy of the system (eqn (3)). Free energetic equilibria (dF_tot/dΛ = 0) correspond to the intersection of the void pressure and the gas pressure as a function of void deformation (Fig. 3), essentially amounting to a force balance. (The gas is replaced by a spring in the case of incompressible fluid loading, Fig. S3, ESI.†) In order for elastic cavitation to occur, there must be at least two equilibrium deformations such that it is thermodynamically favorable for the system to undergo a transition from the smaller to the larger deformed bubble while passing through a thermodynamically unstable (d²F_tot/dΛ² < 0) range of deformation.¹⁸

Fig. 3 illustrates how an increasingly compressed reservoir encounters two types of limits for cavitation that satisfy these criteria. First, the system encounters equilibrium (the absolute stability limit), at which the free energy of the smaller and larger bubble become equal (d²F_tot/dΛ² > 0, dF_tot/dΛ = 0, F_tot(Λ_small) = F_tot(Λ_large)), but are separated by an energetic barrier (d²F_tot/dΛ² < 0, dF_tot/dΛ = 0). Further compression leads to the metastability limit, in which the energetic barrier between the two states vanishes (d²F_tot/dΛ² = dF_tot/dΛ = 0). We use the latter to predict P_c as thermal fluctuations provide little perturbation to overcome an energetic barrier in a macroscale system. Thus, we determine the ability of the system to elastically cavitate for any given criteria, γ/μA, P₀/μ, and V_c0/(4/3πA³) (or kA³/μ for the incompressible fluid).

Non-zero surface energy is required in order to produce an elastic cavitation response in both loading condition cases. The minimum amount of surface energy is determined by the dimensionless ratio, γ/μA, which must be greater than ≈1.4953, the value at which a peak in the pressure (P)–deformation (Λ) curve (dP/dΛ = 0) no longer exists for real values of Λ (for an incompressible, neo-Hookean solid). (Examples of P vs. Λ for which a peak exists (black dashed) and at the limit where the peak disappears (gray dashed) are shown in Fig. 3a.)

When elastic cavitation is possible, we calculate P_c, Λ_c (the deformation at the onset of elastic cavitation), and Λ_cf (the final deformed state after cavitation). (The latter are shown with gray arrows in Fig. 3 where Λ_cf ≈ 11.) This model also allows us to determine the equilibrium void radius as a function of increasing compression, CF, as shown by the progression of increasing compression fraction curves in Fig. 3a. With this series of equilibirium deformation values, we replicate the experimentally observed pressure vs. compression curves displaying a sudden, discontinuous drop in pressure at P_c (Fig. 1 and Fig. S1, ESI†). Also notable from this analysis is the fact that P_c for elastic cavitation is well approximated by the simple ‘cavitation equation’,

proposed by Kundu and Croby² due to the fact that Λ_c for the metastability limit occurs on or near the peak in the pressure–deformation curve.

3 Fracture model

As a first approximation, we employ a fracture model proposed by Gent and Wang¹⁴ and later revisited by Diani²³ and Lin and Hui¹⁵ in which Griffith criteria are applied to a hyperelastic spherical void (neo-Hookean in their case). Following Diani's²³ treatment, the crack surface is taken to be the change in surface area of the void after the void has reached some critical radius at which fracture occurs,²³i.e., a “spherical” crack like that originally proposed by Williams and Schapery.²⁴ This approach is not intended to imply the actual formation of a spherical crack, but rather we aim to approximate failure of the bubble as an inner-shell damaged zone surrounding the surface of the sphere and having the complete loss of any elastic contribution. We ignore effects that the embedded needle might have on fracture geometry based on observations of the fractured surface in very soft solids. In very soft gels (μ < 10 kPa), fissures form a web throughout the ‘bubble’ surface with the needle typically centered within the bubble. We take this to imply that fracture must be initiated from several locations at once. We should note, however, that scattered observations seem to indicate fracture geometry is dependent on the stiffness of the material. Stiffer gels and elastomers typically have more penny-shaped ‘bubbles’ (μ ∼ MPa) or what appears to be delamination up the sides of the needle based on the final tip position at the bottom of a fractured void. As we are focused on the behavior of softer gels, we utilize a spherical approximation here. However, Lin and Hui¹⁵ found that the analytical spherical crack geometry relation actually provides a better estimate of fracture energy for the growth of a penny shaped crack than does Gent's penny-shaped crack approximation. Thus, the spherical crack geometry provides both a simple first approximation and would seem to apply to two elastic moduli extrema for CR (<10 kPa and ∼1 MPa). A more detailed study of the crack geometries that form and how they affect the measured P_c is part of an ongoing study.

The model we implement here, following Gent and Wang¹⁴ and Williams and Schapery,²⁴ is similar to Griffith's analysis, with a strain energy release rate G for an incompressible, neo-Hookean solid:

calculated from the derivative of the free energy of the void (F_surf + F_elas) with respect to a change in the unstrained void radius, A₀. A small change in A₀ is equivalent to a small change in the thickness of the damaged zone. Setting G equal to the material's critical strain energy release rate G_c, we calculate a critical stretch ratio at which fracture occurs, Λ_F, for a given hyperelastic material model and normalized critical strain energy release rate, G_c/μA. As expected for a neo-Hookean model (dashed lines, Fig. 4), increasing G_c/μA (increasing color saturation) corresponds to increasing stretch at fracture, Λ_F. Substituting Λ_F into the pressure–stretch relation for the bubble yields P_c for a fracture event. The final size of the fractured void would depend on the free energy of the system as a whole and likely the dynamics of the fracture itself; thus it is left for future work.

Fig. 4 Stretch ratio at fracture, Λ_Fvs. γ/μA for increasing G_c and varying constitutive response. Blue dashed lines: neo-Hookean, G_c = 2, 20, 200, 2000γ (in order of increasing color saturation). Blue solid lines: modified Yeoh, C₃ = 0.1, G_c = 2, 20, 200, 2000γ (in order of increasing color saturation). Black lines: modified Yeoh, C₃ = 0.1, G_c ∼ 2μ², A = 10, 100 μm (dotted and dash-dotted, respectively). Gray lines: modified Yeoh, C₃ = 0.1, G_c ∼ 4μ⁴, A = 10, 100 μm (dotted and dash-dotted, respectively).

Λ _F depends on the constitutive response of the material and G_c. We address the motivation behind varying these elements in Sections 7 and 8, however we give a brief overview of their effect on Λ_F here. The incorporation of strain hardening, decreases Λ_F overall as a function of γ/μA (solid lines). Fig. 4 also demonstrates that the monotonic dependence of Λ_F on γ/μA disappears when G_c is no longer assumed to be constant over order of magnitude changes in γ/μA (due to a dependence on μ). Assuming that G_c increases with increasing μ, it follows that Λ_F eventually increases with decreasing γ/μA as shown.

4 Mechanism map

Combining the elastic and fracture models outlined above, we build a predicted mechanism map for any given CR test or any series of tests (e.g., Fig. 5). This is accomplished by first considering whether or not an elastic cavitation event is possible. Within this regime fracture may also occur. The relative magnitude of Λ_F with respect to the stretch at cavitation, Λ_c determines the predicted mechanism. The magnitude of Λ_F falls into one of three categories: Λ_F < Λ_c (fracture before cavitation); Λ_c < Λ_F < Λ_cf (cavitation followed by fracture); or Λ_F > Λ_cf (cavitation only). The vertical dashed lines in Fig. 3 provide examples of these mechanism categories for three hypothetical Λ_F values. Outside of the cavitation regime, fracture is the only possible instability mechanism. The region in Fig. 5 where no symbols are shown indicates that the energy storage capacity of the system (either gas reservoir or spring) is too small to either overcome the peak in the pressure–deformation curve and/or too small to deform the void to Λ_F. Following this categorization procedure, we determine the dominant deformation mechanism for different combinations of materials properties, void sizes, and compression volumes. These results are represented as a deformation mechanism map for an incompressible, neo-Hookean solid in Fig. 5. A mechanism map must hold at least one length scale constant, either elastocapillary (by holding μ constant since γ varies little) or needle/bubble initial radius A. This arises from the fact that the dimensionless groups in eqn (3) and (6) have non-separable dependencies on both μ and A. Thus, a single map corresponds to a series of experimental tests either of varying material composition/condition (μ) and constant needle radius or varying needle radius with constant material composition/condition (μ) (Fig. S4–S6, ESI†). Therefore all deformation scenarios may not be present in all mechanism maps.

Fig. 5 Map detailing the mechanism leading to observation of a P_c during a CR experiment. ‘×’ corresponds to fracture, ‘●’ to cavitation, ‘■’ to cavitation followed by fracture, and ‘★’ to a negligible pressure drop during cavitation such that P_c is governed by fracture. The magnitude of P_c is shown via color gradient corresponding to the colorbar on the right. This map was generated assuming a constant needle radius of 50 μm, G_c = 20γ, γ = 27.7 mN m⁻¹ (air-2-ethyl hexanol interface), and P₀ = 1 atm for a neo-Hookean solid. (Grey line corresponds to V_c0 = 1 mL.)

Most CR experiments have been performed with relatively large gas volumes, putting them toward the top of Fig. 5, above the gray line that corresponds to V_c0 = 1 mL. As discussed above, the elastic cavitation limit that exists at finite γ/μA due to the loss of a peak in the P vs. Λ response (Fig. 3a) clearly divides the mechanism governing P_c. Needle radii and/or μ must be small for P_c values to have a chance to fall within the elastic cavitation regime.

An implication for failure not only in CR but in soft solids in general that follows from this work is that even when deformation is accommodated by elastic cavitation events, it is likely that these elastically initiated events lead almost instantaneously to localized fractures. Whether elastic cavitation can occur without fracture depends on the compliance of the loading environment, e.g., the chamber stiffness k or volume V₀, sample size (in the case of dead-loading),¹⁸ or machine compliance,¹⁹ with less compliance/smaller samples leading to an increased likelihood of microfracture. The extent of the fracture that occurs is dependent on dynamics and fracture geometries not studied here.

5 Correlation with experiments

To compare the predictions of Fig. 5a and the associated model with experiments, a series of CR experiments were performed in a thermally reversible triblock co-polymer gel using a single needle radius. The temperature dependent (15–40 °C) P_c data collected from the gel (7 wt% PMMA–PnBA–PMMA with block lengths of 25, 116, and 25 kg mol⁻¹ (Kuraray) dissolved in 2-ethyl-1-hexanol) using a 109 ± 3.4 μm radius needle are plotted in Fig. 6. To demonstrate both the sensitivity of the CR method and the variability that arises from changes in microstructure, we also include two sets of constant needle radius data taken from the literature, HSA Organogels¹¹ (low molecular mass gelators formed via non-covalent assembly, 30–80 °C, 0.25 wt%) and PAAm hydrogels² (cross-linked, volume fraction 0.022–0.088). Modulus values are obtained from shear rheology (triblocks, HSA) or contact mechanics (PAAm) data. The solid black line in Fig. 6 denotes the cavitation equation (eqn (5)). An isochor from Fig. 5 corresponding to V₀ = 1 mL leads to the prediction of the gray connected line in Fig. 6 (small filled symbols indicating cavitation and fracture as in Fig. 5). Note how closely the data falls to the cavitation equation, although scatter in the data (illustrated with standard deviation error bars) clearly show the cavitation equation is not an exact fit. Thus, while the cavitation equation can generally predict μ (to within 15%) for a soft solid, the CR technique is clearly sensitive to other mechanical properties governing material behavior and P_c response.

Fig. 6 P _c vs. γ/μA for the three datasets: (purple-diamonds) triblock copolymer gels (PAA–PnBA–PAA), (orange-inverted triangles) HSA organogels,¹¹ and (blue-left triangles) PAAm hydrogels.² A neo-Hookean model prediction taken from Fig. 5a is shown with a gray line and symbols representing the governing mechanism as outlined in the Fig. 5 caption. The cavitation equation (eqn (5)) is shown with a solid black line.

While the neo-Hookean constitutive relation provides a starting place for exploring hyperelastic behavior, its entropic spring assumption breaks down near material failure, as the network is fully stretched. Fig. 6 demonstrates a clear deficiency in the neo-Hookean model. Several experimental values lie above both the cavitation equation and the isochor prediction for V_c0 = 1 mL. (The gap between these two predictions at high γ/μA arises purely from geometry; the cavitation equation accounts for a spherical-cap like bubble whereas our model is entirely spherical.) The neo-Hookean prediction (Fig. 6, gray line) uses G_c = 20γ, however, even for increasingly tough materials P_c can never exceed the cavitation equation whose magnitude aligns with the peak in the bubble's pressure–stretch relation (Fig. 3a). Thus, the data clearly support the necessity of accounting for strain hardening, especially for the triblock gels,²⁵ prior to failure. Further, as written, the equations in Sections 3 and 4 assume fixed G_c over the set of conditions spanned by a mechanism map. Such an assumption of constant G_c is accurate for a constant μ map, but in the case of order of magnitude variations in μ in a constant A map, G_c will not remain constant, as becomes clear on an Ashby plot for soft solid materials.²⁶ The next two sections detail the incorporation of both strain hardening and G_c variation.

6 Strain hardening

We utilize a modified Yeoh constitutive relation in order to include strain hardening at larger stretch.Here I is the first strain invariant of the right Cauchy–Green tensor and C₃ is a phenomenological parameter that, when increased, leads to an increased strain hardening slope at decreasing stretch values (Fig. S2, ESI†). This model is convenient in that the pressure–stretch relationship for a pressurized void may be expressed in a closed analytical form (as opposed to the Gent relation). One consequence of the addition of strain hardening is that the limit in γ/μA at which elastic cavitation becomes possible (typically leading to cavitation-governed P_c values) is sometimes superceded by fracture due to a decreasing Λ_F, complicating interpretation of experimental data (Fig. 8, star region).

Further, Fig. 7a demonstrates that strain hardening leads to P_c predictions that lie above the cavitation equation threshold. Fig. 7a is generated using the assumption of constant G_c but spans a range of γ/μA values while ‘needle radius’ (corresponding to initial void radius) is held constant at A = 100 μm. The result of this assumption is that a negative slope in P_cvs. γ/μA (for constant A) cannot be predicted. In order to do so, and therefore capture the behavior observed in the triblock data shown in Fig. 6, we must consider the possibility of a G_c dependence on μ. (Experimental conditions are critical for interpreting P_c data. If μ, rather than A is held constant, no negative slope in P_cvs. γ/μA is observed, Fig. S10, ESI.†)

Fig. 7 Normalized P_c values predicted for varying constitutive behavior and G_c for a parameter series at constant needle radius (varying modulus). The cavitation equation (eqn (5)) is the black solid curve in both plots. (a) Effect of increasing strain-hardening corresponding to an increase in C₃/μ. Light and dark blue curve illustrates associated effect of increases in G_c (increased P_c) combined with strain hardening. (A = 100 μm, V_c0 = 1 mL, γ = 27.7 mN m⁻¹). (b) Predictions after incorporation of the power law dependence for G_c (eqn (9)). Increases in α tend to shift the P_c within the fracture regime up. Increases in β tend to make the slope within the fracture regime negative. (A = 100 μm, V_c0 = 1 mL, γ = 27.7 mN m⁻¹, C₃/μ = 0.01.)

Fig. 8 Maps detailing the mechanism leading to observation of a P_c during a CR experiment using the modified Yeoh Model. ‘×’ corresponds to fracture, ‘●’ to cavitation, ‘■’ to cavitation followed by fracture, and ‘★’ to a negligible pressure drop during cavitation such that P_c is governed by fracture. The magnitude of P_c is shown via color gradient corresponding to the colorbar on the right. This map was generated for a modified Yeoh material with C₃ = 0.01 and the same conditions as in Fig. 5.

7 Determining G_c

Measurements of P_c for fracture events should provide a means of determining G_c. If we assume linear fracture mechanics, then the relation between G_c and P_c for this spherical geometry is:²Using this approximation, G_c determined for the triblock, PAAm, and organogel systems are estimated as 0.4–5 J m⁻², 0.7–4 J m⁻², and 6–10 J m⁻², respectively (Fig. S11, ESI†). In all cases, these numbers appear to deviate substantially from those reported elsewhere. For the triblock gels, the maximum G_c,lin (at the lowest temperature) is half that of values reported in the literature²⁷ for similar weight fraction materials (∼10 wt%). Similarly, we would expect the PAAm system to have a G_c almost 5 times larger in order to be comparable with values obtained using traditional methods.³⁰ There are no literature values to compare to in the case of organogels, however, while a value of 10 J m⁻² puts organogels on par with typical hydrogels in the literature, those hydrogels typically have 10 or more times the polymer concentration, making this estimate unrealistic. A possible reason for this inability to quantitatively correlate P_c and G_c is that deformation at failure is highly non-linear and changes in a polymer network structure often leading to changes to both μ and G_c. In determining the predictions in Fig. 6 and 7a, we assumed that G_c would remain constant as μ changed, but this assumption is insufficient to capture experimental observations.

To account for this co-dependency, we can assume a relationship between G_c and μ, based on experimental data (e.g.ref. 26 and 31):

G_c = αμ^β (9)

This power law form also follows from idealized scaling arguments for crosslinked gels using a Lake–Thomas³² theory for G_c (see ESI†). We apply this relationship to determine G_c in the following way. First, we assume a particular material constitutive relation, in this case the modified Yeoh model (eqn (7)). For each of these systems, a fit of the stress–strain response would be necessary to get a precise value for the model parameter, C₃. However, as this data is not available in all cases, we instead approximate C₃ from average stretch at break (λ_F) using the fact that C₃ determines the onset of strain hardening (Fig. S2, ESI†). Thus C₃ is not a fitting parameter, but a way to account for the onset of failure observed experimentally. (Here we neglect changes in C₃ associated with changes in μ since none of the mechanisms by which μ is changed in each of these systems has a strong effect on the polymer chain length between crosslinks, the physical origin for finite extensibility.) Finally, we fit α and β to experimental data as shown in Fig. 9. G_c is calculated using the fit values of α and β and the independently determined μ. These C₃ values, corresponding two fitting parameters, and the G_c values obtained in this way are summarized in Table 1. G_c values obtained in this way either compare favorably to values obtained elsewhere in the literature (triblock, PAAm) or provide a more reasonable estimate of G_c in the case of the HSA organogels. The fit value of β = 0.1 for the PAAm gels is further supported using the scaling arguments of Lake–Thomas and experimental observation of the dependence of μ on polymer volume fraction (see ESI,† Section 5).

Fig. 9 Model fits to experimental data incorporating both strain hardening and power law G_c dependency. Fixed parameters in each case were obtained from either experimental conditions* or estimates from the literature^†: triblock – A* = 54.5 μm, V_c0* = 1 mL, γ^† = 27.7 mN m⁻¹, C₃/μ^† = 0.01; HSA – A* = 206 μm, V_c0* = 10 mL, γ^† = 40.9 mN m⁻¹, C₃/μ^† = 10; and PAAm – A* varied, V_c0* = 10 mL, γ^† = 72.0 mN m⁻¹, C₃/μ^† = 0.1. Power law parameters α and β were fit by eye.

Table 1 G _c predictions arising from Fig. 9

System	C 3 (λF)	α	β	G c [J m^−2] range	G c from the literature
Triblock (7 wt%)	0.01 (2–4^7,25)	3 × 10^−7	3.0	2–20	For 10 wt%, 10–20 J m^−2 (ref. 27)
HSA organogels (0.25 wt%)	10 (1.06^28)	0.8	0.08	∼1	None
PAAm (vol. frac. 0.027–0.088)	0.01 (5^29)	4	0.1	7–10	For vol. frac. ∼0.08, 10–20 J m^−2 (ref. 30)

---

It is important to note that for the modified Yeoh relation used here, the fit value α and thus G_c is dependent on C₃; while the parameter β is relatively insensitive to C₃. This makes the accurate choice of constitutive relation an important part of the interpretation of a quantitative G_c. For this reason our reported G_c values should still be taken as estimates, though, as stated above, we believe these estimates to be more accurate than using the linear approximation. To clarify the role that β plays in replicating experimental results, one must understand its effect on the deformation at fracture, Λ_F. For increasing power, β, Λ_F becomes increasingly large at low γ/μA (Fig. 4). This increased Λ_F corresponds to an increased P_c in a strain hardening material model as shown in Fig. 7b. For very large values of β, Λ_FA³ becomes large with respect to the initial chamber size V_c0 and our model predicts that no fracture and thus no P_c will be observed.

As an aside, the unexpected cubed relation between G_c and μ that is responsible for the large deviation of the triblock P_c data (Fig. 9) from the cavitation equation (the initial driving force for this work) is likely due to the thermally reversible nature of this material's crosslinks. The glass transition temperature, T_g for the PMMA endblocks is just above room temperature, likely resulting in the large change in P_c as a function of temperature right around 20–25 °C (Fig. S8, ESI†) that leads to β = 3.0. We interpret this as the result of a sharp decrease in chain pullout forces over a region of constant elastic response due to relatively unchanged chain configuration between crosslinks (Fig. S9, ESI†). The G_c dependence on temperature is significantly stronger than has been predicted via simulations of pull-out from a purely glassy state,³³ suggesting a transition from pull-out in a condensed state near the gel point to pull-out from a glassy state at lower temperatures.

8 Conclusions

In summary, we outline a framework for understanding the physical origins of CR and, more importantly, provide the means by which to acquire quantitative fracture energies from that data. Our side-by-side comparison of three different materials systems also illustrates the sensitivity of the technique to microstructural variation when series of datum (of varying composition or at varying environmental conditions) are acquired. In these three systems, ‘cavitation equation’ estimates for modulus have less than 15% error, making CR useful for a broad range of modulus measurement in geometrically unconfined samples. However, the importance of knowing which governing mechanism leads to an unstable pressure drop and when non-linearities are likely to affect mechanical response becomes critical when determining G_c. This sensitivity of CR to both scale and microstructure points toward the physical intuition that might be gained from using this technique for in vivo measurements of living tissues.

Nomenclature

a	Deformed bubble radius
A	Initial bubble radius (equivalent to needle radius)
A 0	Zero strain energy bubble radius
α	Power law coefficient
β	Power law power
C 3	Phenomenological parameter in Yeoh constitutive relation
CF	Gas compression fraction, 1 − Vc/Vc0
F tot	Total free energy of the bubble + system + surroundings
G	Strain energy release rate
G c	Critical strain energy release rate
γ	Surface energy
I	First strain invariant of the right Cauchy–Green tensor
k	Stiffness of the fluid pressurizing setup
λ F	Experimentally measured strain at failure
Λ	Deformation of the bubble: a/A
Λ 0	Hoop stretch; zero strain energy bubble deformation: a/A0
Λ c	Deformation of the bubble at Pc due to cavitation
Λ cf	Final deformation of the bubble after elastic cavitation
Λ eq	Zero strain energy bubble deformation due to surface energy: A/A0
Λ F	Deformation of the bubble at Pc due to fracture
μ	Linear elastic shear modulus
P	Pressure
P 0	Pressure of the environment
P c	Critical pressure at which sudden bubble expansion occurs
V c	Chamber volume
V c0	Initial chamber volume
W	Strain energy function

Acknowledgements

The authors would like to thank the Army Research Laboratory for funding through Cluster M of the Center for University of Massachusetts Industry Research on Polymers and NSF-DMR 1304724.

References

1. J. A. Zimberlin, N. Sanabria-DeLong, G. N. Tew and A. J. Crosby, Soft Matter, 2007, 3, 763.
2. S. Kundu and A. J. Crosby, Soft Matter, 2009, 5, 3963.
3. J. A. Zimberlin, J. J. McManus and A. J. Crosby, Soft Matter, 2010, 6, 3632.
4. J. A. Zimberlin and A. J. Crosby, J. Polym. Sci., Part B: Polym. Phys., 2010, 48, 1423–1427.
5. J. Cui, C. H. Lee, A. Delbos, J. J. McManus and A. J. Crosby, Soft Matter, 2011, 7, 7827.
6. A. J. Crosby and J. J. McManus, Phys. Today, 2011, 64, 62.
7. A. Delbos, J. Cui, S. Fakhouri and A. J. Crosby, Soft Matter, 2012, 8, 8204–8208.
8. M. S. Chin, B. B. Freniere, S. Fakhouri, J. E. Harris, J. F. Lalikos and A. J. Crosby, Plast. Reconstr. Surg., 2013, 131, 303e–305e.
9. L. Pavlovsky, M. Ganesan, J. G. Younger and M. J. Solomon, Appl. Phys. Lett., 2014, 105, 114105.
10. A. Blumlein and J. J. McManus, J. Mater. Chem. B, 2015, 3, 3429–3435.
11. P. Fei, S. J. Wood, Y. Chen and K. A. Cavicchi, Langmuir, 2015, 31, 492–498.
12. A. N. Gent and P. B. Lindley, Proc. R. Soc. London, Ser. A, 1959, 249, 195.
13. J. M. Ball, Philos. Trans. R. Soc., A, 1982, 306, 557–611.
14. A. N. Gent and C. Wang, J. Mater. Sci., 1991, 26, 3392–3395.
15. Y. Y. Lin and C.-Y. Hui, Int. J. Fract., 2004, 126, 205–221.
16. V. Lefevre, K. Ravi-Chandar and O. Lopez-Pamies, Int. J. Fract., 2013, 192, 1–23.
17. A. Hamdi, S. Guessasma and M. Nat Abdelaziz, Mater. Des., 2014, 53, 497–503.
18. J. Dollhofer, A. Chiche, V. Muralidharan, C. Creton and C.-Y. Hui, Int. J. Solids Struct., 2004, 41, 6111–6127.
19. V. Muralidharan, C.-Y. Hui, J. Dollhofer, C. Creton and Y. Lin, Int. J. Adhes. Adhes., 2006, 26, 117–124.
20. J. Zhu, T. Li, S. Cai and Z. Suo, J. Adhes., 2011, 87, 466–481.
21. S. B. Hutchens and A. J. Crosby, Soft Matter, 2014, 10, 3679–3684.
22. A. N. Gent and D. A. Tompkins, J. Polym. Sci., Part A-1: Polym. Chem., 1969, 7, 1483–1487.
23. J. Diani, Int. J. Fract., 2001, 112, 151–161.
24. M. L. Williams and R. A. Schapery, Int. J. Fract., 1965, 1, 64–72.
25. S. M. Hashemnejad and S. Kundu, Soft Matter, 2015, 11, 4315–4325.
26. J. Li, Z. Suo and J. J. Vlassak, J. Mater. Chem. B, 2014, 2, 6708–6713.
27. M. E. Seitz, D. Martina, T. Baumberger, V. R. Krishnan, C.-Y. Hui and K. R. Shull, Soft Matter, 2009, 5, 447.
28. J. L. Li, R. Y. Wang, X. Y. Liu and H. H. Pan, J. Phys. Chem. B, 2009, 113, 5011–5015.
29. J.-Y. Sun, X. Zhao, W. R. K. Illeperuma, O. Chaudhuri, K. H. Oh, D. J. Mooney, J. J. Vlassak and Z. Suo, Nature, 2012, 489, 133–136.
30. Y. Tanaka, K. Fukao and Y. Miyamoto, Eur. Phys. J. E: Soft Matter Biol. Phys., 2000, 401, 395–401.
31. H. J. Kong, E. Wong and D. J. Mooney, Macromolecules, 2003, 36, 4582–4588.
32. G. J. Lake and A. G. Thomas, Proc. R. Soc. London, Ser. A, 1967, 300, 108–119.
33. M. Bulacu and E. van der Giessen, EPL, 2011, 96, 49901.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5sm02055g

---