The role of a pre-wetting layer in resolving the contact line paradox of thin films in infinite elastic domains

Abstract

Hydraulic fracturing for oil and gas production from shale formations, as well as natural geological phenomena, involve the propagation of thin viscous films within elastic media. For viscous fluids, stress diverges as the thickness of the film tends to zero, arresting the propagation of the film, and thus implying the contact line paradox. For free-surface films, this paradox is resolved by considering a precursor film, leading to Tanner's law. This approach was extended recently for viscous films between a thin elastic plate and a rigid solid, allowing calculation of the film propagation rate. In this work, we examine the effect of a pre-wetting layer on the rate of propagation of a viscous flow within an infinitely deep and long domain. We analyse the linear and nonlinear dynamic problems, and perform a self-similarity analysis. We find that peeling front propagation scales as time to the power of $1/9$ and $1/3$ for thin and thick pre-wetting layer limits, respectively. Our results contribute to the understanding of the contact line paradox in elastic media and the crucial role of the pre-wetting layer in resolving it.

1. Introduction

Natural geological phenomena, such as the formation of laccoliths by the lateral flow of lava beneath an elastic sediment layer (Michaut 2011; Lister, Peng & Neufeld 2013; Pihler-Puzović et al. 2015) and gravity-driven surface lava flows under solidified crusts (Pihler-Puzović et al. 2015), as well as human-made applications, such as hydraulic fracturing for production of oil and gas from shale formations that releases fluid waste (Dana et al. 2018), formation and growth of blisters (Pihler-Puzović et al. 2015), manufacture of flexible electronics and micro-electromechanical systems, the reopening of airways, and the suppression of viscous fingering in a deformable Hele-Shaw cell (Lister et al. 2013), involve viscous spreading of fluids beneath elastic surfaces.

It is known that in surface-tension-driven problems, in the limit of length scales smaller than the capillary length, an assumption that the thickness of a droplet tends to zero at the contact line leads to divergent viscous stresses, and thus to the arresting of the contact line propagation (Huh & Scriven 1971). This apparent paradox, which conflicts with everyday experience of spreading droplets, was resolved by considering the development of a thin precursor film due to intermolecular interactions (e.g. van der Waals) in advance of the contact line (de Gennes 1985). Thus a local balance between viscous dissipation and the rate of change of surface energy gives rise to Tanner's law (Tanner 1979; Bonn et al. 2009), in which the droplet radius advances with speed ${\rm d}R/{\rm d}t\propto \theta ^3$ for apparent contact angle $\theta$, thus using droplet volume conservation, one obtains readily that $R$ increases as $t^{1/10}$ (de Gennes, Brochard-Wyart & Quéré 2004).

This approach was extended recently by employing the pre-wetting layer to study a propagation of a viscous flow between elastic plates and a rigid solid. For example, Lister et al. (2013) examined an elastic peeling problem in the distinct limits of peeling by bending and peeling by pulling, and applied their results to the radial spread of a fluid blister over a thin pre-wetting film. Pihler-Puzović et al. (2015) investigated single- and two-phase displacement flows in which the injection of fluid is accommodated by the inflation of the sheet and the outward propagation of an axisymmetric front beyond which the cell remains approximately undeformed. Peng & Lister (2020) investigated the spreading of viscous fluid injected under an elastic sheet, where initially the viscous liquid fills the narrow gap, and is driven by gravity and by elastic bending and tension forces, resisted by viscous forces. By identifying all of the possible asymptotic combinations, they revealed a rich variety of different behaviours. Hewitt, Balmforth & De-Bruyn (2015) considered a nonlinear diffusion equation describing the planar spreading of a viscous fluid injected between an elastic sheet and an underlying rigid plane, where the dynamics was assumed to depend on the physical conditions at the contact line where the sheet is lifted off the plane by the fluid. Bunger & Detournay (2007) proposed a small-time asymptotic solution for a penny-shaped fluid-driven fracture, which was obtained semi-analytically. In particular, they concluded that the similarity solution is unusual since the two length scales of the problem – the radius of the fracture and the radius of the fluid front – evolve according to two different power laws of time.

Other works concentrated on the propagation of cracks in infinite domains, in the two limits of a thin or thick elastic plate. For example, Lai et al. (2015) performed an experimental investigation of a fluid-driven crack in a gelatin matrix, and verified the influence of different experimental parameters such as the injection flow rate, Young's modulus of the matrix, and fluid viscosity on the shape of crack. Different governing physical mechanisms of the peeling process were also studied extensively. Considerable work was done taking into account the adhesive interaction between the solid and the fluid (McEwan & Taylor 1966; Hosoi & Mahadevan 2004). Some studies concentrated on investigating systems with gravitation-governed peeling processes (Buckmaster 1977; Huppert 1982; Momoniat 2006; Howell, Robinson & Stone 2013; Balmforth, Craster & Hewitt 2015; Hewitt et al. 2015).

