Concurrent factors determine toughening in the hydraulic fracture of poroelastic composites

TL;DR

This paper investigates how hydraulic coupling in poroelastic composites can enhance toughness by promoting distributed cracking, with implications for designing flaw-insensitive, tougher materials.

Contribution

It extends a previous model to analyze how material and loading parameters influence toughness through crack distribution mechanisms.

Findings

Rapid loading promotes distributed cracking and increases toughness.

Hydraulic coupling enables flaw-insensitive toughening in composites.

Theoretical results suggest new material architectures for higher toughness.

Abstract

Brittle materials fail catastrophically. In consequence of their limited flaw-tolerance, failure occurs by localized fracture and is typically a dynamic process. Recently, experiments on epithelial cell monolayers have revealed that this scenario can be significantly modified when the material susceptible to cracking is adhered to a hydrogel substrate. Thanks to the hydraulic coupling between the brittle layer and the poroelastic substrate, such a composite can develop a toughening mechanism that relies on the simultaneous growth of multiple cracks. Here, we study this remarkable behaviour by means of a detailed model, and explore how the material and loading parameters concur in determining the macroscopic toughness of the system. By extending a previous study, our results show that rapid loading conveys material toughness by promoting distributed cracking. Moreover, our theoretical findings may suggest innovative architectures of flaw-insensitive materials with higher toughness.

Tables

Table 1. Table 1: Values of the dimensionless groups used in the numerical simulations.

Parameter	Value
$h/a_o$	$10$
$h/s$	$5/3$
$h/H$	$0.04$
$G_e/\Lambda$	$0.02$
${\Gamma }_R/{\Gamma }_L$	$0.9$
$\chi$	$0.46$
${\mu }_o/ℛT$	$0$
$G_e{a}_o/{\Gamma }_L$	$6.7$
${\Gamma }_L{a}_o/\eta D$	$2.5$

Equations

$$\delta (X_2,t)=[[f_1(0,X_2,t)]]=2f_1^{\mathrm{c}}(X_2,t),$$

$${\mathbf{v}}_{\mathrm{c}}^{\pm }(x_2,t)=\pm \frac{1}{2}{\dot{\delta }}_{\mathrm{s}}(x_2,t){\mathbf{e}}_1+w(x_2,t){\mathbf{e}}_2,$$

$$\begin{array}{rl}{\dot{\delta }}_{\mathrm{s}}(x_2,t)& =\dot{\delta }(X_2,t)\circ (f_2^{\mathrm{c}}{)}^{-1}(x_2,t)=\\ & =\frac{\partial {\delta }_{\mathrm{s}}(x_2,t)}{\partial t}+w(x_2,t)\frac{\partial {\delta }_{\mathrm{s}}(x_2,t)}{\partial x_2}.\end{array}$$

$${\int }_{\partial {\mathcal{V}}_t}\mathbf{v}\cdot \mathbf{n}\text{mbox}da=0,$$

$${\mathbf{n}}_{\mathrm{c}}^{\pm }=\frac{1}{\Vert {\mathbf{n}}_{\mathrm{c}}\Vert }\left(\pm {\mathbf{e}}_1-\frac{1}{2}\frac{\partial {\delta }_{\mathrm{s}}}{\partial x_2}{\mathbf{e}}_2\right),$$

$${\int }_{{\bar{x}}_2}^{{\bar{x}}_2^{\prime }}\left({\dot{\delta }}_{\mathrm{s}}-w\frac{\partial {\delta }_{\mathrm{s}}}{\partial x_2}\right)\text{mbox}d{x}_2,$$

$$Q({\bar{x}}_2^{\prime },t)-Q({\bar{x}}_2,t)={\int }_{{\bar{x}}_2}^{{\bar{x}}_2^{\prime }}\frac{\partial Q}{\partial x_2}\text{mbox}d{x}_2,$$

$$Q(x_2,t)={\int }_{-{\delta }_{\mathrm{s}}(x_2,t)/2}^{{\delta }_{\mathrm{s}}(x_2,t)/2}{v}_2(x_1,x_2,t)\text{mbox}d{x}_1$$

$${\int }_{{\bar{x}}_2}^{{\bar{x}}_2^{\prime }}\left({\dot{\delta }}_{\mathrm{s}}-w\frac{\partial {\delta }_{\mathrm{s}}}{\partial x_2}+\frac{\partial Q}{\partial x_2}\right)\text{mbox}d{x}_2=0.$$

$$\frac{\partial w(x_2,t)}{\partial x_2}=\left(\frac{1}{\lambda (X_2,t)}\frac{\partial \lambda (X_2,t)}{\partial t}\right)\circ (f_2^{\mathrm{c}}{)}^{-1}(x_2,t),$$

