Power law creep and delayed failure of gels and fibrous materials under stress

Motivated by recent experiments studying the creep and breakup of a protein gel under stress, we introduce a simple mesoscopic model for the irreversible failure of gels and fibrous materials, and demonstrate it to capture much of the phenomenology seen experimentally. This includes a primary creep regime in which the shear rate decreases as a power law over several decades of time, a secondary crossover regime in which the shear rate attains a minimum, and a tertiary regime in which the shear rate increases dramatically up to a finite time singularity, signifying irreversible material failure. The model also captures a linear Monkman-Grant scaling of the failure time with the earlier time at which the shear rate attained its minimum, and a Basquin-like power law scaling of the failure time with imposed stress, as seen experimentally. The model furthermore predicts a slow accumulation of low levels of material damage during primary creep, followed by the growth of fractures leading to sudden material failure, as seen experimentally.


I. INTRODUCTION
Materials subject to loading sustained over a lengthy duration will often catastrophically fail after a long induction time, without apparent warning.This has important implications in myriad contexts, including the shelf life of products, the integrity of construction materials, and the prediction of geophysical phenomena such as avalanches, mudslides and earthquakes.In the context of soft gels and fibrous network materials, delayed failure under constant load has been observed in transient gels [1,2], polymeric gels [3], collagen networks [4], protein gels [5], colloidal gels [6][7][8][9][10][11][12][13], hydrogels [14], biopolymer gels [15], fibre composites [16] and paper [17].In some cases yielding is reversible, with the material reforming and recovering elasticity once the load is removed.In others, failure is irreversible and the original network is permanently broken.
In this work, we focus on the irreversible failure of gels and fibrous materials, motivated in particular by recent experiments demonstrating that the behaviour of a protein gel under sustained loading is strikingly similar to that of brittle solids [5,18].In these experiments, an initial regime of slow creep over several decades of time t later gave way to catastrophic material failure at a time τ f .In the primary creep regime, the shear strain γ increased up to O(1) in a sublinear way, γ ∼ t 1−α with α < 1, with the shear rate decreasing as γ ∼ t −α , reminiscent of Andrade creep in solids [19,20].In a secondary crossover regime for times 0.1 ≲ t/τ f ≲ 0.9, the shear rate attained a minimum γmin at time τ min .The shear rate then increased dramatically as γ ∼ (τ f − t) −1 in a final tertiary regime, diverging in a finite time singularity at τ f , signifying catastrophic material failure.Across all three regimes, the shear rate fit the single master relation The same experiments demonstrated the failure time to decrease with imposed stress as τ f ∼ Σ −β 0 , reminiscent of the Basquin law of fatigue in solids [21], and to scale linearly with τ min via a Monkman-Grant relation [22], τ f = τ min /c with c < 1. Importantly, this suggests that sudden failure at time τ f might be forecast once a shear rate minimum has earlier been observed at τ min .
By combining rheology with optical and ultrasonic imaging [23], these experiments also tracked the strain field within the gel, in tandem with the global strain γ.In the primary creep regime, no macroscopic strain localisation was seen, although the authors noted that this does not rule out local rearrangements on scales below those imaged.Indeed, separate recovery tests suggested low levels of irreversible material damage during primary creep.Fractures then developed during the secondary regime, and grew to cause failure in the tertiary regime.
Motivated by these experiments, we introduce a simple mesoscopic model for the irreversible breakup of gels.Simulations in the step-stress protocol show it to capture the experimental phenomenology just described: (i) primary creep with decreasing shear rate γ ∼ t −α , (ii) a secondary cross-over regime in which the shear rate attains a minimum at time τ min , (iii) a tertiary regime in which the shear rate increases dramatically as (τ f − t) −1 , leading to (iv) catastrophic failure in a finite time singularity at time τ f , (v) a linear Monkman-Grant scaling of τ f with τ min , (vi) a Basquin-like scaling of τ f with Σ 0 , (vii) accumulating low levels of damage during primary creep, followed by (viii) the growth of macroscopic fractures leading to catastrophic material failure.
Although we shall mostly compare our results with the particular experiments of Ref. [5], other experiments [18] report the same phenomenology, with quantitative values of the exponents α, β however being non-universal across materials.Therefore, our focus will be on robustly capturing all the qualitative behaviour just described, while also characterising how the values of exponents such as α, β, etc. depend on model parameters.