In this work, we examine the effect of a pre-wetting layer on the rate of propagation of a viscous flow within an infinitely deep and long domain. Our model, which aims to explain fluid dynamics that occurs when a crack already exists, is commonly denoted in the context of hydraulic fracturing as ‘flowback’ or ‘produced water’ (see e.g. Dana et al. 2018). Specifically, our model problem is governed by the equations of elasticity, where the lubrication theory is used to describe the fluid dynamics in the whole region (close to the peeling front and away from it). In this sense, considerable work was done on infinite elastic spaces and half-spaces in the context of hydraulic fracturing. Various models of hydraulically driven crack propagation were suggested (Spence & Sharp 1985; Lister 1990; Savitski & Detournay 2002; Lai et al. 2015, 2016a,b). Most of the analytic works implement known solutions for general axisymmetric problems of stress functions, e.g. Love's bi-harmonic equations (Sneddon 1995). The integral relation between the fluid film pressure and the surface deformation, derived by Spence & Sharp (1985), together with Reynolds equation obtained from the fluid analysis, completes the definition of the integro-differential problem and allows us to gain insight into the behaviour of the system under consideration. More specifically, in § 2, we present the model problem, and in § 3, we discuss the self-similarity solutions in the two limits of thick and thin pre-wetting layers, where a thin pre-wetting layer represents a precursor film. In § 4, we discuss our results, including a comparison between the solutions for the dynamic linear and nonlinear problems, and the corresponding solutions obtained by self-similarity analysis in these two limits. We conclude our findings in § 5, give the details of our numerical algorithms in the various cases in Appendix A, and present additional details and comparisons between our results and the results known from the literature in Appendix B.

2. Problem formulation

In this work, we consider the configurations that are shown in figure 1, where a fluid layer with initial thickness $h_0$ is situated between infinitely deep and long elastic and rigid solids, or between two elastic solids. Hereafter, we will focus only on the configuration shown in figure 1(a), since the two configurations are equivalent. The configuration in figure 1 is assumed to be axisymmetric, where the radial coordinate, denoted by $r$, is in the plane of the solids’ surfaces, and the $z$-axis is traversing the solids perpendicularly. An inlet is located at $r=0$, which allows the introduction of a controlled fluid flux into the film, and the peeling front is located at $r_F(t)$. The fluid is a non-compressible Newtonian fluid with viscosity $\mu$, and the solid is a Hookean isotropic material with shear modulus $G$ and Poisson's ratio $\nu$.

Figure 1. Illustrations of the cross-sectional view of the considered system, where a thin layer of fluid is situated (a) between a rigid solid and a semi-infinite elastic domain, and (b) in the centre of an infinite elastic domain. The two configurations are equivalent since if $k_1$ and $k_2$ are the material constants of $h_1$ and $h_2$, respectively, defining $\tilde {h}=h_1+h_2$ and the effective material constant $k$ to be $\tilde {k}=1/(1/k_1+1/k_2)$, we get in (b) the same problem formulation as in (a), but with $\tilde {h}$ and $\tilde {k}$ instead of $h$ and $k$, respectively. For simplicity, hereafter we focus on configuration (a).

We denote the overall thickness of the fluid by $h$. Using the lubrication assumption – namely, assuming a small ratio between the characteristic film thickness and the length scale ($h^*/r^*\ll 1$) – we get the Reynolds equation, describing the dynamics of the fluid,

(2.1)\begin{equation} \frac{\partial h}{\partial t}=\frac{1}{12\mu}\,\frac{1}{r}\,\frac{\partial}{\partial r}\left( rh^3\,\frac{\partial p}{\partial r} \right), \quad (r,t)\in(0,r^*)\times(0,T_{fin}), \end{equation}

where our assumption that there exists a pre-wetting layer of thickness $0< h_0\ll \max \{h(r,t)\}$ implies that there exists a function $d(r,t)$, so that

(2.2)\begin{equation} h(r,t)=h_0+d(r,t). \end{equation}

In addition, in (2.1), the following notations are used: $p(r,t)$ denotes the pressure of the fluid, $r^*\gg 1$ is the ‘far’ boundary of our domain, and $T_{fin}$ is some final time. We neglect the shear stress exerted by the fluid on the solid, since the ratio between the shear stress and the pressure is very small due to the lubrication assumption: $\tau /p\approx (d^*/r^*)^2\ll 1$, where $d^*$ can be defined as $d^*:=\max _{r\in (0,r^*)}\{d(r,0)\}$.

Moreover, we use the general relation between the deformation of the surface, $d$, and the fluid pressure, $p$, as derived by Spence & Sharp (1985), which may be expressed as

(2.3)\begin{equation} p(r,t)=-k\,\int_0^\infty M\left(\frac{r}{s}\right)\frac{\partial}{\partial s}d(s,t)\,\frac{{\rm d}s}{s}, \end{equation}

where

(2.4)\begin{equation} k=\frac{G}{1-\nu}=\frac{E}{2(1-\nu^2)} \end{equation}

is a material constant, and $M(m)$ is a kernel function that is given by

(2.5)\begin{equation} m\,M(m) = \begin{cases} \dfrac{2m\,E(m)}{{\rm \pi}(1-m^2)}, & m<1, \\[10pt] \dfrac{2m^2\,E(1/m)}{{\rm \pi}(1-m^2)} + K\left(\dfrac{1}{m}\right), & m>1. \end{cases} \end{equation}

In (2.5), $K(m)$ and $E(m)$ denote the complete elliptic integrals of the first and second kind, respectively.

