Coupled flow and deformation fields due to a line load on a poroelastic half space: effect of surface stress and surface bending

1. Introduction

Gels are ubiquitous and have diverse applications in clinical devices [1,2], tissue engineering [3–5] and soft robotics [6–8]. They are essentially a polymer network swollen in solvent (e.g. water). The three-dimensional network is formed by cross-linked polymer chains, and the small solvent molecules can transport through the network. For example, when a solvent-saturated gel is suddenly brought into contact with a rigid indenter, the solvent cannot immediately flow out of the network, so the gel initially behaves as an incompressible material. This gives rise to a pressure gradient in the solvent driving the solvent flow until the pressure goes to zero everywhere, and all the stresses are transferred to the polymer network.

Classical theory on solvent diffusion in elastic solids dates to Terzaghi [9] who studied one-dimensional consolidation of clay soils. The general theory of poroelasticity was developed by Biot [10–13]. Solutions of the linearized theory can be found in the works of McNamee & Gibson [14,15], Verruijt [16] and Rice & Cleary [17]. Recent interests in hydrogels have motivated many researchers to study poroelasticity in elastic gels using modern continuum mechanics methods; see for example, Hong et al. [18], Doi [19], Duda et al. [20], Chester & Anand [21,22], Wang & Hong [23]. In particular, Bouklas & Huang [24] provide the relationship between the linear theory of poroelasticity by Biot [10] and the nonlinear continuum theory for hydrogels proposed by Hong et al. [18]. However, none of these methods considers the effect of surface stress on solvent flow.

Surface stress effects are typically felt over a characteristic length scale, σ/E, where E is Young's modulus of the bulk and σ is the magnitude of the surface stress. For soft solids, such as gels with elastic modulus in the kPa range and surface stress in the mN m⁻¹ range (e.g. E ≈ 35.6 kPa, σ ≈ 70 mN m⁻¹ for a gelatin-based organogel sample [25], and E ≈ 5.6 kPa, σ ≈ 20 mN m⁻¹ for a silicone gel at zero surface strain [26]), the corresponding characteristic length is of the order of micrometres. Thus, surface stress could flatten sharp features by smoothing corners and undulations [25,27–29]; drive instabilities [30]; stiffen fluid–solid composites [31]; and significantly affect the opening of a crack in a soft material [32–34]. The resistance to indenting a surface by a rigid particle can be governed more by surface stress than by the elasticity of the material. Thus, the contact and adhesive theory by Hertz and Johnson–Kendall–Roberts must be modified [26,35–39]. Surface stress can also play an important role in cell mobility. Recent experiments strongly suggest that biological cells respond to differences in surface stresses more than substrate elasticity [40]. The focus of the research community is how surface stress affects mechanical response of soft elastic solids. However, many gels are poroelastic, and to the best of our knowledge, there has been no theoretical analysis on how surface stress affects solvent flow. This work is a small step in this direction.

In this work we consider a canonical problem in mechanics—the transient response of a poroelastic half space subjected to a line load. The surface of the half space can resist deformation since it has a constant surface stress and bending stiffness. To motivate this problem, we name two examples. In wetting, the vertical component of a droplet's surface stress (locally a line load) pulls the surface of the soft substrate upwards, leading to the formation of a ridge at and near the applied line load [41,42]. The dynamics of contact line motion has been examined by modelling the substrate as viscoelastic [43–45]. The axisymmetric version of this problem corresponds to loading by a circular line load. A second example is the experiment of Berman et al. [46]. In their experiment, a hard sphere (silica bead) is brought into adhesive contact with a silicone gel substrate. The sphere is then retracted quasi-statically until the contact line becomes unstable and slides toward a ‘point contact’ on the bottom of the sphere (figure 1a), from which the gel detaches. They studied the dynamics of the gel's recovery to the flat geometry (figure 1b). They suggested that the recovery of the peak is consistent with a relaxation process driven by surface stresses and slowed by solvent flow through the gel network. In the discussion we use our line load solution to study a simplified two-dimensional (2D) version of this problem.

