Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

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

Henry A. Lockwood , Molly H. Agar and Suzanne M. Fielding
Department of Physics, Durham University, Science Laboratories, South Road, Durham DH1 3LE, UK

Received 28th November 2023 , Accepted 14th February 2024

First published on 22nd February 2024


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.

1 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–13 hydrogels,14 biopolymer gels,15 fibre composites16 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, γt1−α with α < 1, with the shear rate decreasing as [small gamma, Greek, dot above]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 [small gamma, Greek, dot above]min at time τmin. The shear rate then increased dramatically as [small gamma, Greek, dot above] ∼ (τft)−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

image file: d3sm01608k-t1.tif(1)

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 τminvia 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 [small gamma, Greek, dot above]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 (τft)−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 experiments18 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.

2 Model

Our approach will be to model a macroscopic sample of gel by considering the underlying physics at the mesoscopic 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 with combining dot above] = [small gamma, Greek, dot above], giving elastic response. To model the strain-induced breakage of gel strands, when any element's elastic energy image file: d3sm01608k-t2.tif exceeds a local yielding threshold Ey, it fails plastically with stochastic rate τ0−1. Activated strand breakage is also captured via a yielding rate image file: d3sm01608k-t3.tif for elements even below threshold. In materials for which the energy barriers to strand breakage are comparable with thermal energies, we take x = kBT. 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” −tw < t < 0 prior to the switch-on 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/xg) in which xg is a model parameter. As an initial condition at t = −tw, we assume a distribution of thresholds ρ(E). For x < xg, 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 [small gamma, Greek, dot above](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 Ly, with M elastoplastic elements on each streamline.26 Imposing force balance then allows us to calculate the shear rate profile [small gamma, Greek, dot above](y, t), which becomes highly heterogeneous as the sample yields. We do this via two different algorithms. The first considers a solvent stress of small viscosity η = 0.05, additive to the elastoplastic stress,27 giving a total stress Σ = σ + η[small gamma, Greek, dot above]. The results shown are converged to the limit η → 0, suited to the relevant experimental regime, η0. The second has η = 0.26 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: twtw[1 + ε[thin space (1/6-em)]cos(2πy/Ly)]. 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 Ly = 1. We also set xg = 1, thereby setting typical local yield strains of order unity. Unless otherwise stated we use S = 10 streamlines with M = 100[thin space (1/6-em)]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 tw, and the noise temperature x.

3 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 [small gamma, Greek, dot above](t) in Fig. 1(a) and (b). The basic physics seen experimentally5 is clearly evident: primary creep with increasing strain γγ0t1−α and decreasing strain rate [small gamma, Greek, dot above]tα is followed by a sudden near-divergence of these quantities, signifying catastrophic failure. (True divergence is averted in Fig. 1via the solvent viscosity, which ensures a finite [small gamma, Greek, dot above] = Σ0/η even after all gel strands break.) The creep exponent α = x for (noise) temperatures x < 1. Experiments find α = 0.85 in casein gels5 and α = 0.7 in milk gels.18
image file: d3sm01608k-f1.tif
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. (a) 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. (b) Corresponding strain rate [small gamma, Greek, dot above](t). Dotted line shows the power −0.3 = −x. Crosses indicate the minimum value of the strain rate [small gamma, Greek, dot above]min in each simulation, and the time τmin at which it occurs. (c) 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 = 105. η = 0.05.