For the boundary conditions, which depend on the specific case (linear problem, nonlinear problem, self-similar problem in the linear case, or self-similar problem in the nonlinear case), see § 3.

To render the problem dimensionless, we employ the transformations

(2.6)\begin{align} \begin{gathered} R=\frac{r}{r^*}, R_F=\frac{r_F}{r^*}, T=\frac{t}{t^*}, H(R,T)=\frac{h(r,t)}{h^*}, D(R,T)=\frac{d(r,t)}{d^*}, P(R,T)=\frac{p(r,t)}{p^*}, \end{gathered} \end{align}

where $r_F=\min _{r\in (0,r^*)}\{d(r,t)\}$ denotes the fluid front at time $t$. Lowercase letters, uppercase letters and superscript asterisks denote dimensional variables, dimensionless variables and characteristic values, respectively. Requiring that the dimensionless problem will be free of dimensionless numbers implies the following relations between the characteristic values:

(2.7a,b)\begin{equation} p^*=\frac{k d^*}{r^*} \quad \text{and} \quad t^*=\frac{\mu(r^*)^{3}}{k (h^*)^2\,d^*}. \end{equation}

In our numerical computations for the dynamic problem, ‘$\infty$’ in the upper boundary of the integral in (2.3) will be replaced by the finite characteristic value $r^*$, which is assumed to satisfy $r^*\gg r_F(T_{fin})$, which is sufficiently large, so that the pressure gradients vanish for all $r>r^*$ and thus the thickness $h(r,t)$ is constant for all $r>r^*$ and $t\leqslant T_{fin}$. Thus, after rendering the problem dimensionless, the upper boundary of (2.3), which is scaled by $r^*$, becomes 1. Moreover, in the self-similar cases, which will be discussed in the next section, the upper boundary of the integral in (2.3) is the contact line position $R_F$, which serves also as the rescaling factor for obtaining the self-similar variable $\eta$, namely $\eta =r/R_F$. Thus also in the case of self-similar solutions, the upper boundary of this integral in dimensionless form is equal to 1. This allows us to avoid singularity in computations by following Peck et al. (2018), so that in the continuation of this study, we will use the following inverted version of (2.3) in dimensionless terms:

(2.8)\begin{equation} D(R,T)=\frac{8}{\rm \pi}\int_0^1 \frac{\partial P}{\partial Y}\,\mathcal{K}(Y,R)\,{\rm d}Y+\frac{8\sqrt{1-R^2}}{\rm \pi}\int_0^1\frac{Y\,P(Y,T)}{\sqrt{1-Y^2}}\,{\rm d}Y, \end{equation}

where the kernel $\mathcal {K}(Y,R)$ is given by

(2.9)\begin{equation} \mathcal{K}(Y,R)=Y\left[E\left(\arcsin(Y) \left|\frac{R^2}{Y^2}\right.\right)-E\left(\arcsin(\chi)\left|\frac{R^2}{Y^2}\right.\right)\right], \end{equation}

with $\chi :=\min (1,Y/R)$, and where $E(\phi |m)=\int _0^{\phi }\sqrt {1-m\sin ^2(\theta )}\,\textrm {d}\theta$ denotes the incomplete elliptic integral of the second kind. Note that contrary to the kernel $M$, defined in (2.5), the kernel $\mathcal {K}$ is not singular for $(Y,R)\in [0,1]^2$ (Peck et al. 2018).

We present our numerical algorithm for the solution of this dynamic problem in the linear and nonlinear cases in §§ A.1 and A.3, respectively. In § 3, we present the self-similar solution for both relatively thick and relatively thin pre-wetting layer limits, and then in § 4, we compare the results of our dynamic simulations with the corresponding self-similar solution in each of the two limits. Moreover, in § 4, we show and discuss our results for the contribution of the pre-wetting layer thickness to the peeling front propagation rate.

3. Self-similarity analysis

The pre-wetting layer thickness, namely the magnitude of $h_0$ in $h(r,t)=h_0+d(r,t)$, has two distinct limits that lead to different regimes of the solution. If the pre-wetting layer thickness $h_0$ is much greater than the characteristic deformation, then the equation in (2.1) can be approximated at leading order by a linear partial differential equation (PDE). Otherwise, if the pre-wetting layer thickness $h_0$ is much smaller than the characteristic deformation, then the equation in (2.1) becomes strongly nonlinear. Note that in a general case of the equation in (2.1), the answer to the question of whether a similarity solution does or does not exist depends on the differential operator in the pressure $p$. For tension or bending ($p=-h_{xx}$ or $p=h_{xxxx}$), there is no similarity solution to connect to the pre-wetted film (Hewitt et al. 2015). Indeed, tending $h_0$ to zero in these cases results in a singular limit, so that in order to get an advancing contact line, a regularisation, such as a pre-wetted film, is required (Flitton & King 2004). On the other hand, in the case of gravity ($p=h$), a similarity solution exists (Huppert 1982). Moreover, also for an advancing contact line in an elastic half-space, it has been shown previously that solutions converge as $h_0$ tends to zero (Touvet et al. 2011), so that a similarity solution exists.

