Growth of gas-filled penny-shaped cracks in decompressed hydrogels

Abstract

We report that the decompression of soft brittle materials can lead to the growth of internal gas-filled cracks. These cracks are oblate spheroids (‘penny shape’), whose major radius grows linearly in time, irreversibly fracturing the surrounding material. Our optical measurements in hydrogels characterise and quantify the three-dimensional crack geometry and growth rate. These results are in good agreement with our analytical model coupling fracture mechanics and gas diffusion, and predicting the dependence on the mechanical properties, gas diffusivity and super-saturation conditions (gas pressure, solubility, temperature). Our results suggest a new potential mechanism for decompression sickness in scuba diving and for indirect optical measurements of the fracture properties of hydrogels.

---

1 Introduction

It is well-known that dissolved gases can form bubbles following a reduction in solubility caused by a reduction in ambient pressure (decompression). This general phenomenon spans a variety of materials (liquids and solids) and applications, ranging from fizzy drinks (e.g. champagne¹), industrial high-pressure gas containers (cavitation damage of rubber seals²), and volcanic eruptions (fragmentation of rising magma³). A related phenomenon occurs in polymer networks saturated by liquid–liquid mixtures. Liquid droplets may form within the polymer following a reduction in miscibility caused by temperature^4–6 or composition,⁴ thus causing the polymer to deform in order to accommodate the new phase. However, perhaps the best-studied application of the gas–liquid phenomenon is the decompression sickness experienced by human scuba‡ divers. Gases breathed while at depth dissolve in body tissues (blood, skin, muscles, joints, nerves, etc.) at concentrations exceeding saturation at atmospheric pressure, and can form bubbles upon resurfacing. These bubbles can, for inert gases (chiefly N₂) and under certain ‘extreme’ circumstances, grow numerous and large enough to be painful, seriously pathological, and even fatal (see ref. 7 Section 2.2 for a short review of the complex pathophysiology).

Significant research and debate exist on the conditions under which bubble populations can form and grow; in particular on the questions of bubble nucleation at liquid/solid interfaces,^8–10 the competition between perfusion and diffusion in soft tissues,^11–13 and the role of tissue elasticity.^14,15 However, these studies have in common to assume that gas bubbles are spherical in shape, and that when they grow in soft materials, they do so by elastic (and therefore reversible) deformation. These assumptions are broadly consistent with the limited number of observations of small (≲ 1 mm) bubbles in hydrogels mimicking artificial tissues.^15,16 However, the forced syringe injection of liquid^17,18 or gas¹⁹ into hydrogels is known to cause large (≫ 1 mm) cracks that are non-spherical and that grow by irreversible fracture of the material.

In this paper we demonstrate and explain that the growth of large, non-spherical, gas-filled cracks can also generically occur as a result of decompression.

We study hydrogels super-saturated in dissolved CO₂ by swelling cross-linked polyacrylamide beads in carbonated water at modest absolute pressures (≈2–4 bar). Upon decompression to atmospheric pressure, no gas-filled bubble or crack forms spontaneously; instead we drop a hydrogel bead (weighing ≈10 g) directly onto a solid surface (from a modest height ≈0.1–1 m) in order to initiate or enlarge microscopic internal fractures (invisible to the naked eye). For sufficiently high concentrations of dissolved gas (proportional to the gas pressure during swelling) and impact heights, a number of gas-filled cracks (typically one to five) start growing inside the bead.

What follows is illustrated in Fig. 1. The spherical hydrogel bead is imaged almost fully immersed in water, against a high-contrast grid background (the edge of the bead is outlined in black). At time 32 s after impact, two millimetre-sized, circular thin cracks are clearly visible near the centre of the bead. Until time t = 82 s (second panel) these cracks grow independently of one another at similar speeds (≈2 mm radius gain in 50 s). Until t = 132 s, these ‘penny-shaped’ (or ellipsoidal) cracks grow in increasing proximity, slightly distorting each other. Over the next 10 s, they partially merge, before the ‘left’ crack deflates after reaching and breaking through the surface of the bead. As a result, the ‘right’ crack partially deflates as well, and both cracks undergo a few oscillations in size before coming to a relative rest. The whole process is evidently irreversible. The corresponding movie is available as ESI.†

Fig. 1 Growth and interaction of two gas-filled cracks in a spherical bead of polyacrylamide hydrogel upon decompression from ≈2.5 bar absolute CO₂ pressure to atmospheric pressure. The bead is held still and almost fully submerged in water (see blurred level on top) in order to highlight the gas-filled cracks by refraction. The very faint refraction between bead and water corresponds to the edge of the bead, outlined in black.

In the remainder of this paper, we focus on the more elementary problem of the growth of a single gas-filled crack in an apparently infinite elastic brittle material following sudden decompression and impact.

Our objectives and outline are as follows. In Section 2 we explain how we used optical measurements to accurately measure the three-dimensional shape of a crack during its growth. In Section 3 we derive a predictive model for the crack shape and growth from first principles. In Section 4 we validate this model with our experimental measurements, covering a range of pressures and mechanical properties. In Section 5 we discuss implications of our work and conclude.

2 Experiments

2.1 Protocol

Our experiments were conducted as follows.

First, bottled water was carbonised with pure CO₂ using a home soda maker, at variables absolute pressures in the range 2–4 bar (the water was first chilled to reach the highest pressures).

Second, a few commercial-grade cross-linked polyacrylamide dry beads (Soil Moist^TM, JRM Chemical, USA, as in ref. 20) were swollen at room temperature with this carbonated water in a sealed bottle for at least 48 hours.