• Download figure
• Open in new tab
• Download PowerPoint

The key assumptions of this work are: (i) the substrate is linear poroelastic, (ii) surface stress is constant and isotropic, and (iii) surface bending stiffness is constant. Despite the linearity assumption, poroelasticity problems are very difficult to solve and exact, closed-form solutions are rare. An advantage of including surface bending in our surface constitutive model is that it regularizes the singular strain and stress fields due to the line load [47,48]. In addition, there are many situations where surface bending can be important [25,49–51], e.g. the surface of cells consists of lipid bilayers which are known to resist deformation by bending [52,53]. However, in this work we focus on the limit where surface bending stiffness is small.

The plan of the paper is as follows. Section 2 states and formulates the problem. The solution is presented in §3, and numerical results are discussed in §4. In §5 we revisit the experiments and scaling arguments by Berman et al. [46], and in §6 we conclude with a short summary.

2. Problem statement and formulation

The geometry is shown in figure 2 where an isotropic, homogeneous, fully saturated, elastic, porous medium occupies the half space |x| < ∞, 0 ≤ z < ∞. A vertical compressive line force with magnitude N > 0 is suddenly applied at the origin at time t = 0⁺ (figure 2). We assume plane strain deformation where the out-of-plane displacement vanishes, and the in-plane horizontal and vertical displacements u and w depend only on the in-plane coordinates x and z. The undrained gel is assumed to be incompressible and the drained gel has Poisson's ratio ν. The shear modulus of the gel is denoted by μ. The in-plane stresses are denoted by σ_xx, σ_zz and σ_xz, and the pore pressure by p. The in-plane strains are denoted by ε_xx, ε_zz and ε_xz. The stresses acting on the network, often referred to as effective stresses, are denoted by ${\sigma }_{xx}^e,$${\sigma }_{zz}^e$ and ${\sigma }_{xz}^e.$ The relation between stresses, effective stresses and pore pressure p is

• Download figure
• Open in new tab
• Download PowerPoint

We adopt the sign convention of McNamee & Gibson [14] in which compressive stresses and strains are positive. For the fluid phase, the decrement of fluid content, defined as the change of fluid volume per unit reference volume, is denoted by ξ where ξ = ε_xx + ε_zz due to the plane strain assumption. The coefficient of permeability in Darcy's law is denoted by k, and the fluid viscosity is η. Prior to loading the gel is assumed to be in equilibrium with the environment. We assume the gel does not dry or further swell if left at such an environment. The pore pressure p can be taken as zero everywhere in this initial state.

In contrast to standard poroelasticity, the surface of the half space can resist deformation by stretching and bending. Stretching of the surface is resisted by a constant isotropic surface stress σ which is independent of the surface stretch. The bending stiffness of the surface is denoted by D. The traction discontinuity across the interface is ${\sigma }_{zz}(x,z=0,t)-N{\delta }^D(x),$ where δ^D(x) is the Dirac delta function. Static equilibrium requires this traction jump to be balanced by surface bending moment and surface stress [47,54], i.e.

The first term on the r.h.s. of (2.2a) represents the pressure supported by the bending moment, whereas the second term is the curvature induced Laplace pressure due to surface stress. The surface z = 0 is assumed to be permeable and shear stress free, that is, for all t > 0,

and The initial condition is

The surface properties, surface stress σ and surface bending stiffness D, introduce two length scales over which effects of surface stress or bending can be felt. The characteristic length induced by surface stress is the elasto-capillary length and is l_c = σ/2μ. The characteristic length induced by surface bending is l_b = (D/2μ)^1/3; we call l_b the capillary-bending length. The ratio of these two lengths κ ≡ l_c/l_b is a dimensionless parameter that controls the relative importance of surface stress to surface bending. When κ ≡ l_c/l_b ≫1, surface bending effect can only be felt over a very small region underneath the line load and outside this region surface stress dominates. We will focus on this case in this work.