Alongside the global strain signal γ(t), we also plot in Fig. 1(c) the degree to which the strain field becomes heterogeneous, i.e., shear banded, as a function of the flow-gradient 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. 1(a). (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 tw at a fixed noise temperature x. Collapse of the data when normalised as τf/tw indicates a scaling τftw, 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–32 The log–lin axes in (b) reveal τf ∼ exp(−0) at larger Σ0.

image file: d3sm01608k-f2.tif
Fig. 2 Failure time τfvs. imposed stress Σ0 at a noise temperature x = 0.3. (a) On a log–log scale, showing a power law dependence τfΣβ0 at low Σ0. Dashed line shows power β = 3.3. Inset: Exponent β vs. x. (b) On a log–lin scale, showing an exponential dependence τf ∼ exp(−0) at larger Σ0. Dashed line shows exponent m = 3.0. Inset: Exponent m vs. x. Waiting time tw = 10n with n = 3.0, 4.0, 5.0, 6.0 for the red circles, green triangles, blue squares and cyan pentagons respectively. η = 0.05.

A notable achievement of the experiments in ref. 5 was to fit the shear rate across all three time regimes – primary, 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 [small gamma, Greek, dot above](t) in three different ways in Fig. 3. Panel (a) shows the normalised strain rate [small gamma, Greek, dot above]/[small gamma, Greek, dot above]minvs. normalised time t/τf on log–log axes, to emphasise the primary creep regime. The dotted line shows the power [small gamma, Greek, dot above]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 [small gamma, Greek, dot above]/[small gamma, Greek, dot above]minvs. normalised time to failure, (τft)/τf. In this figure, for consistency with the counterpart figure in the experimental ref. 5, we have inverted the time axis so as to progress towards failure from left to right. In this way, the limit τft → 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.

image file: d3sm01608k-f3.tif
Fig. 3 Strain rate [small gamma, Greek, dot above] normalised by its minimum value [small gamma, Greek, dot above]min. (a) Plotted versus time t normalised by the failure time τf, with log–log axes to emphasise the primary creep regime, [small gamma, Greek, dot above]tα. Dotted line shows the power α = x = 0.5. (b) Plotted likewise, with lin–lin axes to emphasise the secondary crossover regime. (c) Plotted as a function of time remaining to failure τft, 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 = 109. M = 10[thin space (1/6-em)]000. η = 0.0.

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, tw and Σ0. The dashed line shows the linear relationship τf = τmin/c with the same prefactor c = 0.56 reported experimentally.

image file: d3sm01608k-f4.tif
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 = 10n 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.

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 linear, corresponding to a near-uniform strain profile γ(y) = ∂yU(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.

image file: d3sm01608k-f5.tif
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 = 105. η = 0.05.

4 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 Basquin-like 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 models33 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 work37 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 rearrangements16 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.


This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 885146). We thank Andrew Clarke for discussions and SLB (Schlumberger Cambridge Research Ltd) for support.

Notes and references

  1. W. C. Poon, L. Starrs, S. Meeker, A. Moussaid, R. M. Evans, P. Pusey and M. Robins, Faraday Discuss., 1999, 112, 143–154 RSC.
  2. P. J. Skrzeszewska, J. Sprakel, F. A. de Wolf, R. Fokkink, M. A. Cohen Stuart and J. van der Gucht, Macromolecules, 2010, 43, 3542–3548 CrossRef CAS.
  3. D. Bonn, H. Kellay, M. Prochnow, K. Ben-Djemiaa and J. Meunier, Science, 1998, 280, 265–267 CrossRef PubMed.
  4. F. Gobeaux, E. Belamie, G. Mosser, P. Davidson and S. Asnacios, Soft Matter, 2010, 6, 3769–3777 RSC.
  5. M. Leocmach, C. Perge, T. Divoux and S. Manneville, Phys. Rev. Lett., 2014, 113, 038303 CrossRef CAS PubMed.
  6. V. Gopalakrishnan and C. Zukoski, J. Rheol., 2007, 51, 623–644 CrossRef CAS.
  7. T. Gibaud, D. Frelat and S. Manneville, Soft Matter, 2010, 6, 3482–3488 RSC.
  8. V. Grenard, T. Divoux, N. Taberlet and S. Manneville, Soft Matter, 2014, 10, 1555–1571 RSC.
  9. J. Sprakel, S. B. Lindström, T. E. Kodger and D. A. Weitz, Phys. Rev. Lett., 2011, 106, 248303 CrossRef.
  10. J. H. Cho and I. Bischofberger, Soft Matter, 2022, 18, 7612–7620 RSC.
  11. S. Aime, L. Ramos and L. Cipelletti, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 3587–3592 CrossRef CAS.
  12. E. Moghimi, A. B. Schofield and G. Petekidis, J. Phys.: Condens. Matter, 2021, 33, 284002 CrossRef CAS.
  13. B. J. Landrum, W. B. Russel and R. N. Zia, J. Rheol., 2016, 60, 783–807 CrossRef CAS.
  14. S. N. Karobi, T. L. Sun, T. Kurokawa, F. Luo, T. Nakajima, T. Nonoyama and J. P. Gong, Macromolecules, 2016, 49, 5630–5636 CrossRef CAS.
  15. A. Pommella, L. Cipelletti and L. Ramos, Phys. Rev. Lett., 2020, 125, 268006 CrossRef CAS PubMed.
  16. H. Nechad, A. Helmstetter, R. El Guerjouma and D. Sornette, Phys. Rev. Lett., 2005, 94, 045501 CrossRef CAS PubMed.
  17. J. Rosti, J. Koivisto, L. Laurson and M. J. Alava, Phys. Rev. Lett., 2010, 105, 100601 CrossRef.
  18. J. Bauland, M. Leocmach, M.-H. Famelart and T. Croguennec, Soft Matter, 2023, 19, 3562–3569 RSC.
  19. E. N. D. C. Andrade, Proc. R. Soc. Lond. Ser. A-Contain. Pap. Math. Phys. Character, 1910, 84, 1–12 CAS.
  20. M.-C. Miguel, A. Vespignani, M. Zaiser and S. Zapperi, Phys. Rev. Lett., 2002, 89, 165501 CrossRef.
  21. O. H. Basquin, Proc., Am. Soc. Test. Mater., 1910, 625–630 Search PubMed.
  22. F. C. Monkman, Proc., Am. Soc. Test. Mater., 1956, 593–620 Search PubMed.
  23. T. Gallot, C. Perge, V. Grenard, M.-A. Fardin, N. Taberlet and S. Manneville, Rev. Sci. Instrum., 2013, 84, 045107 CrossRef.
  24. A. Nicolas, E. E. Ferrero, K. Martens and J.-L. Barrat, Rev. Mod. Phys., 2018, 90, 045006 CrossRef CAS.
  25. P. Sollich, F. Lequeux, P. Hébraud and M. E. Cates, Phys. Rev. Lett., 1997, 78, 2020 CrossRef CAS.
  26. S. Fielding, M. Cates and P. Sollich, Soft Matter, 2009, 5, 2378–2382 RSC.
  27. H. J. Barlow, J. O. Cochran and S. M. Fielding, Phys. Rev. Lett., 2020, 125, 168003 CrossRef CAS PubMed.
  28. R. L. Moorcroft and S. M. Fielding, J. Rheol., 2014, 58, 103–147 CrossRef CAS.
  29. T. Sentjabrskaja, P. Chaudhuri, M. Hermes, W. Poon, J. Horbach, S. Egelhaaf and M. Laurati, Sci. Rep., 2015, 5, 11884 CrossRef CAS PubMed.
  30. P. Ballesta and G. Petekidis, Phys. Rev. E, 2016, 93, 042613 CrossRef CAS.
  31. M. Siebenbürger, M. Ballauff and T. Voigtmann, Phys. Rev. Lett., 2012, 108, 255701 CrossRef.
  32. T. Divoux, C. Barentin and S. Manneville, Soft Matter, 2011, 7, 8409–8418 RSC.
  33. S. Pradhan, A. Hansen and B. K. Chakrabarti, Rev. Mod. Phys., 2010, 82, 499 CrossRef.
  34. E. Jagla, Phys. Rev. E, 2011, 83, 046119 CrossRef CAS PubMed.
  35. F. Kun, H. Carmona, J. S. Andrade Jr and H. J. Herrmann, Phys. Rev. Lett., 2008, 100, 094301 CrossRef CAS PubMed.
  36. Z. Halász, Z. Danku and F. Kun, Phys. Rev. E, 2012, 85, 016116 CrossRef PubMed.
  37. R. L. Moorcroft and S. M. Fielding, Phys. Rev. Lett., 2013, 110, 086001 CrossRef PubMed.

This journal is © The Royal Society of Chemistry 2024