II. MODEL
Our approach will be to model a macroscopic sample of gel by considering the underlying physics at the arXiv:2311.16778v2[cond-mat.soft]15 Feb 2024 mescoscopic level of a few gel strands and their associated crosslinks.In the spirit of existing elastoplastic models [24], each local collection of strands and crosslinks is represented by a single elastoplastic element, with local shear strain l and stress kl (discarding normal stresses).The macroscopic elastoplastic stress σ is the average of these local elemental ones.Before any local failure event, each element affinely follows the macroscopic shear strain, l = γ, giving elastic response.To model the strain-induced breakage of gel strands, when any element's elastic energy E = 1 2 kl 2 exceeds a local yielding threshold E y , it fails plastically with stochastic rate τ −1 0 .Activated strand breakage is also captured via a yielding rate τ −1 0 exp ( 1 2 kl 2 − E y )/x for elements even below threshold.In materials for which the energy barriers to strand breakage are comparable with thermal energies, we take x = k B T .For much larger barriers, x is taken to be a mean field temperature modeling strand breakage due to the propagation of stresses from breakages elsewhere.On the timescales of creep and yielding, the breakage of any strand is assumed irreversible, with an infinite reformation time in comparison, τ ref = ∞.The associated element carries no stress thereafter.
To model the much slower process of gelation during a "waiting time" −t w < t < 0 prior to the switchon of stress at t = 0, we must of course model the (re)formation of gel strands.During gelation, therefore, we include an additional model ingredient: that after any element locally yields, it reforms with local strain l = 0 and a new energy threshold drawn from a distribution ρ(E) ∼ exp(−E/x g ) in which the model parameter x g is the glass transition temperature.As an initial condition at t = −t w , we assume a distribution of thresholds ρ(E).For x < x g , the model captures ageing, with dynamics that become slower over time as the system gels.We take element reformation to be infinitely fast compared with this slow process of gelation, with τ ref = 0.This is oversimplified: a fuller description might consider a finite reformation time τ ref consistently during gelation, creep and yielding.However, our focus here is not to model gelation itself, but the subsequent creep and yielding of a gel, however formed, under stress.
During gelation, with immediate element reformation after breakage via an assumed τ ref = 0, compared with the much slower timescale of gelation, the model is the same as the soft glassy rheology model [25].During creep, the model differs from the SGR model in having irreversible element breakage via an assumed τ ref = ∞, compared with the much faster timescale of yielding.
So far, we have described the model as though the strain rate γ(t) were the same for all elements.In fact, to allow for strain localisation during material failure, we consider an array of S streamlines stacked across the flow gradient direction y of a shear cell of gap width L y , with M elastoplastic elements on each streamline [26].Imposing force balance then allows us to calculate the shear rate profile γ(y, t), which becomes highly heterogeneous as the sample yields.We do this via two dif- ferent algorithms.The first considers a solvent stress of small viscosity η = 0.05, additive to the elastoplastic stress [27], giving a total stress Σ = σ + η γ.The second has η = 0 [26].The results shown are converged to the limit η → 0, suited to the relevant experimental regime, η ≪ kτ 0 .The η = 0 algorithm was needed only to obtain the scaling of shear rate with time in the final tertiary regime, as the gel finally breaks completely.To seed the formation of strain bands, we initialise our creep simulations with a small heterogeneity: t w → t w [1 + ϵ cos(2πy/L y )].Indeed, in physical practice, most rheometers seed slight heterogeneity due to device curvature.For a fuller discussion of other possible sources of heterogeneity that might seed band formation, see Ref. [28].
We choose units in which the attempt time τ 0 = 1, local modulus k = 1 and gap width L y = 1.We also set x g = 1, thereby setting typical local yield strains of order unity.Unless otherwise stated we use S = 10 streamlines with M = 100, 000 elements per streamline, a perturbation amplitude ϵ = 0.1, and numerical time-step dt = 0.01.Parameters explored are then the imposed stress Σ 0 , the waiting time t w , and the noise temperature x.