Before diving into our analysis, let us consider two limiting situations: the response of the poroelastic substrate at t = 0⁺ (incompressible limit) and at times that are very long compared with the characteristic diffusion time (drained time limit). The two cases are also known as the undrained and drained states of the poroelastic substrate, respectively. In the first case, the substrate behaves as an incompressible linear elastic solid since solvent flow is a time-limited process. In this limit, the solution is known [55]. If the surface is governed by a constant surface stress with no bending stiffness (κ → ∞), then the surface kinks with an angle determined by ‘Neumann's triangle’ [56]. The strains in the substrate are bounded but some of the stress components are unbounded with a logarithmic singularity. If bending stiffness is non-zero, then the stresses and strains are bounded and continuous. On the other hand, in the long-time limit where the fluid is drained and the pore pressure vanishes, the substrate behaves as a compressible elastic solid. Again, the solution for this case is known: for the case of zero bending stiffness (κ → ∞), the stresses and strains have a weak logarithmic singularity [57]. However, if bending stiffness is non-zero, then the stresses and strains are continuous and bounded. Recall that if surface stress and bending are both absent, the stresses and strains in the substrate will exhibit a non-integrable 1/r singularity in both limits where r is the distance from the origin.

3. Methods

We use the formulation of McNamee & Gibson [14] where the displacements u and w, decrement of fluid content ξ, pore pressure p and total stresses σ_xx, σ_zz and σ_xz are determined by two displacement potentials e and s:

where and ${\mathrm{\nabla }}^2()$ stands for the Laplacian operation ${\mathrm{\nabla }}^2()={\mathrm{\partial }}^2()/\mathrm{\partial }{x}^2+{\mathrm{\partial }}^2()/\mathrm{\partial }{z}^2$. The displacement potentials, e and s, satisfy and where

To expedite the analysis, we use the following normalization: the normalized coordinates are X = x/l_b, Z = z/l_b; the normalized displacements are U = 2μu/N, W = 2μw/N; the normalized stresses are P = pl_b/N, Σ_XX = σ_xxl_b/N, Σ_ZZ = σ_zzl_b/N, Σ_XZ = σ_xzl_b/N; and the normalized decrement of fluid content is $\mathrm{\Xi }=2\mu l_{\text{b}}\xi /N$. Also, we introduce a time scale $l_{\text{b}}^2/c$, and rescale time using a normalized time $T=ct/l_{\text{b}}^2.$ The displacement potentials are normalized as $E=2\mu e/l_{\text{b}}N,$S = 2μs/N. Equations (3.1a–g) in these normalized variables are given in the electronic supplementary material.

We define the Fourier transformation (FT) of a function f(X, Z, T) by

where Λ is the transform variable and the tilde denotes FT. We also define the Laplace transformation (LT) of the function $\tilde{f}(\Lambda ,Z,T)$ by where Θ is the transform variable for the LT and the bar indicates LT.

The bounded solution of the normalized, transformed displacement potentials is found to be (derivation is provided in the electronic supplementary material)

and where the C_i's (i = 1, 2, 3) are constants to be determined. The coefficients C_i and the FT–LT of displacements and stresses are obtained by using the transformed boundary conditions. Details are given in the electronic supplementary material. We obtain the displacements or stresses in the physical domain by inverting the LT and FT, i.e. where Br stands for the Bromwich–Wagner contour. Due to the complexity of the transformed displacements and stresses (see electronic supplementary material (A5a–g)), we have not been able to determine these inverses in closed form. Instead, we invert these transforms numerically using the improved Talbot's method by Abate & Whitt [58] for the inverse LT and Gaussian quadrature for the inverse FT.

We are particularly interested in the vertical surface displacement. The vertical displacement W is given by (A1b) in the electronic supplementary material and on the surface the first two terms vanish; this results in