In this work, we investigate the dependence of the peeling front velocity on the thickness of the pre-wetting layer, $h_0$, where as we will discuss in this section, the self-similar analysis provides us with the power law of time for the peeling front propagation, with the different powers for both these limits. In this section, both limits will be analysed separately. The self-similar solutions that will be obtained in the current section will be found in the range $R\in (0,R_F(T))$.

3.1. Linear case – obtained for $d^*\ll h_0$

In this case, let us set the characteristic film thickness $h^*$ to be the thickness of the pre-wetting layer, namely $h^*=h_0$, and let us denote by $\varepsilon$ the ratio

(3.1)\begin{equation} \varepsilon=\frac{d^*}{h_0}\ll 1. \end{equation}

Using (2.6) to transform dimensional variables to the dimensionless notation, and substituting this ratio into the definition of $H$, we get, in the case of a relatively large pre-wetting layer thickness, that

(3.2)\begin{equation} H(R,T)=1+\varepsilon\,D(R,T). \end{equation}

It is easy to verify that substituting the dimensionless variables in (2.6) into (2.1) and (2.3), and using the definitions in (2.7a,b) as well as the assumption that $\varepsilon \ll 1$, we get at the leading order with respect to $\varepsilon$ the following linear equation for $D$:

(3.3a)$$\begin{gather} \frac{\partial D}{\partial T}=\frac{1}{12 R}\,\frac{\partial}{\partial R}\left( R\,\frac{\partial P}{\partial R}\right), \quad R\in(0,1), \end{gather}$$

(3.3b)$$\begin{gather}P(R,T)=-\int_0^{1}{ M\left(\frac{R}{S}\right)\frac{\partial D}{\partial S}\,\frac{{\rm d}S}{S} }. \end{gather}$$

These equations are solved subject to the following boundary and initial conditions:

(3.4a)$$\begin{gather} \left.\frac{\partial P}{\partial R}\right|_{R=0}=0, \quad T\in(0,T_{fin}), \end{gather}$$

(3.4b)$$\begin{gather}P(R=1,T)=0, \quad T\in(0,T_{fin}), \end{gather}$$

(3.4c)$$\begin{gather}D(R,T=0)=\left\{\begin{array}{@{}ll} H_{max}(1-q^2 R^2), & \text{for $R\in[0,1/q],$}\\[2pt] 0, & \text{for $R\in(1/q,1],$} \end{array}\right. \end{gather}$$

where $H_{max}$ and $q\gg 1$ are some prescribed constants. Note that $R=1$ is located far from the peeling front, so that the interface is flat there, and the conditions $P(R=1,T)=0$, $\partial P/\partial R|_{R=1}=0$, $D(R=1,T)=0$, and $\partial D/\partial R|_{R=1}=0$ are equivalent. We prescribed $P(R=1,T)=0$ in (3.4b) for numerical convenience. Moreover, note that the problem in (3.3) satisfies volume conservation, namely, there exists $\Delta V>0$ that is independent of time such that

(3.5)\begin{equation} 2{\rm \pi}\int_0^{R_F(T)} D(R,T)\,R\,{\rm d}R\approx\Delta V. \end{equation}

For details regarding the numerical algorithm for solving the dynamic problem in (3.3)–(3.4), see § A.1.

In order to obtain a self-similar formulation for this problem up to the boundary $R=R_F(T)$, rather than up to $R=1$ (which is based on an assumption that the film thickness is undisturbed for $R>R_F(T)$ and implies in particular that the upper boundary of the integral in (3.3b) is also replaced by $R_F(T)$), let us follow Spence & Sharp (1985) and assume that the velocity of the front propagation is proportional to a power law of time. More specifically, let us assume that

(3.6)\begin{equation} R_{F}(T)=L T^{\alpha}, \end{equation}

where $L$ and $\alpha$ are some constants. Furthermore, reasonable candidates for the self-similarity variable and the solution are

(3.7a–c)\begin{equation} \eta=L^{-1}RT^{-\alpha}, \quad P=L^{3} T^\beta\,\tilde{P}(\eta), \quad D=L T^\gamma\,\tilde{D}(\eta), \end{equation}

where $\tilde {P}(\eta )$ and $\tilde {D}(\eta )$ are self-similar functions to be determined. Substituting these transformations into (3.3), and requiring the elimination of explicit time dependency in the integro-differential system, we obtain the following system of equations for $\alpha$, $\beta$ and $\gamma$:

(3.8)\begin{equation} \left.\begin{gathered} \beta-2\alpha-\gamma+1=0,\\ \gamma-\alpha-\beta=0,\\ \gamma+2\alpha=0, \end{gathered}\right\} \end{equation}

which yields the powers

(3.9a–c)\begin{equation} \alpha= \tfrac{1}{3}, \quad \beta=-1, \quad \gamma=-\tfrac{2}{3}. \end{equation}

Thus we obtain the following integro-differential problem for $\tilde {P}$ and $\tilde {D}$:

(3.10a)$$\begin{gather} 8\tilde{D}+4\eta\,\frac{{\rm d}\tilde{D}}{{\rm d}\eta}+\frac{1}{\eta}\,\frac{{\rm d}}{{\rm d}\eta} \left(\eta\,\frac{{\rm d}\tilde{P}}{{\rm d}\eta}\right)=0, \end{gather}$$