III. RESULTS
Following the switch-on of a stress of magnitude Σ 0 at time t = 0, we report the global shear strain γ(t) (relative to the elastic strain γ 0 = Σ 0 that arises immediately upon stress switch-on) and shear rate γ(t) in Fig. 1a,b).The basic physics seen experimentally [5] is clearly evident: primary creep with increasing strain γ − γ 0 ∼ t 1−α and decreasing strain rate γ ∼ t −α is followed by a sudden near-divergence of these quantities, signifying catastrophic failure.(True divergence is averted in Fig. 1 via the solvent viscosity, which ensures a finite γ = Σ 0 /η even after all gel strands break.)The creep exponent α = x for (noise) temperatures x < 1. Experiments find α = 0.85 in casein gels [5] and α = 0.7 in milk gels [18].
Alongside the global strain signal γ(t), we also plot in Fig. 1c) the degree to which the strain field becomes heterogeneous, i.e., shear banded, as a function of the flowgradient coordinate y across the shear cell.Specifically, we characterise this degree of strain localisation, i.e., of shear banding, by defining the quantity ∆γ(t).This is the standard deviation in strain across the streamlines stacked across y at any time t, normalised by the globally average strain γ(t).This remains small during primary creep, consistent with the low levels of irreversible damage evidenced via strain recovery experiments [5].The level of strain localisation then increases dramatically as the material yields, as seen experimentally [5].
As indicated by the dashed line, we define the failure time τ f as the time at which the global strain attains a value equal to unity in Fig. 1a).(In any simulations with η = 0, we instead define τ f as the time at which the strain rate actually diverges.We have checked that these definitions give the same scaling of τ f with Σ 0 and τ min , with only a small fractional quantitative difference.)This is plotted vs the imposed stress Σ 0 in Fig. 2, for several values of the waiting time t w at a fixed noise temperature x.Collapse of the data when normalised as τ f /t w indicates a scaling τ f ∼ t w , consistent with the aging predicted by our model during gelation.The log-log axes of panel a) reveal the Basquin-like scaling τ f ∼ Σ −β 0 at low Σ 0 , as seen experimentally with β = 5.5 in Ref. [5] and β = 2.0 in Ref. [18].This scaling predicts that any stress, however small, will always eventually cause failure.This has important implications for the shelf life of gel-based products (which are subject to gravity at least), and contrasts notably with the physics of soft glassy materials such as dense emulsions, colloids and microgels in which the timescale of yielding diverges at a non-zero yield stress Σ y , with indefinite creep for Σ 0 < Σ y [29][30][31][32].The log-lin axes in b) reveal τ f ∼ exp(−mΣ 0 ) at larger Σ 0 .
A notable achievement of the experiments in Ref. [5] was to fit the shear rate across all three time regimesprimary, secondary and tertiary -to the single master scaling form of Eqn. 1.To explore whether such a fit holds here, we plot the strain rate γ(t) in three different ways in Fig. 3. Panel a) shows the normalised strain rate γ/ γmin vs normalised time t/τ f on log-log axes, to emphasise the primary creep regime.The dotted line shows the power γ ∼ t −α , with α = x.Panel b) shows the same quantity, now on lin-lin axes to emphasise the secondary crossover regime.As seen experimentally, this secondary regime is seen for times 0.1 ≲ t/τ f ≲ 0.9, with the shear rate attaining a minimum during it at t = τ min .Finally, panel c) plots the normalised shear rate γ/ γmin vs normalised time to failure, (τ f − t)/τ f .In this figure, for consistency with the counterpart figure in the experimental reference [5], we have inverted the time axis so as to progress towards failure from left to right.In this way, the limit τ f − t → 0 at the right hand side of the figure corresponds to the limit in which t → τ f and the sample finally catastrophically fails.The +1.0 power seen experimentally is indicated by the dotted line.The dashed line in all three panels shows the single experimental master scaling of Eqn. 1.This can be seen to perform well in the primary and tertiary regimes, but significantly less well in the secondary crossover regime.
Another notable achievement of Ref. [5] was to establish a linear Monkman-Grant relationship between the time τ min at which the shear rate attains its minimum in the secondary regime, with the later time τ f of catastrophic material failure in the tertiary regime.Such a relation suggests that the approaching time of material failure might be forecast in a given experiment, once the earlier shear rate minimum has been observed.Fig. 4 shows the same relation to hold across simulations performed for a wide range of x, t w and Σ 0 .The dashed line shows the linear relationship τ f = τ min /c with the same prefactor c = 0.56 reported experimentally.
Finally in Fig. 5 we plot the displacement profile U (y) at several times during creep and yielding.During early creep, the displacement profile remains essentially lin-ear, corresponding to a near-uniform strain strain profile γ(y) = ∂ y U (y), consistent with the lack of strain localisation in the imaging experiments in this regime [5].In contrast, as the material later yields, strong strain localisation arises, as seen experimentally.