$$\frac{\partial {\delta }_{\mathrm{s}}(x_2,t)}{\partial x_2}=\left(\frac{1}{\lambda (X_2,t)}\frac{\partial \delta (X_2,t)}{\partial X_2}\right)\circ (f_2^{\mathrm{c}}{)}^{-1}(x_2,t).$$

$$\frac{\partial \left(\delta \lambda \right)}{\partial t}+\frac{\partial q}{\partial X_2}=0,$$

$$q_{\mathrm{s}}(x_2,t)=Q(x_2,t)-w(x_2,t){\delta }_{\mathrm{s}}(x_2,t)$$

$${\mathbf{e}}_2\cdot {\int }_{\partial {\mathcal{V}}_t}{\mathbf{T}}_{\mathrm{f}}\mathbf{n}\text{mbox}da=0.$$

$${\mathbf{T}}_{\mathrm{f}}=-p_{\mathrm{s}}\mathbf{I}+2\eta \mathrm{sym}(\text{mbox}\mathrm{grad}\mathbf{v}),$$

$$v_2(x_1,x_2,t)=4(w-v_2^{\max }){\left(\frac{x_1}{{\delta }_{\mathrm{s}}}\right)}^2+v_2^{\max },$$

$$\begin{array}{rl}{\mathbf{t}}_{\mathrm{f}}^{\pm }& ={\mathbf{T}}_{\mathrm{f}}{\mathbf{n}}_{\mathrm{c}}^{\pm }\approx \\ & \approx -p_{\mathrm{s}}{\mathbf{n}}_c^{\pm }-6\eta \frac{q_{\mathrm{s}}}{{\delta }_{\mathrm{s}}^2}({\mathbf{e}}_1\otimes {\mathbf{e}}_2+{\mathbf{e}}_2\otimes {\mathbf{e}}_1){\mathbf{n}}_{\mathrm{c}}^{\pm },\end{array}$$

$$-12{\int }_{{\bar{x}}_2}^{{\bar{x}}_2^{\prime }}\frac{q_{\mathrm{s}}}{{\delta }_{\mathrm{s}}^2}\text{mbox}d{x}_2+{\int }_{{\bar{x}}_2}^{{\bar{x}}_2^{\prime }}{p}_{\mathrm{s}}\frac{\partial {\delta }_{\mathrm{s}}}{\partial x_2}\text{mbox}d{x}_2.$$

$$\displaystyle-p_{\rm s}(x_{2},t)\delta_{\rm s}(x_{2},t)\Big{|}^{\bar{x}^{\prime}_{2}}_{\bar{x}_{2}}=-\int_{\bar{x}_{2}}^{\bar{x}^{\prime}_{2}}{\frac{\partial(p_{\rm s}\delta_{\rm s})}{\partial x_{2}}\,\mbox{d}x_{2}}\,.$$

$$q=-\frac{{\delta }^3}{12\eta \lambda }\frac{\partial p_{\mathrm{f}}}{\partial X_2},$$

$${\mathbf{s}}_{\mathrm{c}}^{\pm }=\mp \sigma (\delta ){\mathbf{e}}_1.$$

$$\sigma (\delta )=(1-d)\frac{{\sigma }_{\text{o}}}{{\delta }_{\text{o}}}\delta ,$$

$$d=\frac{{\delta }_{\text{c}}({\delta }^{\text{m}}-{\delta }_{\text{o}})}{{\delta }^{\text{m}}({\delta }_{\text{c}}-{\delta }_{\text{o}})},$$

$${\mathbf{s}}_{\mathrm{f}}^{\pm }=\pm p_{\mathrm{f}}{\mathbf{F}}^{\star }{\mathbf{e}}_1+6\eta \frac{q}{{\delta }^2}[({\mathbf{F}}^{\star }{\mathbf{e}}_1\cdot {\mathbf{e}}_1){\mathbf{e}}_2+({\mathbf{F}}^{\star }{\mathbf{e}}_1\cdot {\mathbf{e}}_2){\mathbf{e}}_1],$$

$$\text{mbox}\mathrm{div}\mathbf{S}=0,$$

$$\psi (\mathbf{F})=\frac{G_{\mathrm{e}}}{2}(\mathbf{F}\cdot \mathbf{F}-2\log J-3)+\frac{\Lambda }{2}(\log J{)}^2,$$

$$\mathbf{S}=\frac{\partial \psi }{\partial \mathbf{F}}=G_{\text{e}}\left(\mathbf{F}-\frac{1}{J}{\mathbf{F}}^{\star }\right)+\Lambda \frac{\log J}{J}{\mathbf{F}}^{\star }.$$