(3.10b)$$\begin{gather}\tilde{P}(\eta)=-L^{-3}\int_0^1{M\left(\frac{\eta}{\xi}\right)\frac{{\rm d}\tilde{D}}{{\rm d}\xi}\, \frac{{\rm d}\xi}{\xi}}, \end{gather}$$

where $\xi =L^{-1}ST^{-\alpha }$. As mentioned previously, for computational purposes, we will replace (3.10b) by its inverted form

(3.10c)\begin{equation} \tilde{D}(\eta)=L^3\left[\frac{8}{\rm \pi}\int_0^1 \frac{\partial \tilde{P}}{\partial \xi}\,\mathcal{K}(\xi,\eta)\,{\rm d}\xi+\frac{8\sqrt{1-\eta^2}}{\rm \pi}\int_0^1\frac{\xi\, \tilde{P}(\xi)}{\sqrt{1-\xi^2}}\,{\rm d}\xi\right], \end{equation}

where the kernel $\mathcal {K}$ is given in (2.9).

The equations in (3.10a) and (3.10c) are solved subject to the boundary conditions and constraints

(3.11a)$$\begin{gather} \left. \frac{{\rm d}\tilde{D}}{{\rm d}\eta}\right|_{\eta=0}=0, \end{gather}$$

(3.11b)$$\begin{gather}\left.\frac{{\rm d}\tilde{D}}{{\rm d}\eta}\right|_{\eta=1}=0, \end{gather}$$

(3.11c)$$\begin{gather}L=\begin{pmatrix}\displaystyle\frac{\Delta V}{\displaystyle 2{\rm \pi}\int_{0}^1 \tilde{D}\eta \,{\rm d}\eta}\end{pmatrix}^{1/3}, \end{gather}$$

where $\Delta V$ is some given volume, and in order to compare the self-similar solution with the solution to the dynamic problem, it is chosen according to (3.5) (calculated with e.g. the initial condition for the dynamic problem). Note that although the boundary condition in (3.11a) differs from (3.4a), they are equivalent since both of them are obtained from the corresponding integral equations if sufficient regularity of the pressure and the thickness are assumed. This follows from the observation that the kernels $M(R/S)$ and $\mathcal {K}(Y,R)$, which were defined in (2.5) and (2.9), respectively, satisfy $\textrm {d}\mathcal {K}/\textrm {d}R|_{R=0}=\textrm {d}M/\textrm {d}R|_{R=0}=0$. Moreover, it can be observed that the condition in (3.11b) differs from the condition in (3.4b). The reason for this discrepancy is that the domains for the dynamic and self-similar problems are different. While in the dynamic problem the end of the domain is sufficiently far from the centre of symmetry, where the pre-wetting film remains unperturbed, the end of the domain in the self-similar problem is at the peeling front itself (which in the case $h_0\gg d$ is a region, rather than a specific point), so that the condition in (3.11b) follows from the definition of the position of the peeling front. More specifically, we found that the self-similar solution oscillates around $0$ as it reduces in magnitude, and we thus chose that first location in which $\partial \tilde {D}/\partial \eta = 0$ as the location of the front. The power-law relations are not affected by this choice, and the magnitude of $\tilde D$ after this point is approximately $1\,\%$ of the maximal value.

For details regarding the numerical algorithm for finding the solution for the self-similar problem in (3.10)–(3.11), see § A.2.

3.2. Nonlinear case – obtained for $d^*\gg h_0$

In the limit as $h_0\rightarrow 0$, the viscous resistance to flow in the thin pre-wetting region tends to infinity, thus the characteristic time scale diverges as well. This inhibits information from passing beyond the peeling front. These solutions, often referred to as compactly supported solutions, are common to these peeling problems (see discussion in Lister et al. 2013).

In this case, let us set the characteristic film thickness $h^*$ to be $d^*$, and let us define

(3.12)\begin{equation} H_0=\frac{h_0}{d^*}\ll 1, \end{equation}

so that now we get that

(3.13)\begin{equation} H(R,T)=H_0+D(R,T). \end{equation}

In this case, substituting the dimensionless variables in (2.6) into (2.1) and (2.3), and using the definitions in (2.7a,b) as well as the assumption that $H_0 \ll 1$, we get at the leading order with respect to $H_0$ the following nonlinear PDE for $D$:

(3.14a)\begin{equation} \frac{\partial D}{\partial T}=\frac{1}{12 R}\,\frac{\partial}{\partial R}\left( RD^3\,\frac{\partial P}{\partial R}\right), \quad R\in(0,1), \end{equation}

which is accompanied as previously by the integral equation

(3.14b)\begin{equation} P(R,T)=-\int_0^{1} M\left(\frac{R}{S}\right)\frac{\partial D}{\partial S}\,\frac{{\rm d}S}{S} . \end{equation}

We solve the problem in (3.14) subject to the following boundary and initial conditions:

(3.15a)$$\begin{gather} \left.\frac{\partial P}{\partial R}\right|_{R=0}=0, \quad T\in(0,T_{fin}), \end{gather}$$

(3.15b)$$\begin{gather}\left.\frac{\partial P}{\partial R}\right|_{R=1}=0, \quad T\in(0,T_{fin}), \end{gather}$$