As is well known, the displacement field of an infinitely extended elastic body due to a line load is not well defined, since at infinity it diverges logarithmically. To bypass this limitation, a reference point on the surface is used as a datum where the displacement there is set to zero. Of course, this has no effect on the stress and strain distribution.

4. Results

We first consider the limiting case where there is no surface effect. Surprisingly, we have not been able to find the solution of this special case in the literature. Surface stress effect can be eliminated by setting κ = 0. Surface bending effect can be made insignificant by setting l_b small and considering distances away from the line load much larger than l_b. In other words, if we observe at the scale of a new characteristic length l_n = 10l_b, then the surface bending effect is neglectable. For this case, we renormalize all positions, stresses and time by l_n, N/l_n and $l_{\text{n}}^2/c$, respectively. In this way, the solution obtained from inversion can be expressed using these new normalized quantities:

(a) Vertical surface displacement

The vertical surface displacement is plotted in figure 3a for different times T_n using a drained Poisson's ratio ν = 0.22 [59]. Unless otherwise stated, in this and the following calculations we use X_n = 20, Z_n = 0 as our reference point where the vertical displacement is set to zero. As expected, the surface displacement increases over time due to solvent flowing away from the origin. We also compare this result with the predictions of small strain linear elasticity (SSLE) theory for homogeneous materials. Recall that the instantaneous response of the porous medium behaves like an incompressible solid, while at sufficiently long times the porous medium behaves as a compressible solid with the drained Poisson's ratio ν = 0.22. The normalized surface displacement predicted by SSLE is given by [55]

• Download figure
• Open in new tab
• Download PowerPoint

The SSLE predictions for undrained and drained (ν = 0.22) states are plotted as dotted and dashed lines, respectively, in figure 3a. Clearly, our numerical results at short and long times agree very well with (4.2) except for a very tiny region close to the line load where surface bending dominates.

Next, we investigate the surface stress effect on the poroelastic half space. To ensure that the effect of surface stress is visible at the length scale of l_n and dominates over surface bending moment, we set κ = 10 so that l_c = 10l_b = l_n. Figure 3b shows the normalized vertical surface displacement. Comparing with the special case where there is no surface effect (figure 3a), the surface displacement is significantly reduced. As expected, the displacement forms a kink at the origin, consistent with previous studies [41,42]. Wu et al. [57] have provided an exact solution for this ‘kink’ for an elastic solid where the surface stress is constant. Using our notation, the normalized vertical surface displacement is

where Ci( ) is the Cosine integral and si( ) is the Sine integral minus π/2. We plot the instantaneous and equilibrium solutions as dotted and dashed lines, respectively, in figure 3b. Again, they match our numerical results perfectly in these two limits.

(b) Effect of surface stress on pore pressure and flow

Here we study how surface stress alters the pore pressure and solvent flow in gels. In figure 4 we plot contours of normalized pore pressure within a rectangular region $[-5,5]\times [0,5]$ underneath the line load. Solutions with surface stress (κ = 10) and without surface stress (κ = 0) are represented by red solid and blue dashed lines, respectively. At short times T_n = 10⁻⁸ and T_n = 10⁻² (figure 4a,b), the contours congregate near the origin, indicating that the pore pressure close to the line load changes rapidly at short times. At these very short times, the maximal pore pressure occurs roughly at X_n = 0, Z_n = 5.4 × 10⁻⁴. The peak normalized pore pressures P_n are 1.1 (κ = 10) and 3.8 (κ = 0) respectively. Thus, surface stress reduces the peak pore pressure. On the other hand, figure 4a,b shows that surface stress tends to spread the pore pressure distribution—there is less variation in pore pressure in the horizontal direction. As expected, the pore pressure is less affected by the surface stress at distances much larger than the elasto-capillary length. As time increases, the maximal pore pressure occurs at larger values of Z_n; for example, at T_n = 1 (figure 4c), the maximal pore pressures occur at Z_n = 2.2 (κ = 10) and 2.0 (κ = 0). The peak values are 0.09 (κ = 10) and 0.13 (κ = 0). For sufficiently large time T_n = 100 (figure 4d), the pore pressure is negligible everywhere.