Third, the seal cap was broken, and swapped with a pressure gauge cap. The bottle was vigorously shaken in order to stabilise the pressure of the small pocket of CO₂ above the water, which we take as the equilibrium ‘super-saturation’ pressure p_s to determine the dissolved gas concentration in the beads c_s = k_hp_s (where k_h is Henry's law solubility constant).§ The water temperature inside the bottle was measured within 0.5 °C and used as the representative experimental temperature T of the hydrogel.

Fourth, a bead was taken out of the bottle, and its diameter, mass, and elastic Young modulus E were measured. For the latter we used the standard Hertz contact model and fitted the force/displacement curve obtained by compressing the spherical bead between two rigid glass plates using an accurate translation stage and balance (for more details about this method, see the ESI,† Section 2). Inherent variability in the polymer/water volume fraction among the beads caused a range of E = 12–28 kPa, which we will use to our advantage in order to validate our model in Section 4.

Fifth, directly afterwards, the bead was individually dropped from ≈0.2 m onto a flat glass plate (as described in the introduction).¶ The bead was then inspected by eye, and, if no visible gas-filled crack had formed, it was dropped again from incremental heights, until a single crack formed or two cracks formed far enough from one another to remain largely independent for sufficiently long times. If the crack was too close to the surface of the bead, or if multiple cracks formed too close to one another, the bead was discarded, another bead was taken out of the bottle, and these steps were repeated.

Sixth, the ‘successful’ bead was swiftly positioned in a (non-carbonated) water bath at the same temperature T as the bead. Imaging started typically at t ≈ 5–20 s (recorded from the moment of the successful impact) and continued until the crack reached the surface of the bead or became visibly influenced by it.

2.2 Imaging

Image acquisition. Our imaging setup, sketched in Fig. 2a, consisted of four mirrors arranged in a ‘theatre’ using a 3D-printed stand such that they are at 45° from the horizontal, and at 45° from one another.

Fig. 2 (a) Experimental setup to capture four views at 45° of a single crack growing inside the bead. The camera (not displayed here) is placed directly over the bead. The water level (not displayed here) is such that the mirrors are fully submersed, and the bead almost submersed, as in Fig. 1 (the bead is not fully submersed to prevent it from floating). (b) Snapshot illustrating the raw data obtained by the camera, ready for analysis in Fig. 3 (the edge of the bead is outlined in black).

This allowed a single video camera, which had a top view of these mirrors, to capture simultaneously four different and complementary views of the growth of a single ‘isolated’ crack, as shown in Fig. 2b. These four views were then combined during image analysis to reconstruct the time-resolved, three-dimensional geometry of the crack.

A bright and high contrast colour background was displayed by LCD screens placed directly across each mirror (see Fig. 2a). To avoid undesirable reflections and loss of background brightness, an octagonal water tank was used to immerse the bead and mirrors.

Images were captured by an 8.3 MPixels video camera with good depth of field and negligible parallax. The resulting spatial resolution was 40 μm Pixel⁻¹ (although the combination of four views could achieve sub-pixel resolution). The frame rate of 25 Frame s⁻¹ was sufficient given our growth rates in the range 0.5–4.5 Pixel s⁻¹, corresponding to 0.02–0.18 Pixel Frame⁻¹.

Image analysis. The raw images were analysed according to the following steps, summarised in Fig. 3.

Fig. 3 Key steps of the image analysis, from (a) the acquisition of raw images (four views at 45° angles) to (f) the three-dimensional ellipsoid best fitting the crack (more details in the ESI,† Section 3).

The raw images (panel a) first underwent background subtraction (panel b), binarisation (panel c), and contour detection (panel d).

These raw contours were invariably noisy, and were thus filtered in two steps (panel e). In step 1, their convex hull was computed, shrunk to 80% (shown in red shade), and any points inside this region were discarded (shown in red). In step 2, two-dimensional ellipses were fitted to the remaining points, and the points furthest away from the fitted ellipse were in turn discarded (shown in red). A further set of two-dimensional ellipses were fitted to the final set of remaining points, which were then centred based on this last two-dimensional fit (this centering is crucial in order to ensure that all views have a common origin for the next step).

Finally, a three-dimensional ellipsoid was fitted (panel f) by minimising the distance between the four projections of this ellipsoid (at 45°) and the four sets of remaining points centred on a common origin. Our choice of fitting a three-dimensional ellipsoid to the crack was influenced by our observations (see Fig. 1), as well as theoretical expectations (as will become clear in Section 3.1). The adequacy of this choice is confirmed by the excellent agreement between the two-dimensional projections of the fitted ellipsoid and the four filtered contours (see panel f, right hand side).

The absolute position of the crack within the bead could be and was computed, but it is not discussed in this paper since we focus exclusively on the crack geometry. More details about the image acquisition and analysis are given in the ESI,† Section 3.

2.3 Qualitative observations

Fig. 4 shows the results obtained from image processing in a typical experiment (quantitative results for this particular experiment will be discussed in Section 4.1).

Fig. 4 Results of the image analysis showing (a) snapshots of the three-dimensional ellipsoids best fitting the crack (at t = 24.5, 49.0, 73.5, 98.0 s); (b) continuous-time two-dimensional projections (ellipses) in the plane of the major axes (top view) and perpendicular to it (side view). The four snapshots of (a) are highlighted in white in (b). The straight dashed lines show that linear growth in time is a good approximation along the major axes (major radii) but not along the minor axis (thickness).

In panel a, we plot four snapshots (i = 1–4) of the three-dimensional ellipsoid best fitting the crack geometry at t = iΔt where Δt = 24.5 s (for visualisation purposes, the ‘base’ of successive cracks is arbitrarily set as the z = −3i plane). In panel b, we reveal further details by plotting two-dimensional projections of the crack: a ‘top view’ (x–y plane) and a ‘side view’ (y–z plane). These ellipses are plotted successively so that they create a continuous spatio-temporal trace, whose edge allows us to infer the growth rate (the four snapshots of panel a are also highlighted in white in panel b for comparison).

First, we confirm the observations of Fig. 1 that the crack is not spherical; it is circular but thin. In other words it has two approximately equal major axes (here along x,y by convention) and a much smaller minor axis (here along z).

Second, sets of straight dashed lines prove to be excellent guides to the growth of the crack along its major axes (panel a and top view in panel b), suggesting linear growth of the radius in time. However, straight lines prove to be poor guides to the growth of the crack along its minor axis (side view in panel b), suggesting that the thickness grows with sub-linear scaling (the growth slows down with time). In other words, the aspect ratio is not conserved; the crack becomes relatively thinner with time.

3 Model

In this section we introduce an analytical model to account for the above observations and provide quantitative predictions. We establish the link between crack geometry and gas pressure within the crack in Section 3.1, add the diffusion of gas into the crack in Section 3.2, and finally deduce its temporal growth in Section 3.3. The model is sketched in Fig. 5 and the key constants and variables used in the remainder of the paper are summarised in Table 1.

Table 1 Nomenclature used in the model. Some variables depend on the crack radial and axial coordinates (r,z) and/or on time t. In the last column (value), we distinguish between direct experimental measurements and indirect estimations using these measurements together with model assumptions. Note that the growth rate α is both measured and estimated, allowing for the model validation in Section 4. For conciseness, we denote values of indirect interest by −. The values of k_h, D within the gel are assumed to be identical to those in pure water, and their temperature dependence near T = 298.15 K are obtained respectively from ref. 21 (p. 4488 and eqn (19)) and ref. 22 (pp. 6–261, after fitting)

Name	Description/definition	SI unit	Value
p a	Atmospheric pressure	Pa	1.013 × 10^5
p s	Partial gas pressure in vessel during the swelling of the hydrogel	Pa	Varied (2.3–4.1 × 10^5)
p b(t)	Absolute pressure inside growing bubble (uniform)	Pa	—
Δp(t)	Pressure in bubble relative to atmospheric = pb(t) − pa	Pa	—
R(t)	Radius of the penny-shaped (spheroidal) crack (R ≪ R∞)	m	Measured (Fig. 4 and 6)
w(r,t)	Width of the crack	m	Measured (Fig. 4)
A(t)	Area of the crack	m^2	Measured (Fig. 4)
V(t)	Volume of the crack	m^3	Measured (Fig. 6)
E	Young modulus of the gel	Pa	Measured (Fig. 7)
ν	Poisson ratio of the gel	—	0.5 (incompressible)
K I(t)	Stress intensity factor at the crack tip (‘opening’ mode I loading)	Pa m^1/2	—
K IC	Critical stress intensity factor (mode I), also called material toughness	Pa m^1/2	Estimated by V, R, E (Fig. 7)
G C	Critical strain energy release rate (mode I) = 3KIC^2/(4E)	Pa m	Estimated by V, R, E (Fig. 7)
[small script l]	Elasto-fracture length of the gel = 9KIC^2/(4πE^2) = 3GC/(πE)	m	Estimated by V, R (Fig. 7)
	Ideal gas constant	J K^−1 mol^−1	8.314
T	Temperature during swelling of the gel and crack growth	K	Measured (295–299)
k h	Henry's solubility constant in water (here for CO2)	mol m^−3 Pa^−1	3.4 × 10^−4 exp{2400(T^−1 – 298.15^−1)}
D	Diffusivity of the gas in water (here for CO2)	m^2 s^−1	1.9 × 10^−9+ 4.8 × 10^−11(T – 298.15)
c(r,z,t)	Dissolved gas concentration in the hydrogel, i.e. (r,|z|) > (R,w)	mol m^−3	—
c b(r,t)	Dissolved gas concentration at the crack boundary = khpb(t)	mol m^−3	—
c s	Dissolved gas concentration far from the crack = khps	mol m^−3	—
Δc	Gas concentration difference driving the diffusion = kh(ps − pa) > 0	mol m^−3	—
S	Non-dimensional super-saturation ‘strength’ = [scr R, script letter R]TΔc/pa	—	varied through ps (1.07–2.67)
α	Growth rate of the crack radius = dR/dt	m s^−1	measured and estimated by D, S, [small script l] (Fig. 8)

---

3.1 Fracture mechanics

Assumptions. Based on our previous observations, we assume that the crack is growing purely through ‘mode I’ fracture, i.e. through opening of the surrounding material under a tensile stress equal to the net pressure difference Δp(t) = p_b(t) − p_a. The absolute gas pressure inside the crack is p_b(t) and p_a is the atmospheric pressure outside. The crack is assumed to be axisymmetric (with cylindrical coordinates (r, z) centred on the crack), and grows in a homogeneous, isotropic material with Young modulus E and Poisson ratio ν. The surrounding material is assumed to be infinite, i.e. the crack remains sufficiently far from the edge of the hydrogel bead such that it does not feel its influence. The crack is also assumed to be sufficiently thin such that the stress field at its boundary can be approximated by the uniform axial tensile stress |σ_zz| = Δp. We further assume that the crack propagates slowly, in quasi-static mechanical equilibrium.

Geometry. The problem described by the above assumptions was first solved in ref. 23 (although we prefer the more recent and concise treatment in ref. 24 Section 3.6.2). The ‘penny-shape’ crack solution is an axisymmetric|| ellipsoid of major radius R(t) and width w(r,t) < R given by

Mathematically, this ‘penny’-shaped crack is an oblate spheroid (sketched in Fig. 5a and b). Since our hydrogel is primarily made of (incompressible) water, we use ν = 0.5 throughout the remainder of the paper.

Fig. 5 Model sketch. (a and b) Top and side view of the general crack geometry. (c) Axial diffusion (along z) of the gas concentration between the ‘far-field’ c_s and the crack flat boundary c_b ≈ c_s − Δc, the gas flux being proportional to the gradient ∂c/∂z at the growing boundary. (c) Principle of the time integration of the diffusive flux during crack growth as explained in (13)–(16). In (c and d) the bottom half of the crack z < 0 is implicit by mirror symmetry.

Propagation. The stress field caused by the crack pressure Δp and displacement R, w concentrates at the crack tip (r = R, z = 0), which controls fracture propagation. Linear elastic fracture mechanics uses the local stress intensity factor K_I (the subscript ‘I’ indicates ‘mode I’ fracture), which for this penny-shape solution is equal to

The crack propagates when K_I reaches the critical stress intensity factor (or toughness) K_IC, a material property.** This quasi-static propagation criterion allows us to link, at all times, crack pressure and radius as

Inserting (3) into (1), we link the crack geometry to the mechanical properties

Volume. From (4), the volume of the crack is

We recognise that this fundamental V(t) ∝ R^5/2(t) scaling relation predicted for penny-shaped cracks can be equivalently written as

In other words, the maximum width (thickness) of the crack is the geometric mean of R and [small script l], the elasto-fracture length²⁷ of the material defined as

using either the critical stress intensity factor (toughness) K_IC or the critical strain energy release rate G_C (see previous footnote).

This penny-shaped crack geometry has been used extensively, in particular in the geophysical literature on the hydraulic fracturing of rocks.^17,18,28

Validity. We note that (6) is only consistent with a penny shape (oblate spheroid) provided that the minimum radius considered is large enough

We interpret this as the main condition under which the above fracture mechanics model is valid, and under which ‘large’ (R > [small script l]) gas-filled penny-shaped cracks are expected, as opposed to the ‘small’ (R < [small script l]) spherical bubbles considered elsewhere in the literature.^8–15

3.2 Gas diffusion

We now add the diffusion of dissolved gas, as sketched in Fig. 5c, to obtain a dynamical equation for the volume V(t) and the radius R(t).

Assumptions. Far away from the crack, we assume that the dissolved gas maintains a constant equilibrium molar concentration c_s = k_hp_s, where k_h(T) is the Henry constant for the given gas/liquid couple and temperature T (assumed constant and uniform throughout), and p_s is the constant ‘super-saturation’ pressure of the gas during the swelling of the hydrogel (because we use pure CO₂, this partial pressure is in fact the total absolute pressure, and of course p_s > p_a). At the crack gas–liquid boundary, the dissolved gas has concentration c_b(t) = k_hp_b(t). We assume a constant c_b ≈ k_hp_a, which is a good approximation for our experiments (in other words the excess pressure inside the crack Δp(t) is important for the deformation but negligible for the diffusion).†† The diffusion of gas through the liquid phase of the gel and towards the crack boundary is assumed to be purely axial (along z), driven by the concentration difference Δc = c_s − c_b ≈ k_h(p_s − p_a), with diffusivity D. Inside the crack, we assume the ideal gas law with constant .

Gas flux. We solve the one-dimensional diffusion equation by similarity variables, approximating the boundary at z ≈ 0 (consistent with the thin crack assumption in Section 3.1). We find the following molar flux of gas per unit area, on either side of the boundary:

We convert the above molar flux of dissolved gas into the following volume flux of ideal gas inside the crack:

Importantly, we note that it represents the flux per unit area at a time t after this unit area was created.

As the crack grows, new area is continually created, such that the instantaneous flux (10) across ‘new areas’ close to the crack tip (r ≲ R) is always higher than that across ‘old areas’ close to the centre (r ≈ 0) as sketched in Fig. 5c.‡‡ For a unit area created at t = 0, the total change in volume per unit area after time t is thus given by the integral

We define the convenient effective diffusivity D_eff, which can also be expressed in terms of the non-dimensional ‘super-saturation’ S that we define as

Time integration. To compute the total volume flux through the growing crack area, it is helpful to consider the following history, sketched in Fig. 5d:

• At t₀ = 0 the area is A(t₀) ≈ 0 and the volume is V(t₀) ≈ 0.

• At t₁ = δt (small) the crack surface is an infinitesimal disk of area and the crack volume is

• At t₂ = 2δt the cumulative flux through the ‘old’ area A₁ slowly increases (thicker arrow in Fig. 5d), while we now have an additional flux through the ‘new’ infinitesimal annulus of area The total volume is thus given by

• At t_i = iδt, we generalise this sum by mathematical induction

Note that t_i stays constant in the sum over t_j.

In the limit δt → 0 we recover the following convolution integral

3.3 Deducing the growth dynamics

We now assume that the crack grows with power law R(t) = αtⁿ, and we seek the growth rate α and exponent n.

Combining the V–R scaling relation of fracture mechanics (6), the total area A = 2πR² (thin crack), and the gas diffusion history (16), we find

Using the power law ansatz, we find

Since α is independent of t, the last equation requires n = 1, thus our model predicts linear growth

R(t) = αt. (19)

By evaluating the integral in (18) as = 4/15 we find the growth rate

We conclude that the growth rate is given analytically by the dimensional speed D[small script l]⁻¹ (gas diffusivity/material elasto-fracture length), modulated by S² (square of the non-dimensional super-saturation), and a pre-factor ≈0.81. This is the main result of our model.

Note that the time scaling of diffusion-driven cracks (dV/dt ∝ t^3/2 and R ∝ t) differs from that of fluid-driven cracks created by a constant injection rate (dV/dt = const. and R ∝ t^2/5 in the toughness-dominated limit).^17,18

Validity. We now re-examine our assumption in (9) that the diffusion in the gel is primarily axial, despite the presence of radial gradients due to the crack growth history (see Fig. 5c). The full two-dimensional axisymmetric diffusion problem does not admit an analytical solution because it requires mixed boundary conditions at z = 0 (fixed c = c_b for r < R and no flux for for r > R by symmetry) as well as moving boundaries R(t) = αt.

However, by a simple advection–diffusion scaling argument we consider the ratio of maximum length-scales achieved by advection ατ and diffusion where τ is the total experimental run time. The resulting Péclet number is (in terms of maximum crack radius achieved R_max = R(τ)). Using our experimental values (discussed later), we find a range of Pe ≈ 40–800, which is indeed high enough (≫1) to neglect radial gradients.

Numerical solutions of the full problem confirmed that even for values as low as Pe = 15 (corresponding to a crack growth slower and longer than in any of our experiments) purely axial diffusion under-estimates the actual final gas flux across the crack boundary by less than 1%. For more details on these scaling arguments and numerical solutions, see the ESI,† Section 4.

Using (19) we recast this Pe ≫ 1 validity condition of purely axial diffusion as

R_max ≫ [small script l]S⁻² (21)

In other words, the maximum radius achieved by the crack must greatly exceed [small script l]S⁻², otherwise radial diffusion will noticeably reduce the gas flux, and therefore the pre-factor in (20).

4 Results and model validation

We now turn to experimental results to validate the model, which can be summarised by four equations only: the shape of the crack (6); the relation between shape and mechanical properties (7); the linear growth of its radius (19); and the dependence of the growth rate on system parameters (20). We examine (6) and (19) in Section 4.1, (7) in Section 4.2, and (20) in Section 4.3.

4.1 Crack shape and linear growth

To illustrate our analysis of the data, we focus in this section on a single representative experiment (identical to that of Section 2.3 and Fig. 4).

R–t scaling. In Fig. 6a we plot the crack radius against time, using two different measures. In orange we plot the mean and spread between the minimum and the maximum radii obtained from the four raw two-dimensional contours (see Fig. 3d). The small spread indicates that all views give a consistent measure of the largest dimension of the crack (i.e. the radius R). In black we plot the mean and spread of the two major radii of the final three-dimensional ellipsoid fitted to these four contours (see Fig. 3f). The small spread between these two radii confirm that the crack is very nearly axisymmetric (i.e. spheroidal), as assumed in the model.

Fig. 6 Validation of the crack shape and linear growth in a representative experiment. (a) R–t scaling using the four ‘raw’ 2D contours (in orange) and the two major radii of the 3D ellipsoids (in black), and their respective linear fits. (b) V–t scaling (log–log scale) using the 3D ellipsoids, and its power law fit. Here R is taken as the mean of the two major radii (black symbols in (a)). The symbols and error bars represent the mean and standard deviation over 12 successive frames in which the crack radius grows by less than 1 Pixel. In the inset (in blue) we plot the corresponding values of [small script l]; the solid and dashed lines represent the mean and standard deviation over the range of R.

Both measures of the radius are in good agreement and grow linearly to a very good approximation (as anticipated in Fig. 4), which validates the key model eqn (19). Linear fits (dashed lines) give similar growth rates α = 0.0410 mm s⁻¹ (from here on quoted using the more direct ‘orange’ contour data, but within a few percents of the ‘black’ ellipsoid data). In this and other experiments, the standard error of the linear fit is typically negligible (<0.1%), such that we consider this measured α virtually exact. We defer the comparison between this value and that predicted by the model in (20) to Section 4.3, because the predicted α will suffer from errors in S and [small script l] that we have not yet addressed.

From the intercept of the fit R(t = 0) ≈ 0.4 mm we conjecture that, immediately upon impact, the radius was <0.4 mm and grew more rapidly than αt (i.e. sub-linear in time). However, more detailed observations would be needed to clarify the early-time behaviour, which falls outside the scope of this paper.

V–R scaling. In Fig. 6b we plot the crack volume against radius during the same interval of time t = 18–98 s. We use the volume of the three-dimensional ellipsoid and the mean radius along its two major axis (also shown in black in panel a). We also use the high temporal resolution of our data to provide ‘sub-pixel’ error estimates to V and R. We do so by computing and plotting the average and standard deviation across ≈1 s time windows (corresponding to 12 frames, after analysing only every second frame) with 50% overlap between successive windows. Across each such time windows, this crack grows only by ≈1 Pixel in radius, such that any observed variability in V and R can be chiefly attributed to random noise in the image analysis (shown as error bars).

We see that the crack geometry closely follows the penny shape scaling V ∝ R^2.5 assumed in the model, since the best power law fit gives an exponent 2.54. This is typical of our other experiments, which had an average exponent 2.55 with standard deviation 0.12. This validates the key model eqn (6).

[small script l] Value. The intercept of the fit, V(R = 1 mm) = 1.23 mm³, provides a first straightforward optical measurement of the elasto-fracture length of the material [small script l] = {(3 × 1.23)/(4π)}² ≈ 0.0862 mm (assuming the V ∝ R^2.5 scaling). More rigorously, we apply (6) directly to all the data and plot (shown in the inset in blue). The mean value across the full range of R is [small script l] = 0.0931 mm (solid line), and the standard deviation is 0.0044 mm = 4.7% of the mean (dashed lines). Henceforth we take these values as the measured [small script l] and its error (the error being primarily caused by apparent deviation from the exact scaling V ∝ R^2.5).

4.2 Mechanical properties

In Fig. 7 we present results from 70 independent experiments in which [small script l] was measured from R(t), V(t) data as explained in the previous section, E was measured as explained in Section 2.1 (and ESI,† Section 2), and K_IC, G_C were deduced from [small script l], E using (7). We plot [small script l], K_IC, G_C (respectively in panel a, b, c) against E (common horizontal axis) to investigate potential correlations between elastic and fracture properties.

Fig. 7 Linear elastic fracture mechanical properties in 70 independent experiments: (a) elasto-fracture length [small script l] (determined as in Fig. 6b); (b) critical stress intensity factor (or toughness) K_IC = (2π^1/2/3)E[small script l]^1/2; (c) critical strain energy release rate G_C = (π/3)E[small script l]. All are plotted against the Young modulus E. The grey shaded areas represent their distributions (estimates of the probability densities). Error bars in K_IC, G_C were propagated from those in E, [small script l] (respectively and ).

First, we see that E covers a broad range ≈12–28 kPa (see its distribution in grey on top of panel a), despite the fact that we did not intentionally vary it (only p_s was intentionally varied in these experiments). Beads were left to swell in similar conditions for ≈50–150 hours (but swelled little beyond the first ≈24 hours), and had a polymer/water volume fraction ϕ ≈ 0.3–2% (using extreme values of their dry radius 2.5–3 mm and wet radius 11–17 mm). We attribute theses variations in E to inherent variability in the manufacturing of our commercial-grade polyacrylamide gel beads, resulting in slightly different ϕ and relative weight of cross-linkers/monomers.²⁹ This is qualitatively consistent with two previous studies on fractures in polyacrylamide gels at higher ϕ: ref. 19 reported a range E ≈ 1.5–600 kPa sharply increasing with ϕ ≈ 3–9% (see their Fig. 5), while ref. 30 reported E ≈ 50–700 kPa for ϕ ≈ 8–15% and varying cross-linkers/monomers weight fractions ≈2.5–6% (see their Fig. 4).

Second, [small script l] covers a range of ≈0.05–0.2 mm (see its distribution in grey on the right of panel a). This remains an order of magnitude below the minimum crack radius measured in our experiments (R_min ≈ 1 mm), which confirms the validity condition (8) of the penny-shaped model.

Third, the corresponding range for the fracture properties is K_IC ≈ 100–300 Pa m^1/2 and G_C ≈ 1–4 Pa m (or J m⁻²). The latter is at least an order of magnitude larger than the surface tension of water (≈0.07 Pa m) which justifies that we neglected surface tension in the model.

4.3 Growth rate predictions

Pressure measurements. In order to compare the measured radial growth rates α with the predictions from the model (20) we need sufficiently accurate measurements of the super-saturation pressure p_s (which become squared in S²). As mentioned in Section 2.1, these pressure measurements suffered from a loss of CO₂ during the swapping of a ‘seal cap’ for a ‘pressure gauge cap’. We adopted a method to accurately correct for this (explained in the ESI,† Section 1) for a subset of 37 experiments, i.e. about half of those shown in Fig. 7 (the other half being carried prior to this method, in a way that did not allow for subsequent correction). The range of pressure was p_s = 2.31–4.14 bar, corresponding to super-saturation S = 1.07–2.67. Lower or higher p_s were not practical, because cracks were respectively too difficult to initiate or too numerous and fast-growing to measure accurately.

Results. The results of predicted vs. measured growth rates are shown in Fig. 8. Our range of super-saturations S, together with our range of mechanical properties [small script l] (Fig. 7a), allowed us to cover a range of α ≈ 0.01–0.15 mm s⁻¹ (an order of magnitude). We find that the predicted α is on average 27% below the measured value (see the green dashed linear fit ‘y = 0.733x’). However, the discrepancy between predicted and measured α is highly variable (extreme values are ≈±50%).

Fig. 8 Comparison of growth rates: measured from dR/dt (linear fit of the orange data in Fig. 6a) vs. predicted from (20). The fit of the data has slope 0.733 (green dashed line), lower than the expected unit value (black solid line). Errors in the measured α are smaller than the symbol size, while errors in the predicted α (vertical error bars) were propagated from those in S, [small script l] (the relative error is thus ). The error in the measured temperature T (±0.2%), involved in the calculation of D, S, is small enough to be neglected.

Error bars. Although the measured α is virtually exact, the predicted α has a mean relative uncertainty (or ‘error’) of 14% propagated from errors in our measurements of p_s (representing only an average of 1%), and of [small script l] (representing an average of 13%, primarily due to the imperfect V ∝ R^5/2 scaling). These error bounds are such that the measured and predicted α are consistent for a quarter of our data (i.e. error bars cross the 1 : 1 black solid line). However, it is clear that these uncertainties are not sufficient to account for the observed discrepancy. In other words, our model is essentially validated by our measurements, but it appears to systematically under-estimate the growth rate.

Reasons for the under-estimation of α. We believe that this systematic under-estimation of the growth rate by the model can be explained as follows.

First, an imperfect equilibrium between the concentration of dissolved CO₂ in the beads c_s and the gas pressure p_s measured above the solution (due to the surface tension barrier to bubble nucleation explained in the ESI,† Section 1.3) would cause our application of Henry's law c_s = k_hp_s to slightly under-estimate c_s and therefore S (through (12)).

Second, a more quantitatively plausible reason lies in the fact that some cracks do not grow in a perfectly flat plane as we assumed; instead minor inhomogeneities in the gel cause the crack plane to warp slightly. This warped plane, even when barely detectable by the naked eye, can cause our fitting to over-estimate the true crack thickness, its volume, and therefore our measured [small script l] (through (6)).

Third, note that we may not invoke the finite size of the hydrogel bead as a possible reason. It is true that our finite surrounding medium may be felt by the crack, and oppose slightly less resistance than the assumed infinite medium. This would indeed cause us to under-estimate the true crack thickness, and therefore [small script l], and therefore to over-estimate the predicted α; but only if we obtained [small script l] from the mechanical properties E, K_IC. However, in this paper we obtain [small script l] from the observed crack geometry instead (and deduce K_IC from it), which should effectively account for possible finite-size effects of the surrounding medium already. Finally, it is possible that the CO₂ diffusivity D and solubility k_h within the hydrogel are not identical to those in pure water (as assumed, see Table 1). Limited evidence appears to suggest that the polymer may slightly decrease D³¹ and increase k_h.³² Both effects would thus tend to compensate one another, but details remain speculative.

5 Conclusions

Key result. We have shown and explained how initially microscopic (invisible) defects in decompressed hydrogels can, through diffusion and irreversible fracturing, grow as a penny-shaped crack with linear scaling in time R(t) = αt, where α ∝ DS²[small script l]⁻¹ ∝ D(Tk_hp_sE/K_IC)² (in terms of the basic parameters in Table 1), provided that R ≫ max([small script l], [small script l]S⁻²). This key result has at least two applications.

Optical measurements of fracture properties. It allows fracture properties of hydrogels (the toughness K_IC, or equivalently its fracture energy G_C) to be estimated indirectly by optical measurements of dR/dt = α (assuming that all other quantities D, T, k_h, p_s, E can be measured accurately). This new method represents an interesting alternative to that proposed in ref. 30 who used the forced syringe injection of liquid and optical dye attenuation to measure the penny-shaped crack thickness. This new method also complements previous cavitation rheology methods,³³ based on the growth of spherical elastic (non-fracturing) bubbles to estimate material properties such as E.³⁴

Decompression sickness. Our key result also suggests a new potential mechanism for decompression sickness: the growth of inert gas pockets with a potentially catastrophic linear scaling, and a thin penny-shaped geometry that irreversibly fractures body tissues. Future work is undoubtedly needed to assess the relevance of this mechanism in the very complex pathophysiology of decompression sickness.

Open questions. Our work raises at least three open questions.

First, what is the importance of tissue perfusion and of the potential presence of multiple cracks? Both phenomena would limit the available dissolved gas and therefore the crack growth. In particular, we would only expect to find ‘large’ cracks (say of a given radius R₀) in poorly-perfused tissues in which the ‘flushing’ time-scale is ≫R₀/α ∝ R₀D⁻¹S⁻²[small script l].

Second, what is the range of parameters involved in the physiology of decompression sickness? It is apparent that N₂-filled cracks would grow much more slowly than our analogue CO₂-filled cracks due to the 1 : 50 ratio in solubilities k_h (causing α to be reduced by a factor ≈50² = 2500 for a given temperature). We also expect body temperature (at least 10 °C larger than in any of our experiments) to further decrease α (since k_h decays with T faster than 1/T). However, we also need to consider the potentially large (and poorly known) range of body tissue elasticity and toughness, and in particular their ratio, i.e. the elasto-fracture length [small script l] ∝ (K_IC/E)². Tissues with particularly small [small script l] could still be subject to appreciable growth rates.

Third, what is the early-time behaviour of this crack growth? This includes understanding the conditions in which micro-cracks are formed at t = 0. In these experiments we chose mechanical impact; in the body this could be micro-traumata in muscles, tendons, ligaments, joints, bones, or nerves. It also includes investigating the pathway by which initially small, ‘classical’ spherical bubble micro-nuclei might grow large enough (R ≫ [small script l]) to initiate fracture and bifurcate to the penny-shaped crack regime.

This last direction is also relevant to the growth of liquid cavities in polymer matrices mentioned in the introduction.^4,5 Irreversible cavity formation has been reported in ref. 6 and attributed to a ductile response of the material, although, in principle at least, the brittle regime we described could also occur in these systems. An important difference between gas–liquid and liquid–liquid systems is the compressibility of the gas phase and the connection between gas solubility and pressure within the cavity. This may be the reason why liquid–liquid systems show spontaneous cavity growth after quenching, whereas our system requires an impact. We would also expect much slower crack growth in liquid–liquid systems due to the absence of volume change upon phase separation (compared to gas–liquid systems), perhaps making them ideal to study the early stages of crack formation.

Supporting data

The data and codes associated with this paper can be freely downloaded from the repository doi.org/10.17863/CAM.60212.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Mr Zhiyuan Wei for discussions about the modelling. This research was conducted while Y. Z. was an undergraduate student at Trinity College, Cambridge, supported by a Trinity Overseas Bursary. A. L. was supported by an Early Career Fellowship funded by the Leverhulme Trust and Isaac Newton Trust.

Notes and references

1. G. Liger-Belair, G. Polidori and P. Jeandet, Chem. Soc. Rev., 2008, 37, 2490–2511.
2. A. N. Gent and D. A. Tompkins, J. Appl. Phys., 1969, 40, 2520–2525.
3. K. Ittai, L. Vladimir and N. Oded, Bull. Volcanol., 2011, 73, 39–54.
4. R. W. Style, T. Sai, N. Fanelli, M. Ijavi, K. Smith-Mannschott, Q. Xu, L. A. Wilen and E. R. Dufresne, Phys. Rev. X, 2018, 8, 011028.
5. K. A. Rosowski, T. Sai, E. Vidal-Henriquez, D. Zwicker, R. W. Style and E. R. Dufresne, Nat. Phys., 2020, 16, 422–425.
6. J. Y. Kim, Z. Liu, B. M. Weon, T. Cohen, C.-Y. Hui, E. R. Dufresne and R. W. Style, Sci. Adv., 2020, 6, eaaz0418.
7. C. Walsh, PhD thesis, University College London, 2017.
8. V. Papadopoulou, R. J. Eckersley, C. Balestra, T. D. Karapantsios and M.-X. Tang, Adv. Colloid Interface Sci., 2013, 191-192, 22–30.
9. M. Chappell and S. Payne, Respir. Physiol. Neurobiol., 2006, 153, 166–180.
10. O. R. Enríquez, C. Sun, D. Lohse, A. Prosperetti and D. van der Meer, J. Fluid Mech., 2014, 741, R1.
11. H. D. Van Liew and M. P. Hlastala, Respir. Physiol., 1969, 7, 111–121.
12. S. R. Srinivasan, W. A. Gerth and M. R. Powell, Ann. Biomed. Eng., 2002, 30, 232–246.
13. S. A. Mohammadein and K. G. Mohamed, ESAIM: Math. Modell. Numer. Anal., 2016, 21, 762–773.
14. S. Goldman, J. Chem. Phys., 2009, 131, 184502.
15. C. Walsh, E. Stride, U. Cheema and N. Ovenden, J. R. Soc., Interface, 2017, 14, 20170653.
16. D. E. Yount and R. H. Strauss, J. Appl. Phys., 1976, 47, 5081–5089.
17. C.-Y. Lai, Z. Zheng, E. Dressaire and H. A. Stone, Philos. Trans. R. Soc., A, 2016, 374, 20150425.
18. N. J. OKeeffe, H. E. Huppert and P. F. Linden, J. Fluid Mech., 2018, 844, 435–458.
19. S. Kundu and A. J. Crosby, Soft Matter, 2009, 5, 3963.
20. T. Bertrand, J. Peixinho, S. Mukhopadhyay and C. W. MacMinn, Phys. Rev. Appl., 2016, 6, 064010.
21. R. Sander, Atmos. Chem. Phys., 2015, 15, 4399–4981.
22. in CRC Handbook of Chemistry and Physics: a Ready-Reference Book of Chemical and Physical Data, ed. W. M. Haynes, D. R. Lide and T. J. Bruno, CRC Press, 97th edn, 2016.
23. I. N. Sneddon, Philos. Trans. R. Soc., A, 1946, 187, 229–260.
24. T. Kundu, Fundamentals of Fracture Mechanics, CRC Press, 2008.
25. G. R. Irwin, J. Appl. Mech., 1962, 29, 651–654.
26. M. F. Kanninen and C. H. Popelar, Advanced Fracture Mechanics, Oxford University Press, 1985.
27. S. Rattan and A. J. Crosby, ACS Macro Lett., 2019, 8, 492–498.
28. E. Detournay, Annu. Rev. Fluid Mech., 2016, 48, 311–339.
29. A. K. Denisin and B. L. Pruitt, ACS Appl. Mater. Interfaces, 2016, 8, 893–902.
30. N. J. OKeeffe and P. F. Linden, Exp. Mech., 2017, 57, 1483–1493.
31. L. Liu, A. Chakma and X. Feng, J. Membr. Sci., 2008, 310, 66–75.
32. Y. Zhao, Y. Shen, L. Bai and S. Ni, Appl. Surf. Sci., 2012, 261, 708–716.
33. C. W. Barney, C. E. Dougan, K. R. Mcleod, A. Kazemi-Moridani, Y. Zheng, Z. Ye, S. Tiwari, I. Sacligil, R. A. Riggleman, S. Cai, J.-H. Lee, S. R. Peyton, G. N. Tew and A. J. Crosby, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 9157–9165.
34. J. A. Zimberlin, N. Sanabria-DeLong, G. N. Tew and A. J. Crosby, Soft Matter, 2007, 3, 763.

---

Footnotes

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

‡ SCUBA = Self-Contained Underwater Breathing Apparatus.

§ The cap swap was necessary because the available pressure gauge fittings did not provide a sufficiently good seal for >48 hours. Unfortunately this swap caused the system to lose a small amount of CO₂, which we corrected for as described in the ESI,† Section 1.

¶ This impact proved necessary since cracks did not grow spontaneously upon decompression in the range of pressures p_s investigated here. However, preliminary experiments in which the beads were heated well above room temperature proved that thermal shock could also be sufficient to initiate crack growth. This method was not pursued in this work because of the need to carefully control and measure temperature for the validation of our model.

|| This axisymmetry in isotropic materials is expected theoretically: the energy derivation of ref. 25 (eqn (24)) shows that any crack base that is initially non-circular will quickly become circular, even before the crack grows in thickness.

** Here we present the (local) stress intensity factor approach. In linear elasticity, it is equivalent to the (global) energy approach, originally due to Griffith and improved by Irwin, which considers that fracture occurs when the rate of release per unit area of strain energy G stored in the surrounding material exceeds the material's critical value G_C = 3K_IC²/(4E) (under incompressible, axisymmetric i.e. plane strain conditions, see ref. 24 eqn (3.54) and ref. 26 eqn. (3.3)–(3.18)). Note that the surface tension of water is negligible compared to G_C.

†† We use (3) to estimate the maximum excess pressure inside the crack: using our extreme experimental values (see Fig. 7b). At larger radii R(t) ≫ 1 mm, assuming Δp(t) ≪ p_s is very reasonable for the diffusion problem.

‡‡ This creates a radial gradient in the gas concentration c(r,z,t) inside the gel, which we neglected (this is validated later in Section 3.3).

---