(3.15c)$$\begin{gather}D(R,T=0)=\begin{cases} (H_{max}-H_0)(1-q^2 R^2), & \text{for $R\in[0,1/q],$}\\ 0, & \text{for $R\in(1/q,1],$} \end{cases} \end{gather}$$

where, as in the linear case, $H_{max}$ and $q\gg 1$ are some prescribed constants, and $H_0$ denotes the thickness of the pre-wetting layer whose influence on the power-law dynamics of the peeling front we are interested to investigate. Although the conditions in (3.15) look somewhat different to the conditions in (3.4), they are equivalent, and this difference is only for numerical purposes. Specifically, as stated previously, $P(R=1,T)=\textrm {d}P/\textrm {d}R|_{R=1}=0$, so we may chose any of them according to numerical convenience.

In addition, similarly to the linear case, the volume is conserved also in the nonlinear problem, namely, there exists $\Delta V>0$ such that (3.5) is satisfied. For further details regarding the numerical solution of (3.14) subject to (3.15), see § A.3.

In order to obtain a self-similar solution for the problem in (3.14) up to the boundary $R=R_F(T)$ (so that also the upper boundary of the integral in (3.14b) is $R_F(T)$), we use again the transformation in (3.6) and define

(3.16a–c)\begin{equation} \eta=L^{-1}RT^{-\alpha}, \quad P=T^\beta\,\tilde{P}(\eta), \quad D=LT^\gamma\,\tilde{D}(\eta), \end{equation}

where, in order to eliminate the time dependency in the nonlinear case, we get that $\alpha$, $\beta$ and $\gamma$ should satisfy the system

(3.17)\begin{equation} \left.\begin{gathered} 2\gamma-2\alpha+\beta+1=0,\\ \gamma-\alpha-\beta=0,\\ \gamma+2\alpha=0, \end{gathered} \right\}\end{equation}

which yields the powers

(3.18a–c)\begin{equation} \alpha= \tfrac{1}{9}, \quad \beta=-\tfrac{1}{3}, \quad \gamma=-\tfrac{2}{9}. \end{equation}

Thus in this case, we obtain the following nonlinear integro-differential problem for $\tilde {P}$ and $\tilde {D}$:

(3.19a)$$\begin{gather} 8\tilde{D}+4\eta\,\frac{{\rm d}\tilde{D}}{{\rm d}\eta}+\frac{3}{\eta}\,\frac{{\rm d}}{{\rm d}\eta} \left(\eta\tilde{D}^3\,\frac{{\rm d}\tilde{P}}{{\rm d}\eta}\right)=0, \end{gather}$$

(3.19b)$$\begin{gather}\tilde{P}(\eta)=-\int_0^1 M\left(\frac{\eta}{\xi}\right)\frac{{\rm d}\tilde{D}}{{\rm d}\xi}\,\frac{{\rm d}\xi}{\xi}, \end{gather}$$

or equivalently,

(3.19c)\begin{equation} \tilde{D}(\eta)=\frac{8}{\rm \pi}\int_0^1 \frac{\partial \tilde{P}}{\partial \xi}\,\mathcal{K}(\xi,\eta)\,{\rm d}\xi+\frac{8\sqrt{1-\eta^2}}{\rm \pi}\int_0^1\frac{\xi\, \tilde{P}(\xi)}{\sqrt{1-\xi^2}}\,{\rm d}\xi, \end{equation}

where $\xi =L^{-1}ST^{-\alpha }$.

The equations in (3.19) are solved subject to the boundary conditions

(3.20a)$$\begin{gather} \left.\frac{{\rm d}\tilde{D}}{{\rm d}\eta}\right|_{\eta=0}=0, \end{gather}$$

(3.20b)$$\begin{gather}\tilde{D}(1)=0. \end{gather}$$

Note that the condition in (3.20a) is exactly as in the linear case, namely in (3.11a) (and contrary to the nonlinear dynamic problem, we do not need to prescribe here explicitly the condition for $\textrm {d}P/\textrm {d}\eta |_{\eta =0}$, since we use a different methodology to solve this problem, rather than finite differences.) Moreover, similarly to our discussion in the linear case, the condition in (3.20b) is prescribed at the peeling front, thus it differs from the condition in (3.15b), which is prescribed at the far end of the domain. However, note that the condition in (3.20b) also differs from the boundary condition in (3.11b), which was prescribed at the peeling front in the linear case. The reason for this difference is that in the nonlinear case, the function $\tilde {D}$ is monotone decreasing and the derivative $\textrm {d}\tilde {D}/\textrm {d}\eta$ diverges at $\eta =1$ (and attains a global minimum at the peeling front, rather than a local one, as in the linear case). Since we cannot prescribe now that $\textrm {d}\tilde {D}/\textrm {d}\eta |_{\eta =1}=0$, we use an alternative definition of the peeling front, which is that the thickness of the film must vanish at $\eta =1$. Note that in the linear case, $\tilde {D}(\eta )$ is not a monotone decreasing function, but oscillates around the $\eta$-axis. Thus we cannot prescribe the boundary condition $\tilde {D}(1)=0$ in the linear case, since it is ambiguous.

The volume conservation constraint in the nonlinear case expressed in terms of self-similar variables yields the following formula for the length scale $L$:

(3.21)\begin{equation} L=\begin{pmatrix}\displaystyle\frac{\Delta V}{\displaystyle 2{\rm \pi}\int_{0}^1 \tilde{D}\eta \,{\rm d}\eta}\end{pmatrix}^{1/3}. \end{equation}

For details regarding the solution of the self-similar problem in this case, see § A.4.

4. Results

In figure 2, we show a comparison between the dynamic solution for the linear problem (under the assumption that $h_0\gg d^*$) and the corresponding self-similar solution, according to the discussion in § 3.1. Specifically, in figures 2(a) and 2(c), we show the film thickness $D$ and the pressure $P$, respectively, versus $R$ for various time instances, where the samples of the solution at different time instances obtained by solving the linear dynamic problem are marked by solid lines, and the samples of the solution at the same time instances obtained by solving the linear self-similar problem and transferring the result to the dynamic solution are marked by dashed lines of the corresponding colours. In figures 2(b) and 2(d), we show the film thickness $\tilde {D}$ and the pressure $\tilde {P}$, respectively, versus $\eta$ for various time instances, where the self-similar solution is marked by a dashed black line, and the coloured curves correspond to the samples of the dynamic solution (at the same time instances as in figures 2a,c) that were transferred to self-similar variables. It can be seen that although initially there is a large difference between the dynamic and self-similar solutions, as the time increases, the dynamic solution converges to the self-similar one. In particular, this confirms that at the limit of a thick pre-wetting layer, the front propagation behaves as $R_F(T)\sim L T^{1/3}$, where $L$ is a constant. For an additional comparison between the dynamic and self-similar solutions in the linear case, see figure 5(a) in Appendix B.

Figure 2. Linear case where the initial conditions are as prescribed in (3.4c), with $H_{max}=1$ and $q=400$. (a) The film thickness $D$ versus $R$ for different time instances $T\in [5\times 10^{-6},0.001]$. (b) The self-similar function $\tilde {D}$ versus $\eta$. (c) The pressure $P$ versus $R$. (d) The self-similar function $\tilde {P}$ versus $\eta$. In the insets of (a,c), we show an increased view of the same graphs. The legend in (d) represents the time instances for all of the plots, where the solid lines represent the samples of the solution of the dynamic linear problem (in (a,c) and their transformation to self-similar variables in (b,d)), the dashed lines of the corresponding colour represent the samples of the self-similar solution transferred to the corresponding time levels (in (a,c)), and the black dashed lines (in b,d) show the solution of the self-similar problem. To transform the dynamic solution to self-similar variables, we used $L\approx 1.1987$, which was obtained by using (3.11c) with $\Delta V\approx 1.1027\times 10^{-5}$.

In figure 3, we show a comparison between the dynamic solution for the nonlinear problem (under the assumption that $h_0\ll d^*$) and the corresponding self-similar solution, according to the discussion in § 3.2. Specifically, in figures 3(a) and 3(c), we show the film thickness $D$ and the pressure $P$, respectively, versus $R$ for various time instances, where the samples of the solution at different time instances obtained by solving the nonlinear dynamic problem are marked by solid lines, and the samples of the solution at the same time instances obtained by solving the nonlinear self-similar problem and transferring the result to the dynamic solution are marked by dashed lines of the corresponding colours. In figures 3(b) and 3(d), we show the film thickness $\tilde {D}$ and the pressure $\tilde {P}$, respectively, versus $\eta$ for various time instances, where the self-similar solution is marked by a dashed black line, and the coloured curves correspond to the samples of the dynamic solution (at the same time instances as in figures 3a,c) that were transferred to self-similar variables. It can be seen that although initially there is a large difference between the dynamic and self-similar solutions, as the time increases, the dynamic solution becomes very close to the self-similar one. Note that the agreement (in figure 3d) for the pressure near the contact line is worse for $T=1.6\times 10^{-4}$ than for $T=10^{-5}$ because the pre-wetting layer thickness is approximately 1 % of the maximal film thickness at time $T=10^{-5}$, whereas at time $T=1.6\times 10^{-4}$, it increases up to 2 %. This discrepancy between the pressure obtained via a dynamical simulation and by self-similar solution is larger near the peeling front, since the ratio of the pre-wetting layer thickness to the film thickness increases even more in the vicinity of the peeling front, so that the dynamic solution in this region gets closer to the ‘thick’ pre-wetting layer limit (which is equivalent to a linear problem). In particular, the results shown in figure 3 confirm that at the limit of a very thin pre-wetting layer, the front propagation behaves as $R_F(T)\sim LT^{1/9}$, where $L$ is a constant. For an additional comparison between the dynamic and self-similar solutions in the nonlinear case, as well as a comparison between the dynamic solution in the nonlinear case in the vicinity of the peeling front and the expected power-law asymptotics in that region (Spence & Sharp 1985; Desroches et al. 1994), see figures 5(b) and 8 in Appendix B.