IV. CONCLUSIONS
To summarise, we have introduced a mesoscopic model for the irreversible breakup of gels.Its predictions under stress capture the experimentally observed phenomenology, including a primary regime of power-law creep, a secondary cross-over regime in which the shear rate attains a minimum, and a tertiary regime in which the shear rate diverges in a finite time singularity (for η → 0), signifying material failure.Our results also confirm the Basquinlike and Monkman-Grant scalings of failure time seen experimentally.The model furthermore predicts low levels of material damage during primary creep, later giving way to macroscopic strain localisation and material failure, as seen experimentally.In contrast, theories based on simpler so-called fibre bundle models [33] predict Andrade creep and failure but only above a critical load [34]; or Basquin's law but without creep [35,36].
Despite these successes, our approach has several shortcomings.First, we assumed for simplicity a fast gel strand (re)formation time τ ref = 0.0 compared to the slow timescale of initial gelation, but an infinite τ ref compared to the faster timescale of yielding under stress.Future work might consider a finite τ ref consistently during gelation and yielding.Arguably, however, our approach does not accurately model the process of gelation itself: the dynamics for t < 0 might better be seen as leading to a sensible initial condition for the creep simulations.As noted above, setting τ ref = 0 throughout, including during straining, recovers the original SGR model.This has been used in earlier work [37] to model creep and yielding of soft glassy systems, including gels, that flow after yielding.Second, the experiments of Ref. [5] suggested part of the primary power law creep to arise from reversible viscoelastic squeezing of fluid through the gel matrix.In contrast, in our model all creep arises via local plasticity.Indeed, the origin of power law creep in amorphous materials is controversial, having been variously attributed to plastic rearrangements [16] and linear viscoelasticity [4,5].Finally, we have modelled spatial heterogeneity only in the gradient direction y, whereas significant heterogeneity was also observed experimentally in the vorticity direction z.Future modelling efforts should increase the dimensionality simulated.

CONFLICTS OF INTEREST
There are no conflicts to declare.

FIG. 1 .
FIG. 1. Creep and failure under a shear stress imposed for times t > 0.0.Stress Σ0 = 1.0, 1.2 • • • 2.0 in curves upwards.Top: Strain γ − γ0 as a function of time t since the stress was imposed.Dashed line denotes the strain value γ = 1.0 used to define the time τ f at which material failure occurs.Dotted line shows the power 0.7 = 1 − x.Recall that γ0 is the elastic strain that arises immediately upon stress switch-on.Middle: corresponding strain rate γ(t).Dotted line shows the power −0.3 = −x.Crosses indicate the minimum value of the strain rate γmin in each simulation, and the time τmin at which it occurs.Bottom: corresponding growing strain heterogeneity (degree of shear banding), defined at any time t by the standard deviation of strain across y, normalised by the global strain γ(t).Dotted line shows the power 0.8.Noise temperature x = 0.3.Waiting time tw = 10 5 .η = 0.05.

FIG. 3 .
FIG.3.Strain rate γ normalised by its minimum value γmin.Top: Plotted versus time t normalised by the failure time τ f , with log-log axes to emphasise the primary creep regime, γ ∼ t −α .Dotted line shows the power α = x = 0.5.Middle: Plotted likewise, with lin-lin axes to emphasise the secondary crossover regime.Bottom: plotted as a function of time remaining to failure τ f − t, normalised by τ f , to emphasise the tertiary regime of final failure.Note that the time axis is visually inverted so as to move towards failure from left to right.Dotted line shows the power +1.0.Dashed line in each panel shows the fit to the master scaling form of Eqn. 1 with µ = 0.019 and λ = 0.26.Imposed stress Σ0 = 0.1, 0.2 • • • 1.0 in curves black, red, green, blue, yellow, brown, grey, violet, cyan and magenta.Noise temperature x = 0.5, waiting time tw = 10 9 .M = 10000.η = 0.0.

FIG. 4 .
FIG. 4. Time τ f of final material failure versus the earlier time τmin at which the shear rate attained its minimum.Each data point corresponds to a single simulation.Data for noise temperatures x = 0.30, 0.35 • • • 0.50 are shown by circles, triangles, squares, pentagons, stars and crosses respectively.Data for waiting times tw = 10 n with n = 3.0, 4.0, 5.0, 6.0 are shown by red, green, blue and cyan symbols respectively.Imposed stress Σ0 = 0.05, 0.10, 0.15 • • • 2.0 increases right to left for each symbol shape and colour.Dashed line shows the linear Monkman-Grant relationship observed experimentally: τ f = τmin/c with c = 0.56.η = 0.05.

FIG. 5 .
FIG. 5. Displacement profile U (y) as a function of position across the gradient direction y, normalised by the displacement at y 1.0.Profiles are shown for times t indicated by symbols in the creep curve γ(t) of the inset.The colour of each displacement curve indicates the corresponding time in the creep curve.Imposed stress Σ0 = 1.6.Noise temperature x = 0.3.Waiting time tw = 10 5 .η = 0.05.