$$J=1+\Omega (c-c_{\text{o}}),$$

$$\dot{c}=-\text{mbox}\mathrm{div}\mathbf{h},$$

$$\psi (\mathbf{F},c)={\psi }_{\mathrm{e}}(\mathbf{F})+{\psi }_{\mathrm{m}}(c),$$

$${\psi }_{\mathrm{e}}(\mathbf{F})=\frac{G_{\mathrm{d}}}{2J_{\text{o}}}(J_{\text{o}}^{1/3}\mathbf{F}\cdot \mathbf{F}-3)$$

Full text

Keywords: Hydraulic fracture, Toughening, Multiple cracking, Brittle layer, Hydrogel, Cohesive zone.

1 Introduction

Nature has adopted diverse, remarkable strategies to enhance the flaw-tolerance of biological tissues, such as bone and tooth [7]. As a result, they can sustain relatively high levels of strain while maintaining their integrity. Epithelium is another example of such tissues. Despite the intrinsic brittle behavior of the cell monolayer [9], this can display a high fracture toughness, as it has recently been shown experimentally [3]. In a previous work [11] we have demonstrated that this behavior is determined by the hydraulic coupling between the epithelial layer and the extracellular matrix, which can be regarded as a poroelastic, hydrogel-like material. Specifically, because fracture requires flow into the crack to fill its volume, the kinetics of solvent migration within the hydrogel controls the velocity of crack propagation, so that decreased permeability promotes multiple-cracking at cell-cell junctions. Then, this distributed cracking mechanism maximizes the external work performed on the system before failure (which is a measure of toughness), since cell-cell separation requires a significant amount of work and is accompanied by dissipation due to solvent flow.

The just described behavior sharply contrasts with that of brittle materials, which are highly flaw-sensitive and typically fail catastrophically, by localized fracture [1]. For instance, the equilibrium of two edge cracks with the same length in an elastic brittle layer under tensile load is unstable, so that any perturbation will cause only one of the cracks to propagate dynamically [12]. In general, since toughening of soft materials typically relies on energy dissipation [14], studying the energetics of crack propagation [13] is crucial to understand the interplay of different physics involved in fracture phenomena and for the implementation of toughening strategies in synthetic materials [8]. Therefore, we here reconsider the problem studied numerically in [11], where a brittle layer containing two pre-cracks and bonded to a hydrogel substrate is subject to a remote strain. In particular, we analyze the relative influence of the different dissipative mechanisms on the macroscopic toughness of the system. Moreover, we focus on the effect of the strain rate on the transition from distributed to localized cracking.

The paper is organized as follows. We first establish the model accounting for the poroelasticity of the hydrogel substrate, the elasticity of the brittle layer and its hydraulic fracture. Then, the weak formulation of the model is presented, which is suitable for its implementation into a finite element code. Numerical results are shown and discussed concerning the dissipations for both the regimes of distributed and localized fracture.

2 Model for the hydraulic fracture of a brittle layer bonded to a hydrogel substrate

In this section, we derive the governing equations for the nonlinear model of hydraulically driven crack propagation in a brittle layer bonded to a hydrogel substrate. Both the brittle solid and the hydrogel are modelled as layers of finite thickness and infinite length. To simplify notation, the derivation is presented for the case of a single crack. Then, in the following section, the numerical model is extended to account for two competing cracks.

2.1 Notation and kinematics

We introduce an orthonormal basis $\{\mathbf{e}_{i}\}$, $i=\{1,2,3\}$, for the three-dimensional Euclidean space $\mathcal{E}$. We denote by $\mathcal{B}^{\textrm{g}}_{t}\subset\mathcal{E}$ and $\mathcal{B}^{\textrm{e}}_{t}\subset\mathcal{E}$ the current configurations at time $t$ of the hydrogel and the brittle layer, respectively; the corresponding reference configurations are indicated by dropping the subscript $t$. Further, we introduce two Cartesian coordinate systems $\{X_{i}\}$ and $\{x_{i}\}$, for the reference and the current configuration of the system, respectively. A material point $X$ of the reference solid domain $\mathcal{B}=\mathcal{B}^{\textrm{g}}\cup\mathcal{B}^{\textrm{e}}$ with coordinates $X_{i}$ is mapped into a spatial point $x$ with coordinates $x_{i}=f_{i}(X_{j},t)$ in the current configuration $\mathcal{B}_{t}$ (Fig. 1), where $f_{i}$ are the coordinate representations of the motion of the solid domain. We will use the symbol $\mathbf{F}$ for the deformation gradient, and write $J=\det\mathbf{F}$ for its determinant and $\mathbf{F}^{\star}=J\mathbf{F}^{-T}$ for its cofactor. We assume plane strain conditions, such that all the quantities do not depend on $X_{3}$ and thus $x_{\alpha}=f_{\alpha}(X_{1},X_{2},t)$, $\alpha=\{1,2\}$, whereas $x_{3}=f_{3}(X_{3},t)=X_{3}$. We take all of the three-dimensional domains to be of unit-depth, so that, in the ensuing derivation, volume integrals and area integrals over the cross-section corresponding to the cut plane $X_{3}=x_{3}=0$ coincide.