• Download figure
• Open in new tab
• Download PowerPoint

The gradient of pore pressure can be readily computed and using Darcy's law we determine the flux of solvent flow j. We normalize the flux by $Nk/\eta l_{\text{n}}^2$, i.e. $\mathbf{\text{J}}=\mathbf{\text{j}}\eta l_{\text{n}}^2/Nk$. The normalized flux J is

where e₁ and e₂ are the unit basis vectors for Cartesian coordinates. We isolate a small box region $[-1,1]\times [0,1]$ and compute the normal component of the flux on the boundary of this box at different times. Specifically, at the top and bottom boundaries, the normal components are $\mathbf{\text{J}}\cdot {\mathbf{\text{e}}}_2\equiv J_2=-\mathrm{\partial }{P}_{\text{n}}/\mathrm{\partial }{Z}_{\text{n}}$ and $\mathbf{\text{J}}\cdot (-{\mathbf{\text{e}}}_2)\equiv -J_2=\mathrm{\partial }{P}_{\text{n}}/\mathrm{\partial }{Z}_{\text{n}}$, respectively; at the left boundary is $\mathbf{\text{J}}\cdot (-{\mathbf{\text{e}}}_1)\equiv -J_1=\mathrm{\partial }{P}_{\text{n}}/\mathrm{\partial }{X}_{\text{n}}$; and at the right boundary is $\mathbf{\text{J}}\cdot {\mathbf{\text{e}}}_1\equiv J_1=-\mathrm{\partial }{P}_{\text{n}}/\mathrm{\partial }{X}_{\text{n}}$. In figure 5, we plot ±J₁ or ±J₂ along the boundaries at different times. Due to symmetry, we plot the flux for (κ = 0) on the right-half of the box (blue vectors) and (κ = 10) on the left-half (red vectors). For κ = 0 and at very short times, T_n = 10⁻⁸, there is considerable flow out of the substrate surface. The magnitude is very large at the origin (of the order of 10 000 times $Nk/\eta l_{\text{n}}^2$) and decays rapidly along the surface. On the other boundaries, the normal components of flux are relatively small. The simple conclusion is that there is net solvent flow out of the substrate. It is evident that surface stress dramatically reduces the outward flow through the surface; for example, the magnitude of flux at the origin is reduced by nearly 75%. However, the decay of flux along the surface is slower. On the interior boundary, the outward flow is negligible. At T_n = 10⁻², the magnitude of flux reduces to the order of $Nk/\eta l_{\text{n}}^2$. On all the boundaries, liquid flows outward. Again, surface stress reduces flow. At the transition time T_n = 1, the flow decreases further, and a different feature emerges: flow at the bottom boundary is inward. This is because the position of the maximum pore pressure now lies below Z_n = 1, and this drives the flow upwards across the bottom boundary. However, there is still overall fluid lost from this isolated region to the outside, since more fluid flows through the top surface to the environment. For sufficiently long times, the flux moving outwards on the top and inward on the bottom are nearly identical and negligible.

• Download figure
• Open in new tab
• Download PowerPoint

(c) Time-dependent stresses