Figure 3. Nonlinear case where the initial conditions are as prescribed in (3.15c), with $H_0=0.02$, $H_{max}=2.3$ and $q=50$. (a) The film thickness $D$ versus $R$ for different time instances $T\in [10^{-5},1.6\times 10^{-4}]$. (b) The self-similar function $\tilde {D}$ versus $\eta$. (c) The pressure $P$ versus $R$. (d) The self-similar function $\tilde {P}$ versus $\eta$. In the inset of (a), we show an increased view of the same graph. The legend in (d) represents the time instances for all of the plots, where the solid lines represent the samples of the solution of the dynamic nonlinear problem (in (a,c) and their transformation to self-similar variables in (b,d)), the dashed lines of the corresponding colour represent the samples of the self-similar solution transferred to the corresponding time levels (in (a,c)), and the black dashed lines (in b,d) show the solution of the self-similar problem. To transform the dynamic solution to self-similar variables, we used $L\approx 0.0722$, which was obtained by using (3.21) with $\Delta V\approx 1.5\times 10^{-3}$.

In figure 4, we show the dependence of the power $\alpha$ in the approximation of the front propagation, $R_F(T,H_0)\sim LT^{\alpha (H_0)}$, on the normalized pre-wetting layer thickness $H_0=h_0/d^*$. It can be seen that at the limits of thin ($H_0\to 0$) and thick ($H_0\to H(0,0)$) pre-wetting layer thickness, the power $\alpha (H_0)$ tends to the analytically expected (obtained by self-similarity analysis) powers $\alpha =1/9$ and $1/3$, respectively. Our main result is that the dependence of $\alpha$ on $H_0$ is given approximately by

(4.1)\begin{equation} \alpha(H_0)\approx A\exp{(-BH_0)}+C, \end{equation}

where $A=-0.2253$, $B=1.195$ and $C=0.3348$. In the inset of figure 4 we show the dependence of the normalized (divided by its maximal value) velocity at different time instances on the normalized (divided by its maximal value) pre-wetting layer thickness. It can be seen that the velocity of the peeling front $V:=\textrm {d}R_F/\textrm {d}T$, which is approximately a linear function of $H_0$, tends (for all examined time instances) to non-vanishing constants, as $H_0$ tends to 0. Although, the sequence of these constants constitutes a monotone decreasing function of time, in our simulation the constants do not vanish.

Figure 4. The power $\alpha$ in the approximation of the peeling front propagation $R_F(T,H_0)\sim LT^{\alpha (H_0)}$, versus $H_0$. The blue dots represent values of $\alpha$, which were obtained from fitting of $R_F(T,H_0)$ (calculated by our simulations) versus $T$ to the function $L T^{\alpha (H_0)}$ (for a series of examples of this fitting, see figures 6 and 7 in Appendix B). All of the simulations were performed with the initial conditions prescribed in (3.15c), where we have set $H_{max}=2.3$ and $q=50$. The blue curve was obtained by fitting the results obtained by simulation to the function $F=A\exp {(-BH_0)}+C$, which yielded $A=-0.2253$, $B=1.195$ and $C=0.3348$, with $\lim _{H_0\to 0}F=0.1095\approx 1/9$, $\lim _{H_0\to \infty }F=0.3348\approx 1/3$, and with the goodness of fit satisfying R-square $0.9989$. The red and black dashed lines represent the powers $1/3$ and $1/9$ (obtained by self-similarity analysis) at the two limits of a very thick and a very thin pre-wetting layer, respectively. In the insets, we show a simplified sketch of our system that should illustrate the meaning of our results, and a normalized velocity of the contact line $V/\max \{V\}$ versus the normalized thickness of the pre-wetting layer $H_0/\max \{H_0\}$, for various values of time instances $T\in [9.6\times 10^{-6},1.596\times 10^{-4}]$, where the uppermost curve corresponds to $T=9.6\times 10^{-6}$, the lowermost curve corresponds to $T=1.596\times 10^{-4}$, and the time difference between any two adjacent curves (the time corresponding to any curve minus the time corresponding to the curve just above it) is $10^{-6}$.

5. Conclusions

In this work, we studied the effect of the pre-wetting layer thickness on the propagation rate of the peeling front of fluid that is situated between two infinitely deep and long elastic solids or between an elastic solid and a rigid substrate. We solved the dynamic problem numerically in two limits, a relatively thin and a relatively thick pre-wetting layer. Moreover, we found the corresponding self-similar solution in each of these two limits, and obtained excellent agreement between the dynamic and self-similar solutions in both cases.

Our findings show that in the linear limit, the peeling front propagates at rate $t^{1/3}$, while in the nonlinear case, it propagates at rate $t^{1/9}$. When compared to the results of Lister et al. (2013), who found that the peeling front propagates at rates $t^{7/22}$ and $t^{3/8}$ (for pulling and bending limits, respectively), it is clear that the front propagation is slower for the limit of very thick elastic layers. In order to obtain a dependence of the peeling front propagation rate on the pre-wetting layer thickness, we simulated our dynamic nonlinear solver for a range of pre-wetting layer thicknesses, between the limits of a relatively thin and a relatively thick pre-wetting layer. We found that the peeling front propagation scales as time to the powers that constitute a monotone increasing function of the pre-wetting layer thickness, and accept the values within the range $[1/9,1/3]$, where $\lim _{H_0\to 0}\alpha (H_0)\approx 1/9$ and $\lim _{H_0\to H(0,0)}\alpha (H_0)\approx 1/3$, exactly as expected by the self-similarity analysis.