The brittle layer contains an initially closed pre-crack with length $a_{\rm{o}}$, whose faces belong to the referential segments $\mathcal{J}^{\pm}$, being $\mathcal{J}_{t}^{-}$ and $\mathcal{J}_{t}^{+}$ their current counterparts. For the reference configuration $\mathcal{L}$ of the longitudinal crack axis, we choose the segment aligned with the $X_{2}$-axis: $\mathcal{L}=\{X\in\mathcal{E}\,|\,X_{1}=X_{3}=0,X_{2}\in[0,h]=\mathcal{H}\}$, where $h$ is the thickness of the brittle layer. We assume that the axis $\mathcal{L}_{t}$ remains straight and vertical upon deformation, and that the crack undergoes a plane motion symmetric with respect to such axis. Thus, the axis of the crack stretches by the amount $\lambda(X_{2},t)=\partial f^{\rm c}_{2}/\partial X_{2}(X_{2},t)$ to follow the vertical motion $f_{2}^{\rm c}(X_{2},t)=\lim_{X_{1}\rightarrow 0^{\pm}}f_{2}(X_{1},X_{2},t)$, with $X_{2}\in\mathcal{H}$.

We define the material description of the crack opening $\delta(X_{2},t)$ as

[TABLE]

where the symbol $\llbracket\cdot\rrbracket$ denotes the jump operator, whereas $f_{1}^{\rm c}(X_{2},t)=\lim_{X_{1}\rightarrow 0^{+}}f_{1}(X_{1},X_{2},t)$ and the second equality follows from symmetry. By composing $\delta$ with the inverse of the deformation mapping, we obtain the crack opening as a function of the spatial location along $\mathcal{L}_{t}$: $\delta_{\rm s}(x_{2},t)=\delta(X_{2},t)\circ(f_{2}^{\rm c})^{-1}(x_{2},t)$. The velocity of the crack faces may be readily computed as

[TABLE]

where $w(x_{2},t)=\dot{f}^{\rm c}_{2}(X_{2},t)\circ(f_{2}^{\rm c})^{-1}(x_{2},t)$ is its vertical component. Here, a superposed dot denotes the material time derivative (at $X_{2}$ fixed), so that

[TABLE]

We call $\mathcal{B}^{\rm{f}}_{t}\subset\mathcal{E}$ the fluid domain, that is, the spatial volume with longitudinal axis $\mathcal{L}_{t}$ enclosed by the crack faces and filled with solvent. Finally, we denote by $\mathbf{v}(x,t),\ x\in\mathcal{B}_{t}^{\rm{f}}$, the spatial velocity field of the solvent within the crack.

2.2 Solvent flow within the crack

In this section, we establish a reduced one-dimensional model for the solvent flow within the crack. Specifically, we derive a set of equations, defined over the axis $\mathcal{L}$, by localizing the integral balance laws for the fluid expressed in terms of resultant quantities over horizontal cross-sections. A related model was introduced in [5] without accounting for stretching along the axis. We consider the steady flow of an incompressible, viscous fluid under the assumption that the motion is quasi-1D, i.e. the velocity field $\mathbf{v}$ is independent of $x_{3}$ and is dominated by the longitudinal component $v_{2}=\mathbf{v}\cdot\mathbf{e}_{2}$.

We define a control volume $\mathcal{V}_{t}\subset\mathcal{B}^{\textrm{f}}_{t}$, enclosed by the cross-sections that are located at two fixed positions $\bar{x}_{2}$ and $\bar{x}^{\prime}_{2}$ along $\mathcal{L}_{t}$, and the crack faces $\mathcal{J}_{t}^{\pm}$.

2.2.1 Conservation of solvent mass

The equation of conservation of mass for the incompressible fluid in $\mathcal{V}_{t}$ reads

[TABLE]

where $\mathbf{n}$ is the outward unit normal to $\partial\mathcal{V}_{t}$. On the crack faces this is given by the following relation

[TABLE]