Finally, we study the distribution of normalized stress ${\mathrm{\Sigma }}_{ZZ}^{\text{n}}$ and ${\mathrm{\Sigma }}_{XX}^{\text{n}}$ along the surface. In the following κ = 10. Since the surface is permeable (the pore pressure is zero), the stresses equal the effective stresses. Figure 6a plots ${\mathrm{\Sigma }}_{ZZ}^{\text{n}}$ along the surface at four different times. As expected, the maximum of ${\mathrm{\Sigma }}_{ZZ}^{\text{n}}$ occurs right underneath the line load and decreases with increasing distances from the line load. Also, ${\mathrm{\Sigma }}_{ZZ}^{\text{n}}$decreases over time. In figure 6b we plot the distribution of ${\mathrm{\Sigma }}_{XX}^{\text{n}}$ along the surface. Similar to the distribution of ${\mathrm{\Sigma }}_{ZZ}^{\text{n}}$, the maximum of ${\mathrm{\Sigma }}_{XX}^{\text{n}}$ occurs at the origin and decreases with increasing distances. However, the magnitude of the maximum value increases with increasing time. Lastly, note that a linear–log scale is used in figure 6a,b. The region where bending dominates is |X_n| < 0.1. Between 0.1 and 1 surface stress dominates, and the stress is proportional to $\ln X_{\text{n}}$ as predicted by the linearized theory [57]. This is a much weaker singularity compared to the SSLE prediction where the stress scales with 1/X_n.

• Download figure
• Open in new tab
• Download PowerPoint

5. Discussion

Our solution can be used to lend insight in the experiment and analysis by Berman et al. [46]. To do this, we review their experiment and scaling analysis. As illustrated in figure 1b, the detachment distance δ(t) is measured from the final contact point to the surface peak. Their experimental observation indicates a scaling relation between the detachment distance and time, where $\delta \sim t^{1/4}$. This behaviour is explained by Berman et al. [46] using a scaling argument. If we use the same scaling argument assuming that their experiments were performed under plane strain conditions, for example, on a very long cylinder, then the scaling between the detachment distance and time is found to be $\delta \sim t^{1/3}$. The reasoning is as follows: the Laplace pressure p^L associated with the surface curvature at the peak (figure 1b) is estimated by

which drives flow into the gel. The flux j associated with this flow is assumed to be

As stated by Berman et al., they only observe the peak and surroundings losing volume, therefore fluid must be flowing out of a region with a characteristic size L (in [46], L is taken to be the elasto-capillary length). According to (5.1a,b), the magnitude of flow velocity $|\mathbf{\text{j}}|$ should scale as

Note that the argument presented so far does not depend on whether the experiment is done in three dimensions (sphere) or two dimensions (cylinder). The volume rate of fluid flowing out of this region, $\dot{V},$ scales with the flow velocity in (5.2) multiplying by the area. For the 2D cylinder problem, $\dot{V}$ is simply the rate of volume flow per unit length of the cylinder, so

Note that in two dimensions, L cancels out since the area scales with L. The volume of fluid lost per time per unit length of cylinder from the tip of the peak is

Combining equations (5.3) and (5.4), the scaling law is

To check the scaling in equation (5.5), we simulate the detachment process of the gel from the rigid silica bead using the superposition strategy illustrated in figure 7. Without loss in generality, we assume that detachment occurs at T_n = 0⁺. To simulate the experiment, we first applied a constant tension line force −N at the origin at some negative time, say at $-T^{\mathrm{\infty }}$. $T^{\mathrm{\infty }}$ is chosen to be sufficiently large so the peak is formed in a fully drained gel: this step approximates the stable attachment between the silicone gel and the silica bead. We then applied a compressive line force N at T_n = 0⁺ to remove the tension force: this step corresponds to the detachment which occurs at T_n = 0⁺. We superimpose these two line-load solutions to obtain the normalized detachment distance, defined as Δ ≡ 2μδ/N, by

• Download figure
• Open in new tab
• Download PowerPoint

In figure 8 we plot the detachment distance Δ versus T_n in a log–log plot. If the scaling argument were correct, then the curve should have a linear region with slope 1/3 (dashed line in figure 8). While there is a linear region for$T_{\text{n}}\in (1,100),$ the slope is much smaller (about 1/20) compared with 1/3. This suggests that the scaling argument may not be applied, at least for the 2D problem.

• Download figure
• Open in new tab
• Download PowerPoint

Several reasons may account for the discrepancy between their scaling arguments and our model:

(1)	The scaling for the flow velocity given by (5.1b) is difficult to justify. Darcy's law states that the flux is proportional to the pore pressure gradient, $\mathrm{\nabla }p$. However, the ‘pore pressure’ used in equations (5.1a,b) and (5.2) is the Laplace pressure pL induced by the surface stress. The Laplace pressure is defined on the surface only, while the pore pressure is defined in the bulk—they are completely different quantities. Specifically, the Laplace pressure is the traction acting on the surface (it is proportional to the mean surface curvature). On the other hand, the pore pressure exists everywhere inside the gel, and its gradient controls liquid flow. On the surface, both the elastic stress in the network and the pore pressure contribute to the Laplace pressure. Therefore, p ≠ pL unless the stress acting on the network is zero which is not true in general. Indeed, for our case, the pore pressure p on the surface is identically zero at all times due to the permeable boundary condition. On the other hand, the Laplace pressure pL is non-trivial as long as the surface is not flat. This means that the Laplace pressure is exerted by elasticity of the network.
(2)	The assumption that $p^{\text{L}}\sim \sigma /\delta$ is difficult to justify.
(3)	To compute the detachment distance, we need a reference point serving as a datum (here the reference point is at Xn = 20). This is an inevitable feature of 2D deformation of an elastic half space. If we pick a different reference point, Δ will be modified by a constant. Therefore, the log–log plot of Δ versus Tn will be different. For instance, let us pick a new reference point extremely far away and accordingly we add a big number C to Δ, i.e. the new detachment distance ${\mathrm{\Delta }}^{\text{new}}=\mathrm{\Delta }+C,C\gg \mathrm{\Delta }.$ In this case, ${\mathrm{\Delta }}^{\text{new}}$ will not change with time practically—the curve will be just a flat line.

Finally, it should be noted that both our model and the scaling arguments presented by Berman et al. [46] are based on linear poroelastic theory. However, the deformation observed in their experiments is quite large. Also, our history of peak formation is over-simplified. In experiments, peak formation involves adhesive contact and contact line slip, which is difficult to model. All these factors could cause mismatch between experiments and theories.

It is possible to extend our plane strain analysis to an axisymmetric one using Hankel transform. The essential techniques have been demonstrated by McNamee & Gibson [14]. Nonlinear finite-element method (FEM) can be used to study the effect of large deformation on solvent flow. For example, Zhang et al. [60], Chester et al. [61] and Bouklas et al. [62] recently developed a nonlinear, transient FEM to study coupled solvent diffusion and large deformation. One can supplement this model with a surface finite-element model to study nonlinear solvent flow and large deformation. Finally, future work needs to consider more realistic boundary conditions which allow for non-equilibrium drying of the gel.

6. Summary

We study a problem of a suddenly applied line load on a linear poroelastic half space. The surface of the half space can resist deformation since it has a constant surface stress and bending stiffness. We solve this problem using transform methods. We show how the displacement, stress and flow fields are affected by surface stress. We find that surface stress reduces the pore pressure gradient, thus reducing flow near the line load. The stress and strain fields are regularized by the presence of surface bending. However, in our simulations, we deliberately make this a small effect to highlight the role of surface stress. For this case, the line force is mostly supported by the surface stress, forming the well-known Neumann's triangle [56]. Except at distances very close to the line load, the normal stress σ_zz on the surface relaxes as the fluid drains out, while the lateral component σ_xx increases with time. We use our solution to understand the scaling argument by Berman et al. [46] and find discrepancies which require future work.

Data accessibility

This work does not have any experimental data. The Jupyter code for numerical results is provided in the electronic supplementary material.

Authors' contributions

Z.L. and C.-Y.H. formulated the problem and drafted the initial manuscript. Z.L. carried out the numerical calculations in consultation with N.B. and C.-Y.H. N.B. helped draft the manuscript and checked the formulae. All authors gave final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

C.-Y.H. is supported by National Science Foundation, USA MoMS program under grant no. 1903308.

Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.4804197.