with $\|\mathbf{n}_{\rm c}\|=\sqrt{1+1/4(\partial\delta_{\rm s}/\partial x_{2})^{2}}$. Notice that the area element along the crack faces in Eq. (4) transforms according to $\mbox{d}a=\|\mathbf{n}_{\rm c}\|\mbox{d}x_{2}\mbox{d}x_{3}$. Using eqs. (2) and Eq. (5) the volume flux contributions over parts of $\mathcal{J}_{t}^{\pm}$ in Eq. (4) may be evaluated as

[TABLE]

while the net flux through the cross-sections delimiting $\mathcal{V}_{t}$ is

[TABLE]

with

[TABLE]

the solvent volume flux per unit depth. By summing these two contributions, Eq. (4) may be recast as

[TABLE]

To express this equation in the reference configuration $\mathcal{L}$, we first observe that

[TABLE]

Then, upon exploiting Eq. (3) and localizing, we obtain

[TABLE]

where $q(X_{2},t)=q_{\rm s}(x_{2},t)\circ f_{2}^{\rm c}(X_{2},t)$ is the material description of the solvent flux

[TABLE]

relative to the material particles on the crack faces.

2.2.2 Balance of forces

Upon neglecting inertia, the balance of forces for the solvent within the crack along the longitudinal axis of $\mathcal{V}_{t}$ reads

[TABLE]

As a constitutive law for the stress tensor we take the representation for an incompressible, Newtonian viscous fluid, such that

[TABLE]

with $p_{\rm s}$ and $\eta$ the pressure and the viscosity of the solvent, respectively. Following lubrication theory [2], we assume that: i) the velocity profile is parabolic along the cross-sections, i.e.

[TABLE]

where $v_{2}^{\rm max}=(3Q/\delta_{\rm s}-w)/2$; ii) the pressure $p_{\rm s}(x_{2},t)$ is uniform over each cross-section; iii) the components of the velocity gradient along $x_{2}$ are negligible. Under these hypotheses, the traction on the crack faces is

[TABLE]

where in the last equality we have used the representation (16) evaluated at $x_{2}=\pm\delta_{\rm s}/2$. With this, the contribution to the integral in (14) extending over the crack faces may be computed as

[TABLE]

Moreover, the force resultant over the cross-sections at $\bar{x}_{2}$ and $\bar{x}^{\prime}_{2}$ is

[TABLE]

Finally, by expressing the sum of eqs. (18)-(19) in the reference crack axis $\mathcal{L}$ and by localizing we obtain

[TABLE]

where $p_{\rm f}(X_{2},t)=p_{\rm s}(x_{2},t)\circ f_{2}^{\rm c}(X_{2},t)$.

2.3 Cohesive zone model

To model crack propagation, we employ a cohesive zone approach [6]. We take the cohesive tractions $\mathbf{s}^{\pm}_{\rm c}$ per unit reference area acting on the crack faces $\mathcal{J}^{\pm}$ to be orthogonal to the reference longitudinal axis $\mathcal{L}$ of the crack, and such that their magnitude is a function of the crack opening:

[TABLE]

We choose a bilinear traction-separation law for $\sigma(\delta)$, such that

[TABLE]

where $\sigma_{\textrm{o}}$ is the cohesive strength (attained for $\delta=\delta_{\textrm{o}}$) and $d(X_{2},t)\in[0,1]$ is the damage. We prescribe the following evolution law for the damage:

[TABLE]

where $\delta^{\textrm{m}}(X_{2},t)$ is the maximum value attained by the opening during crack evolution and $\delta_{\textrm{c}}$ is the crack opening at failure ($d=1$). The fracture energy $\Gamma$ is related to the parameters of the cohesive zone model through the equation $\Gamma=\sigma_{\textrm{o}}\delta_{\textrm{c}}/2$.

In addition to the cohesive tractions, the crack faces $\mathcal{J}_{t}^{\pm}$ are subject to the pressure and to the shear stress exerted by the fluid, which contributes to the traction with the term $-\mathbf{t}_{\rm f}^{\pm}$, see Eq. (17). The corresponding reference traction is

[TABLE]

such that the total reference traction acting on the crack faces is given by $\mathbf{s}^{\pm}=\mathbf{s}_{\rm f}^{\pm}+\mathbf{s}^{\pm}_{\rm c}$.

2.4 Elasticity of the brittle layer

By assuming inertia negligible, the balances of forces and moments in $\mathcal{B}^{e}$ read

[TABLE]

where $\mathbf{S}$ denotes the first Piola-Kirchhoff stress tensor.

We model the brittle layer as an impermeable, elastic solid characterized by the compressible neo-Hookean free energy density

[TABLE]

where $G_{\textrm{e}}$ and $\Lambda\gg G_{\textrm{e}}$ are the Lamé moduli. The corresponding first Piola-Kirchhoff stress is given by

[TABLE]

Notice that, because of Eq. (27) and the plane strain hypothesis, $\mathbf{S}\mathbf{e}_{1}\cdot\mathbf{e}_{3}=\mathbf{S}\mathbf{e}_{2}\cdot\mathbf{e}_{3}=\mathbf{S}\mathbf{e}_{3}\cdot\mathbf{e}_{1}=\mathbf{S}\mathbf{e}_{3}\cdot\mathbf{e}_{2}=0$ and the out-of-plane balance equation in (25)1 is trivially satisfied. Hence, the balance of forces and moments may be formulated in terms of the plane components of $\mathbf{S}$ only.

2.5 Poroelasticity of the hydrogel layer

In the following, we briefly summarize the non-linear swelling theory for hydrogels introduced in [10]. The state of the hydrogel is described by the motion $f$ of the polymer network and the solvent concentration $c$ per unit reference volume. The chemical potential $\mu$ of the solvent within the hydrogel quantifies the energy carried by the solvent and represents the driving force of solvent migration. The corresponding solvent molar flux $\mathbf{h}$ characterizes the relative motion of the solvent with respect to the polymer matrix. Consistently with the plane strain hypothesis, we assume that solvent migration takes place in the plane $\mathbf{e}_{1}$-$\mathbf{e}_{2}$, so that $c(X_{1},X_{2},t)$, $\mu(X_{1},X_{2},t)$, and $\mathbf{h}\cdot\mathbf{e}_{3}=0$.

The polymer matrix and the solvent are considered to be separately incompressible; hence, the change in volume of the hydrogel is related to the change in solvent concentration:

[TABLE]

where $\Omega$ is the solvent molar volume and $c_{\textrm{o}}=(J_{\textrm{o}}-1)/(\Omega J_{\textrm{o}})$ is the solvent concentration per unit reference volume associated to the initial free swelling from the dry configuration to $\mathcal{B}^{\textrm{g}}$. Here, $J_{\textrm{o}}$ is the volume ratio between the reference and the dry configuration. The constraint of Eq. (28) is enforced through the Lagrange multiplier $p$.

Swelling processes are governed by the equations of balance of forces and moments, which are the same as those for the brittle layer, see Eqs. (25), and by the balance of solvent mass in $\mathcal{B}^{\textrm{g}}$

[TABLE]

subject to the initial condition $c(X_{1},X_{2},0)=c_{\rm o}$.

As concerns the constitutive equations, we prescribe the following Flory-Rehner representation for the free energy density of the hydrogel [4, 10]

[TABLE]

where

[TABLE]

and

[TABLE]

are the neo-Hookean elastic energy and the Flory-Huggins free energy of mixing, respectively. Here, $G_{\textrm{d}}$ is the shear modulus of the dry polymer, $\mathcal{R}$ is the universal gas constant, $T$ is the absolute temperature of the environment, and $\chi$ is the polymer-solvent mixing parameter. For the consistency with thermodynamical principles, the corresponding constitutive equations are given by:

[TABLE]

where $D$ is the diffusivity of the solvent within the hydrogel. In the reference configuration $\mathcal{B}^{\textrm{g}}$, the hydrogel is in equilibrium with the external solvent whose chemical potential is $\mu_{\textrm{o}}$. Specifically, we consider an external solvent in equilibrium with its vapor, such that $\mu_{\textrm{o}}=0\ \mbox{J/mol}$. Thus, the reference configuration is identified by the conditions $\mathbf{F}=\mathbf{I}$, $\mu=\mu_{\textrm{o}}$ and $\mathbf{S}=\mathbf{0}$, which together determine the swelling ratio $J_{\textrm{o}}$ as the solution of following equation

[TABLE]

2.6 Boundary and interface conditions

In the numerical implementation, we refer to a computational domain which is limited along the $\mathbf{e}_{1}$-axis by vertical boundaries of normal $\mathbf{m}=\pm\mathbf{e}_{1}$. In order for the computational model to appropriately approximate the theoretical setting, we choose its length to be sufficiently larger than the thickness of the composite.

Starting from the equilibrium state of free swelling with $\mu=\mu_{\rm o}$, the system is subject to a homogeneous, horizontal strain of constant rate $\dot{\varepsilon}$. Correspondingly, the vertical component $\mathbf{u}\cdot\mathbf{e}_{1}$ of the displacement field $\mathbf{u}(X_{1},X_{2},t)$ is prescribed along the vertical boundaries $\partial_{\mathbf{u}}\mathcal{B}$ of the composite, while the vertical component $\mathbf{u}\cdot\mathbf{e}_{2}$ of the displacement field is constrained at the bottom boundary of the hydrogel only. In order to ease crack opening, continuity of the displacements is relaxed at the interface between the brittle layer and the hydrogel along small segments near the cracks. The remaining parts of the boundary $\partial\mathcal{B}$ are traction-free.

As regards the boundary conditions for solvent migration, we assume zero solvent flux on $\partial\mathcal{B}^{\textrm{g}}$. For short times, no significant exchange of solvent with the exterior can take place. At the crack inlet, we impose the continuity of solvent pressure: $\mu=\Omega p_{\rm f}$, where $\mu/\Omega$ is the pressure within the hydrogel. At the crack tip, we impose zero solvent flux: $q=0$.

2.7 Power balance

As mentioned in the introduction, we are interested in analyzing the breakdown of dissipative contributions upon the toughness enhancement of the system. To this aim, we write the total power balance for the brittle layer–hydrogel composite as

[TABLE]

where

[TABLE]

are the total free energy of the system, the mechanical power expended by the applied tractions, the power expended by the cohesive tractions in fracture processes, the power dissipated by transport of the fluid within the cracks, and the power dissipated by solvent transport within the hydrogel, respectively. Notice that there is no contribution associated to solvent transport across the boundary of the system because we have assumed the hydrogel to be impermeable for short times.

2.8 Weak form of the governing equations

In order to solve the problem set in the previous paragraphs by means of the finite element method, we recast the governing equations (12), (25)1, (28), and (29) in weak form. The continuity of the solvent pressure at the crack inlet (subscript ‘in’) is enforced through the Lagrange multiplier $g$ corresponding to the (concentrated) solvent mass flux per unit length. Then, the weak formulation of the problem reads: find $p_{\rm f}$, $\mathbf{u}$, $c$, $p$ and $g$ such that the following equations

[TABLE]

hold for arbitrary test fields (indicated with a superposed tilde) compatible with the Dirichlet conditions. The weak form of the governing equations is complemented by Eq. (20) and by the constitutive relations (33)-(35). This, and the corresponding boundary conditions, have been implemented into the finite element software COMSOL Multiphysics v5.2. Specifically, quadratic shape functions were used for all the unknown fields, except for the pressure field $p$, which was discretized using linear shape functions to get a reliable approximation of the volume constraint. The implicit, variable-order (from 1 to 5), adaptive step-size BDF solver was used for time-stepping. A quasi-Newton algorithm was employed to solve iteratively the non-linear algebraic system resulting from the finite element discretization at each time step. The direct solver MUMPS was chosen for the solution of the linearized system at each iteration. The mesh consisted of about $2\times 10^{4}$ triangular elements corresponding to about $10^{5}$ degrees of freedom, and was made symmetric with respect to the vertical symmetry axis of the computational domain to avoid introducing any numerical bias in the distribution of solvent flux between the cracks. Local mesh refinement along the crack paths was performed to ensure that the cohesive zones were discretized with at least 20 elements. Notice that the faces $\mathcal{J}^{+}$ and $\mathcal{J}^{-}$ are modeled as distinct (but overlapping) segments in the numerical model where the tractions $\mathbf{s}^{\pm}$ are prescribed.

3 Numerical results

We apply the just described model to the analysis of the ideal case of a brittle layer containing two edge pre-cracks of initial length $a_{\rm o}$ and separation $s$. To break symmetry, we introduce a 10% difference between the toughnesses $\Gamma_{\rm R}$ and $\Gamma_{\rm L}$ of the right and left crack, respectively.

Dimensional analysis dictates that any relevant quantity, such as the total dissipated energy, has to depend on the following set of dimensionless groups: $h/a_{\rm o}$, $h/s$, $h/H$, $G_{\rm e}/\Lambda$, $G_{\rm d}/J^{1/3}_{\rm o}G_{\rm e}$, $G_{\rm d}\Omega/\mathcal{R}T$, $\Gamma_{\rm R}/\Gamma_{\rm L}$, $\chi$, $\mu_{\rm o}/\mathcal{R}T$, $|\dot{\varepsilon}|s^{2}/D$, $\tau$, $G_{\rm e}a_{\rm o}/\Gamma_{\rm L}$, $\Gamma_{\rm L}a_{\rm o}/\eta D$. Here, $H$ is the thickness of the hydrogel substrate, whereas $\tau=|\dot{\varepsilon}|t$ is a dimensionless measure of time.

For the computational study, we set the values of part of these dimensionless groups as reported in Table 1. The first three are geometrical ratios, while the following three are related to the elastic moduli of the brittle layer and the hydrogel. We recall from [11] that crack propagation under compression (tension) requires the ratio $3G_{\rm d}/4J^{1/3}_{\rm o}G_{\rm e}$ between the moduli of the hydrogel and the brittle layer to be higher (lower) than $1/2$. In compression, for instance, this condition allows the solvent pressure within the cracks to overcome the compressive stresses due to the remotely applied strain. Thus, we set $G_{\rm d}/J^{1/3}_{\rm o}G_{\rm e}=0.38$, $G_{\rm d}\Omega/\mathcal{R}T=6\times 10^{-5}$ to study cracking in tension, and $G_{\rm d}/J^{1/3}_{\rm o}G_{\rm e}=3.82$, $G_{\rm d}\Omega/\mathcal{R}T=4\times 10^{-4}$ for the simulations of cracking in compression. Given $\chi$, $G_{\rm d}\Omega/\mathcal{R}T$ and $\mu_{\rm o}/\mathcal{R}T$, the initial swelling ratio $J_{\textrm{o}}$ is computed from Eq. (36). Eventually, by allowing the remaining parameters to vary, we regard a dimensionless, physical quantity as a function of $\Pi=|\dot{\varepsilon}|s^{2}/D$ and $\tau$.

We now focus on the significance of the dimensionless group of $\Pi$ upon the fracture behavior of the system. In a previous study [11], we have demonstrated that decreased diffusivity conveys toughening by promoting multiple-cracking. Here, in light of dimensional analysis, we deduce that the transition from localized to distributed fracture is actually controlled by the dimensionless group $\Pi$. The results in Fig. 2, where we report snapshots of the system at failure, confirm that such a transition is associated to an increase in the value of $\Pi$, both in tension and compression. Further evidence is provided in Figs. 3(a),(d), where the relative length difference $\Delta\tilde{a}=(a_{\rm R}-a_{\rm L})/a_{\rm o}$ between the cracks sharply grows for $\Pi=2\times 10^{-7}$, as rapid, single crack propagation occurs. Therefore, multiple-cracking may be obtained by either an increase in strain rate or separation between the cracks, or by a decrease in diffusivity of the solvent within the hydrogel.

The just described behavior may be rationalized as follows. In a brittle solid, the mechanism of crack tip shielding is responsible for localized, dynamic fracture. Any offset between the crack tips implies a progressive decrease (increase) in the energy release rate of the lagging (leading) crack [12]. In the presence of the hydrogel, crack advance always requires solvent transport and is thereby controlled by its diffusivity. Hence, continued loading that is rapid with respect to the kinetics of solvent transport can favour distributed cracking by sustaining the driving force of both fractures.

Multiple-cracking results in an enhancement of the macroscopic toughness of the system. To quantify such an enhancement, we report in Figs. 3(b),(e) the dimensionless mechanical work $\tilde{W}_{\rm t}=W_{\rm t}/(G_{\rm e}h+G_{\rm d}H)$ performed on the system as computed by time-integration of Eq. (39). Notice that, an increase of two orders of magnitude in the strain rate produces more than a three-fold increase in the work at failure, both in tension and in compression. This trend is motivated by the increase of the dissipation $\tilde{E}_{\rm d}=E_{\rm d}/(\Gamma_{\rm L}+\Gamma_{\rm R})a_{\rm o}$ that accompanies distributed cracking, see Figs. 3(c),(f). Here, $E_{\rm d}$ collects the time-integrals of the powers dissipated by fracture, Eq. (40), and solvent transport, Eq. (42), as the energy dissipated by viscous flow within the cracks is comparably negligible. Interestingly, we notice that dissipation due to solvent transport mainly contributes to such an increase in $\tilde{E}_{\rm d}$. As a consequence, the relative contribution of the work of fracture to $\tilde{E}_{\rm d}$ decreases in the transition from localized to distributed cracking.

4 Conclusions

Motivated by recent experimental results on the fracture of epithelial cell monolayers adhered to an hydrogel substrate, a model has been developed that allows for the analysis of hydraulically driven cracking of this system. We have shown that such a composite can develop a toughening mechanism that relies on the multiple-cracking of the brittle phase. In particular, we have demonstrated that the transition from localized to distributed cracking is determined by either a decrease in solvent diffusivity or an increase in loading rate. A detailed energy analysis has quantified the dissipative processes associated with multiple-cracking. Future work will investigate upon the possibility of implementing this concept in the engineering of flaw-tolerant, biomimetic materials.

Acknoledgements

The authors acknowledge the support of the European Research Council (AdG-340685 MicroMotility) and of National Group of Mathematical Physics (GNFM-INdAM) through the initiative “Progetto Giovani”. The authors also thank prof. Robert M. McMeeking for useful discussions on the subject.