2. Model formulation

We consider a Newtonian liquid film coating on a horizontal fibre of radius $a$ with the $x$-axis along the centreline (figure 1). The initial radius of the film measured from the $x$-axis is $r = h_0$.

Figure 1. Schematic of a liquid film on a fibre.

The incompressible NS equations are employed to predict the dynamics of the flow inside the liquid film. To identify the governing dimensionless parameters, we non-dimensionalise the NS equations with the rescaling variables

(2.1a–d)\begin{equation} (x,r,h)= \frac{(\tilde{x},\tilde{r},\tilde{h})}{h_0}, \quad t =\frac{\gamma}{\mu h_0}\, \tilde{t}, \quad \boldsymbol{u} =\frac{\mu}{\gamma}\,\tilde{\boldsymbol{u}}, \quad (p,\boldsymbol{\tau}) = \frac{h_0}{\gamma}\,(\tilde{p},\tilde{\boldsymbol{\tau}}), \end{equation}

where $\tilde {h}$, $\tilde {t}$, $\tilde {\boldsymbol {u}}$ and $\tilde {p}$ represent the dimensional interface height, time, velocity and pressure, respectively (note that the dimensional material parameters are not given tildes). Also, $\tilde {\boldsymbol {\tau }}$ is the shear stress, which is proportional to the strain rate in Newtonian fluids, $\mu$ is the liquid dynamic viscosity, and $\gamma$ is the surface tension of the liquid–gas interface. The dimensionless NS equations can be written as

(2.2)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} =0 , \end{gather}$$

(2.3)$$\begin{gather}\partial_{t} \boldsymbol{u}+ \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = {Oh}^2\,( \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\tau} - \boldsymbol{\nabla} p ) . \end{gather}$$

The non-dimensional quantity ${Oh}=\mu /\sqrt {\rho \gamma h_0}$ is the Ohnesorge number, which relates the viscous forces to inertial and surface-tension forces. Note that in the previous experiment for the slip-enhanced RP instability (Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015), a high-viscosity polymer liquid, entangled polystyrene, is coated on a fibre with radius $10\unicode{x2013}25\ \mathrm {\mu }\mathrm {m}$. The density of the entangled polystyrene is $\rho =1.05\ \mathrm {g}\ \mathrm {cm}^{-3}$, and its liquid–gas surface tension is $\gamma = 30.8\ \mathrm {mN}\ \mathrm {m}^{-1}$. The liquid viscosity depends on temperature, with various values in the range $0.2\unicode{x2013}10\ \mathrm{kg}\ \mathrm {m}^{-1}\ \mathrm {s}^{-1}$ (Haefner Reference Haefner2015). So the Ohnesorge number of liquid film in the experiment is larger than 10, and ${Oh}^2 \gg 1$ in (2.3). It is reasonable to assume that the inertial terms are negligible compared to the viscous terms. With axisymmetric initial perturbations, (2.2) and (2.3) can be reduced to axisymmetric Stokes equations, written as

(2.4)$$\begin{gather} \frac{\partial u}{\partial x} + \frac{1}{r}\,\frac{\partial (vr)}{\partial r} = 0 , \end{gather}$$

(2.5)$$\begin{gather}\frac{\partial p}{\partial x} = \frac{\partial^2 u}{\partial x^2} + \frac{1}{r}\,\frac{\partial}{\partial r} \left( r\,\frac{\partial u}{\partial r} \right) , \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial p}{\partial r} = \frac{\partial^2 v}{\partial x^2} + \frac{\partial}{\partial r} \left[\frac{1}{r}\,\frac{\partial (vr)}{\partial r}\right], \end{gather}$$

where $u$ and $v$ represent the axial and radial velocities, respectively. Since the density of gas around the film is much smaller than that of liquid, the gas flow outside can be assumed to be dynamically passive to simplify the problem. The liquid–gas interface height $h(x,t)$ satisfies the kinematic boundary condition

(2.7)\begin{equation} \frac{\partial h}{\partial t} + u\,\frac{\partial h}{\partial x} = v. \end{equation}

The normal stress balances at the interface $r=h$ give

(2.8)\begin{equation} p- \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{n} , \end{equation}

where $\boldsymbol {n}$ is the outward normal, and $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {n}$ represents the dimensionless Laplace pressure. The tangential force balance is

(2.9)\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{t}=0 , \end{equation}

where $\boldsymbol {t}$ is the tangential vector. With $\boldsymbol {n}$ and $\boldsymbol {t}$ expressed in terms of the unit vectors in the $x$-direction ($\hat {\boldsymbol {e}}_x$) and $r$-direction ($\hat {\boldsymbol {e}}_r$), i.e.

(2.10a,b)\begin{align} \boldsymbol{n} ={-}\frac{\partial_x h}{\sqrt{1+(\partial_x h)^2}}\,\hat{\boldsymbol{e}}_x + \frac{1}{\sqrt{1+(\partial_x h)^2}}\,\hat{\boldsymbol{e}}_r \quad\textrm{and}\quad \boldsymbol{t} = \frac{1}{\sqrt{1+(\partial_x h)^2}}\,\hat{\boldsymbol{e}}_x + \frac{\partial_x h}{\sqrt{1+(\partial_x h)^2}}\,\hat{\boldsymbol{e}}_r, \end{align}

(2.8) and (2.9) give explicitly

(2.11)\begin{align} &p - \frac{2}{1+ \left(\partial_x h \right)^2} \left[ \frac{\partial v}{\partial r} - \frac{\partial h}{\partial x} \left( \frac{\partial u}{\partial r} + \frac{\partial v}{\partial x}\right) + \left(\frac{\partial h}{\partial x} \right)^2 \frac{\partial u}{\partial x} \right] \nonumber\\ &\quad =\frac{1}{h \sqrt{1+\left(\partial_x h \right)^2}} - \frac{\partial_x^2 h}{\left(1+\left(\partial_x h \right)^2\right) ^{3/2}} \end{align}

for the normal forces, and

(2.12)\begin{equation} 2\,\frac{\partial h}{\partial x} \left( \frac{\partial v}{\partial r} - \frac{\partial u}{\partial x} \right) + \left[ 1- \left(\frac{\partial h}{\partial x} \right)^2 \right] \left( \frac{\partial u}{\partial r} + \frac{\partial v}{\partial x} \right) = 0 \end{equation}

for the tangential forces. Here, $\partial _x$ and $\partial _x^2$ refer to first and second spatial derivatives. In terms of the boundary conditions at the fibre surface $r = \alpha$ (where $\alpha$ is the dimensionless fibre radius, i.e. $\alpha = a/h_0$), we have the slip boundary condition in the tangential direction and the no-penetration boundary condition in the normal direction, such that

(2.13a,b)\begin{equation} u = l_{s}\,\frac{\partial u}{\partial r}, \quad v=0. \end{equation}

Here, $l_{s}$ represents the dimensionless slip length rescaled by $h_0$. The governing equations above can be simplified further under the long-wave approximation to give a one-dimensional lubrication equation, compared with experimental results in the previous work (Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015). This lubrication equation is displayed here for the instability analysis appearing in § 3 (see Appendix A for the derivation). Its dimensionless format is

(2.14)\begin{equation} \frac{\partial h}{\partial t} = \frac{1}{ h}\,\frac{\partial}{\partial x} \left[ M(h)\, \frac{\partial p}{\partial x} \right], \end{equation}

where $M(h)$ is the mobility term

(2.15)\begin{equation} M(h) =\frac{1}{16} \left[{-}3h^4-\alpha^4 + 4 \alpha^2 h^2 + 4 h^4 \ln\left(\frac{h}{\alpha}\right) + 4 l_{s}\,\frac{(h^2-\alpha^2)^2}{\alpha} \right], \end{equation}

and $p$ is the Laplace pressure

(2.16)\begin{equation} p = \frac{1}{h} - \frac{\partial^2 h}{\partial x^2} . \end{equation}

3. Instability analysis

In this section, linear instability analysis based on (2.4)–(2.13a,b) is performed using the normal mode method, which has been used widely for the RP instability in different fluid configurations (Rayleigh Reference Rayleigh1878; Tomotika Reference Tomotika1935; Craster & Matar Reference Craster and Matar2006; Liang et al. Reference Liang, Deng, Nave and Johnson2011). The dimensionless perturbed quantities are expressed as

(3.1)\begin{equation} \left.\begin{gathered} u(x,r,t) = \hat{u}(r) \exp({\omega t+{\rm i}kx}),\quad v(x,r,t) = \hat{v}(r) \exp({\omega t+{\rm i}kx})\\ \textrm{and}\quad p(x,r,t) = 1 + \hat{p}(r) \exp({\omega t+{\rm i}kx}), \end{gathered}\right\} \end{equation}

where $\omega$ is the growth rate of perturbations, and $k$ is the wavenumber. Substituting these quantities into the mass and momentum equations (2.4)–(2.6) leads to

(3.2)$$\begin{gather} {\rm i}k \hat{u} + \frac{1}{r}\,\frac{{\rm d} (\hat{v} r)}{{\rm d} r} = 0 , \end{gather}$$

(3.3)$$\begin{gather}{\rm i}k \hat{p} ={-}k^2 \hat{u} + \frac{1}{r}\,\frac{{\rm d}}{{\rm d} r} \left( r\,\frac{ {\rm d} \hat{u}}{{\rm d} r} \right) , \end{gather}$$

(3.4)$$\begin{gather}\frac{{\rm d} \hat{p}}{{\rm d} r} ={-}k^2 \hat{v} + \frac{{\rm d}}{{\rm d} r}\left[\frac{1}{r}\,\frac{{\rm d} (\hat{v} r)}{{\rm d} r}\right] . \end{gather}$$

With (3.2)–(3.4), we can eliminate $\hat {u}$ and $\hat {p}$ to get a fourth-order ordinary partial differential equation for $\hat {v}$:

(3.5)\begin{equation} \frac{{\rm d} }{{\rm d} r}\left\lbrace\frac{1}{r}\,\frac{{\rm d}}{{\rm d} r} \left[r\,\frac{{\rm d}}{{\rm d} r} \left(\frac{1}{r}\,\frac{{\rm d}(\hat{v}r)}{{\rm d}r} \right) \right] \right\rbrace - 2k^2\, \frac{{\rm d}}{{\rm d}r}\left[\frac{1}{r}\,\frac{{\rm d}(\hat{v}r)}{{\rm d}r} \right]+k^4\hat{v} = 0 . \end{equation}

The general solution for (3.5) can be obtained in terms of Bessel functions, written as

(3.6)\begin{equation} \hat{v} = C_1 r\,{\rm K}_0(kr) + C_2\,{\rm K}_1(kr)+ C_3 r\,{\rm I}_0(kr) + C_4\, {\rm I}_1(kr) , \end{equation}

where I and K are modified Bessel function of the first and second kinds, the subscripts $0$ and $1$ represent the order of the Bessel functions, $C_1\unicode{x2013} C_4$ are four arbitrary constants to be determined by the boundary conditions, and $r\in [\alpha,1]$. Substituting (3.6) into (3.2) and (3.3) gives us the solutions for $\hat {u}$ and $\hat {p}$:

(3.7)$$\begin{gather} \hat{u} = \left[ C_1 ( kr\,{\rm K}_1-2\,{\rm K}_0 )+ C_2 k\,{\rm K}_0 - C_3(2\,{\rm I}_0 + kr\,{\rm I}_1 ) - C_4 k\,{\rm I}_0 \right] / ({\rm i}k) , \end{gather}$$

(3.8)$$\begin{gather}\hat{p} = 2 ( C_1\,{\rm K}_0 + C_3\,{\rm I}_0 ) . \end{gather}$$

In a similar approach, the dimensionless perturbed quantities for $\hat {u}, \hat {v}, \hat {p}$, combined with $h(x,t) = 1 + \hat {h}\exp ({\omega t+\textrm {i}kx})$, are substituted into perturbed boundary equations (2.7)–(2.13a,b). For the boundary conditions at the interface ($r=1$), their linearisation gives

(3.9)$$\begin{gather} \frac{{\rm d} \hat{u}}{{\rm d} r} + {\rm i}k \hat{v} = 0, \end{gather}$$

(3.10)$$\begin{gather}\hat{p} - 2\,\frac{{\rm d} \hat{v}}{{\rm d} r} = \hat{h} (k^2-1 ), \end{gather}$$

(3.11)$$\begin{gather}\omega \hat{h} = \hat{v} . \end{gather}$$

And for the boundary conditions on the fibre surface ($r=\alpha$), their linearised forms are

(3.12)$$\begin{gather} \hat{u} = l_{s}\,\frac{{\rm d} \hat{u}}{{\rm d} r}, \end{gather}$$

(3.13)$$\begin{gather}\hat{v} = 0 . \end{gather}$$

According to (3.11), $\hat {h}$ in (3.10) can be eliminated to give the final four equations of the boundary conditions, i.e. (3.9), (3.10), (3.12) and (3.13). Substituting the Bessel functions (3.6)–(3.8) into these perturbed equations leads to a homogeneous system of linear equations for $C_1\unicode{x2013} C_4$, which has a non-trivial solution only if the determinant of the coefficients vanishes. In this way, we have the final equation

(3.14)\begin{equation} \left|\begin{array}{@{}cccc@{}} F_{11} & F_{12} & F_{13} & F_{14} \\ F_{21} & F_{22} & F_{23} & F_{24} \\ \alpha\,{\rm K}_0(k \alpha) & {\rm K}_1(k \alpha) & \alpha\,{\rm I}_0(k \alpha) & {\rm I}_1(k \alpha) \\ k\,{\rm K}_0(k) - {\rm K}_1(k) & k\,{\rm K}_1(k) & k\,{\rm I}_0(k) + {\rm I}_1(k) & k\,{\rm I}_1(k) \\ \end{array}\right| =0 , \end{equation}

where

(3.15)\begin{equation} \left.\begin{gathered} F_{11} = (k^2-1)\,{\rm K}_0(k) - 2 \omega k\,{\rm K}_1(k), \\ F_{12} = (k^2-1)\,{\rm K}_1(k) - 2 \omega \left[ k\,{\rm K}_0(k) + {\rm K}_1(k) \right], \\ F_{13} = (k^2-1)\,{\rm I}_0(k) + 2 \omega k\,{\rm I}_1(k) , \\ F_{14} = (k^2-1)\,{\rm I}_1(k) + 2 \omega \left[ k\,{\rm I}_0(k) - {\rm I}_1(k) \right], \\ F_{21} = (2-l_{s} k^2 \alpha)\,{\rm K}_0(k \alpha) + k(2 l_{s}-\alpha)\,{\rm K}_1(k \alpha) , \\ F_{22} ={-}k\,{\rm K}_0(k \alpha) - l_{s} k^2\,{\rm K}_1(k \alpha), \\ F_{23} = (2-l_{s} k^2 \alpha)\,{\rm I}_0(k \alpha) -k(2 l_{s}-\alpha)\,{\rm I}_1(k \alpha) , \\ F_{24} = k\,{\rm I}_0(k \alpha) - l_{s} k^2\,{\rm I}_1(k \alpha). \end{gathered}\right\} \end{equation}

Because $\omega$ appears only linearly in the first line of the determinant, the dispersion relation between $\omega$ and $k$ can be expressed explicitly as

(3.16)\begin{align} \omega = \frac{k^2-1}{2}\,\frac{{\rm K}_0(k)\,\varDelta_1 - {\rm K}_1(k)\,\varDelta_2 + {\rm I}_0(k)\,\varDelta_3 -{\rm I}_1(k)\,\varDelta_4}{k\,{\rm K}_1(k)\,\varDelta_1- \left[ k\,{\rm K}_0(k)+{\rm K}_1(k) \right] \varDelta_2 -k\,{\rm I}_1(k)\,\varDelta_3 + \left[ k\,{\rm I}_0(k) -{\rm I}_1(k) \right] \varDelta_4} , \end{align}

where

(3.17)\begin{equation} \left.\begin{gathered} \varDelta_1 = \left|\begin{array}{@{}ccc@{}} F_{22} & F_{23} & F_{24} \\ {\rm K}_1(k \alpha) & \alpha\,{\rm I}_0(k \alpha) & {\rm I}_1(k \alpha) \\ k\,{\rm K}_1(k) & k\,{\rm I}_0(k) + {\rm I}_1(k) & k\,{\rm I}_1(k) \end{array}\right|, \\ \varDelta_2 = \left|\begin{array}{@{}ccc@{}} F_{21} & F_{23} & F_{24} \\ \alpha\,{\rm K}_0(k \alpha) & \alpha\,{\rm I}_0(k \alpha) & {\rm I}_1(k \alpha) \\ k\,{\rm K}_0(k) - {\rm K}_1(k) & k\, {\rm I}_0(k) + {\rm I}_1(k) & k\,{\rm I}_1(k) \end{array}\right|, \\ \varDelta_3 = \left|\begin{array}{@{}ccc@{}} F_{21} & F_{22} & F_{24} \\ \alpha\,{\rm K}_0(k \alpha) & {\rm K}_1(k \alpha) & {\rm I}_1(k \alpha) \\ k\,{\rm K}_0(k) - {\rm K}_1(k) & k\,{\rm K}_1(k) & k\,{\rm I}_1(k) \end{array}\right|, \\ \varDelta_4 = \left|\begin{array}{@{}ccc@{}} F_{21} & F_{22} & F_{23} \\ \alpha\,{\rm K}_0(k \alpha) & {\rm K}_1(k \alpha) & \alpha\,{\rm I}_0(k \alpha) \\ k\,{\rm K}_0(k) - {\rm K}_1(k) & k\, {\rm K}_1(k) & k\,{\rm I}_0(k) + {\rm I}_1(k) \end{array}\right|. \end{gathered}\right\} \end{equation}

To compare the Stokes model (3.16) with the lubrication theory, (2.14) is linearised using $h = 1 + \hat {h} \exp ({\omega t + \textrm {i}kx})$, providing the dispersion relation

(3.18)\begin{equation} \omega = (k^2-1) k^2 M, \end{equation}

where $M =[3+ \alpha ^4 - 4 \alpha ^2 + 4 \ln \alpha - 4 l_{s} (1-\alpha ^2)^2/\alpha ]/16.$

In figure 2, predictions of (3.16) and (3.18) in the limiting case ($l_{s}=0$) are found to agree with the results of Craster & Matar (Reference Craster and Matar2006). Note that the lubrication dispersion relation (dashed lines) matches the Stokes dispersion relation (solid lines) when $\alpha \geqslant 0.6$. However, the agreement deteriorates for the smaller $\alpha$, mainly due to the limitations of lubrication approximations for ‘thicker’ films. Despite the failure in predicting the dispersion relation for small $\alpha$, the Plateau instability criterion with critical wavenumber, $k_{crit}=2{\rm \pi} h_0/\lambda _{crit}=1$ (Plateau Reference Plateau1873), can be predicted by the lubrication model for all $\alpha$ values. The reason is that both the circumferential and tangential curvature for the Laplace pressure are kept in (2.16) to give the term $k^2-1$ in the dispersion relation.

Figure 2. The dispersion relation between the growth rate $\omega$ and the wavenumber $k$ on the no-slip boundary condition ($l_{s}=0$) for various dimensionless fibre radii: (a) $\alpha =0.7$ (black), $\alpha =0.6$ (red), $\alpha =0.5$ (blue); (b) $\alpha =0.3$ (black), $\alpha =0.2$ (red), $\alpha =0.1$ (blue). The predictions of the Stokes model (3.16) are plotted in solid lines, and results from the lubrication model (3.18) are illustrated in dashed lines. The triangles and circles are the long-wave and Stokes results of Craster & Matar (Reference Craster and Matar2006), employed here to verify (3.16) and (3.18), respectively.

Figure 3 illustrates the dispersion relations of different $l_{s}$ values. Both the Stokes model and the lubrication model show that the RP instability is enhanced by slip with faster growing perturbations for larger $l_{s}$, consistent with previous findings based on lubrication models (Liao et al. Reference Liao, Li and Wei2013; Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015; Halpern & Wei Reference Halpern and Wei2017). However, obvious discrepancies between the Stokes dispersion relation and the lubrication dispersion relation are found in most cases in figure 3. Even for thin films ($\alpha \geqslant 0.7$), where the lubrication model has been shown to perform well on the no-slip boundary condition (see the agreement of black lines in figures 3a,b), remarkable discrepancies appear when $l_{s} \geqslant 0.5$, indicating that the lubrication model (2.14) is not suitable for large-slip cases. The main reason is that the lubrication approximation imposes a restriction on $l_{s}$, namely $l_{s} \ll \lambda ^2 -1+\alpha$, where $\lambda$ is the perturbation wavelength (see the derivation of (A30) and its corresponding explanation in Appendix A). Moreover, the discrepancy increases for larger $\alpha$, similar to the trend found in the no-slip cases in figure 2.

Figure 3. The dispersion relation between the growth rate $\omega$ and the wavenumber $k$ on different boundary conditions of various dimensionless fibre radii: (a) $\alpha =0.9$, (b) $\alpha =0.7$, (c) $\alpha =0.5$, (d) $\alpha =0.3$. The solid and dashed lines are the predictions of the Stokes model (3.16) and the lubrication model (3.18), respectively, on three boundary conditions with $l_{s}=0$ (black), $l_{s}=0.25$ (blue) and $l_{s}=0.5$ (red).

The $k_{crit}$ values for the RP instability in figure 3 are found to be not affected by the slip. These values are determined by the competition between the two curvature terms for the Laplace pressure, shown on the right-hand side of (2.11), where the term of the circumferential curvature $1/ (h \sqrt {1+(\partial _x h )^2} )$ is the driving force, and the term of tangential curvature $\partial _x^2 h/ (1+(\partial _x h )^2) ^{3/2}$ is the resisting force. The balance of the two forces gives $k_{crit}$, which is independent of $l_{s}$. The dominant wavenumber $k_{max}$ predicted by the Stokes model in figure 3 is shown to be modified by slip, while $k_{max}$ from the lubrication model is unchanged with an analytical expression $k_{max} = \sqrt {1/2}$, derived from (3.18). The variations of $k_{max}$ with slip are illustrated further in figure 4(a), where $k_{max}$ declines with increasing $l_{s}$, possibly leading to larger drop sizes after the film breakup. This trend holds for different fibre radii with more rapid decrease of $k_{max}$ of thicker films (smaller $\alpha$). These findings indicate the limitations of the classical slip lubrication models (Liao et al. Reference Liao, Li and Wei2013; Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015; Halpern & Wei Reference Halpern and Wei2017; Chao et al. Reference Chao, Ding and Liu2018; Ji et al. Reference Ji, Falcon, Sadeghpour, Zeng, Ju and Bertozzi2019), which can be severe with large slip, while the Stokes model proposed here has the capability of predicting the instability over a much larger range of $l_{s}$. Additionally, the lubrication model shows that corresponding growth rates $\omega _{max}$ of each $k_{max}$ grow linearly with increasing $l_{s}$, which can be described by $\omega _{max} = -M/4$, while the Stokes model presents nonlinear growth with smaller $\omega _{max}$.

Figure 4. Variations of the dominant wavenumber $k_{max}$ and corresponding growth rate $\omega _{max}$ with slip lengths. The solid and dashed lines are the predictions of the Stokes model (3.16) and the lubrication model (3.18), respectively, for different fibre thicknesses, i.e. $\alpha =0.9$ (black), $\alpha =0.7$ (blue), $\alpha =0.5$ (red) and $\alpha =0.3$ (purple).

4. Numerical simulations

In this section, direct numerical simulations are performed to support the theoretical findings in § 3 and provide more physical insights into the slip-enhanced RP instability of the cylindrical films on fibres. The open-source framework Basilisk (Popinet Reference Popinet2014, Reference Popinet2018) is employed to solve axisymmetric incompressible NS equations with surface tension. The interface between the high-density liquid and the low-density ambient air is reconstructed by a volume-of-fluid method, which has been validated by different kinds of interfacial flows such as jet breakup (Deblais et al. Reference Deblais, Herrada, Hauner, Velikov, Van Roon, Kellay, Eggers and Bonn2018), bubble bursting (Berny et al. Reference Berny, Deike, Séon and Popinet2020) and breaking waves (Mostert & Deike Reference Mostert and Deike2020).

The simulation domain is a rectangle (the section of a hollow cylinder in cylindrical coordinates) of size $[\alpha, 2]\times [0, L]$. Here, $L$ is the length of the film/fibre, and $\alpha$ is the radius of the fibre. The liquid film of initial radius $h_0 = 1$ is placed at the bottom of the domain, with small perturbations at the liquid–gas interface. The left and right boundary conditions of the simulation domain are considered periodic. The top is a symmetry boundary and bottom is a slip-wall boundary, as described by (2.13a,b). All simulations are performed in the dimensionless unit with the rescaling variables displayed in (2.1a–d). The non-dimensional parameter ${Oh}=10$ comes from liquid properties displayed in the previous experiment (Haefner Reference Haefner2015) with $h_0 = 10\ \mathrm {\mu }\mathrm {m}$, $\rho =1.05\ \mathrm {g}\ \mathrm {cm}^{-3}$, $\gamma = 30.8\ \mathrm {mN}\ \mathrm {m}^{{-1}}$ and $\mu =0.22\ \mathrm {kg}\ \mathrm {m}^{-1}\ \textrm {s}^{{-1}}$. Two configurations with different film lengths and initial perturbations are simulated here to investigate the influence of slip on the growth rates of perturbations (§ 4.1) and the dominant wavelength (§ 4.2), respectively.

4.1. Influence of slip on perturbation growth rates

To study the influence of slip on the growth rates, we consider simulations of a relatively short film $L=10$ with an initial perturbation $h(x,t) =1+ \varepsilon \cos [2 {\rm \pi}(x/L-1/2)]$, where $\varepsilon =0.01$. High-density grids are employed to capture the interface position and velocity profiles inside the film, with $2^{10}$ grid points along the $x$-axis. We choose three slip lengths, i.e. $l_{s}=0$, $0.25$ and $0.5$, for two fibre radii, $\alpha =0.7$ and $0.5$.

Figure 5 shows the simulation results on a no-slip fibre of the radius $\alpha =0.7$. The upper half of figure 5 shows a time evolution of the interface positions of the film, with the initial perturbations starting to grow due to the RP instability. The axial velocity $u(x,r)$ is crucial to understand influence of the slip on the fibre thinning, shown in the lower half of figure 5, where opposite fluxes are found to be directed towards the left and right boundaries, respectively, as the perturbation increases. More velocity fields are presented in figures 6(a,c,e) for three different slip lengths at $t=100$. According to the variations of the contours in figures 6(a,c,e), the wall slip is found to accelerate the instability with larger $u(x,r)$. The velocity vectors in the lower plots of figures 6(a,c,e) show that slip decreases the velocity gradient $\partial u/ \partial r$ near the fibre. The velocity profiles at different $x$-coordinates (see the dashed lines with arrows in figures 6a,c,e) are also extracted from the contours, shown in figures 6(b,d, f), where the velocity on the fibre $u(x,\alpha )$ and its gradient $\partial _r u(x,\alpha )$ are calculated numerically. According to (2.13a,b), $l_{s} = u(x,\alpha ) / \partial _r u(x,\alpha )$. So we can have numerically predicted slip lengths at different positions, shown by the dashed lines in figures 6(d, f). These values agree with the input $l_{s}$ in the boundary condition of the simulations, indicating that the numerical solutions can capture the flow features on the slip boundary conditions.

Figure 5. Film thinning simulation results of one perturbation wave in an axially symmetric domain. For this case, $L=10$, $l_{s}=0$ and $\alpha =0.7$. Interface profiles $h$ (upper half) and contours of the axial velocity $u(x,r,t)$ inside the film (lower half) are shown at three time instants: (a) $t_1=300$, (b) $t_2=600$, (c) $t_3=900$. Here, $h_{min}$ represents the minimum radius of the film.

Figure 6. Axial velocity fields at $t=100$ for different slip lengths: (a,b) $l_{s} = 0$, (c,d) $l_{s} = 0.25$, (e, f) $l_{s} = 0.5$. Here, the fibre radius is $\alpha = 0.7$. (a,c,e) Contours of the entire configuration (upper plots) and velocity vectors of the local field near $x=6$ (lower plots). (b,d, f) Velocity profiles at four local positions $x$: $5.5$ (black), 6.0 (purple), 6.5 (blue) and 7.0 (red), shown by dashed lines with arrows in (a,c,e).

Figure 7 illustrates the growth rates of perturbations in the six cases simulated here. Because $h(x,t) = 1+\hat {h}\,\textrm {e}^{\textrm {i}kx}\,\textrm {e}^{\omega t}$ is employed in the instability analysis (§ 3), we set the $y$-coordinate as $\ln (1-h_{min})$ to present the linear growth of perturbations. Here, the numerical solutions agree with predictions of the Stokes model, showing the capability of the Stokes model to describe the RP instability of films on slippery fibres, while the predictions of the lubrication model are not accurate for the slip cases. These results are consistent with dispersion relations in figure 3 and confirm the Stokes model numerically.

Figure 7. Linear time evolution of the minimum radii of films $h_{min}(t)$ on two fibre radii: (a) $\alpha =0.7$, (b) $\alpha =0.5$. The numerical solutions (dotted lines with symbols) are compared to the predictions of both the lubrication model (dashed lines) and Stokes model (solid lines) for three slip lengths: $l_{s} = 0$ (black), $l_{s}=0.25$ (red) and $l_{s}=0.5$ (blue). The inset illustrates $h_{min}(t)$ in the uniform grids. Here, $t$ is scaled by $\mu h_0/\gamma$.

Based on the approach proposed by Liao et al. (Reference Liao, Li and Wei2013) and Wei, Tsao & Chu (Reference Wei, Tsao and Chu2019), we have a physical explanation for why large slip causes the failure of lubrication models, i.e. the radial derivative of the axial velocity is no longer much larger than its axial derivative, violating the lubrication approximation $| {\partial ^2 u}/{\partial x^2} | / |({1}/{r}) ({\partial }/{\partial r})(r ({\partial u}/{\partial r}) )| \ll 1$ (see more details in Appendix A). Figure 8 shows the contour inside the film of the ratio $| {\partial ^2 u}/{\partial x^2} | / |({1}/{r})({\partial }/{\partial r})(r ({\partial u}/{\partial r}) )|$, extracted from the numerical solutions in this section. The ratio is found to increase from about $0.01$ for the no-slip case to $0.1$ for the slip case ($l_{s}=0.5$) at the same location, which is consistent with the physical explanation above.

Figure 8. Ratio of the velocity gradients at $t=100$ for different slip lengths: (a) $l_{s} = 0$, (b) $l_{s} = 0.25$, (c) $l_{s} = 0.5$. Here, the fibre radius is $\alpha = 0.7$, and $u_{xx}$ and $(r u_r)_r/r$ represents the second-order axial and radial derivatives of the axial velocity in (A24), respectively.

4.2. Influence of slip on the dominant wavelength

For the influence of slip on the dominant wavelengths, we consider a configuration of a long film $L=160$ with $2^{13}$ grid points along the $x$-axis. Eight different cases are simulated with four slip lengths, i.e. $l_{s}=0$, $1$, $2$ and $3$ on fibres of two radii, $\alpha =0.7$ and $0.5$.

The simulations start from random initial perturbations $h(x,0) = 1+ \varepsilon \,N(x)$, where $\varepsilon = 2 \times 10^{-3}$, and the random number $N(x)$ follows a normal distribution with zero mean and unit variance. The random initial perturbations are designed to mimic numerically the arbitrary disturbances in the experiment (Haefner et al. Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015). Figure 9 presents the spatial and temporal evolutions of interface profiles $h(x,t)$, where these small random perturbations, driven by the surface tension in the RP instability, grow gradually against time to generate significant capillary waves at early stages, plotted in figure 9(b). Figure 9(c) shows that wavelengths of these capillary waves are unchanged at the later stage of perturbation growth.

Figure 9. (a) Space–time plot of interface positions $h(x,t)$ on a no-slip fibre of the radius $\alpha =0.7$. Darker (brighter) colours correspond to the smaller (larger) values of $h(x,t)$. (b,c) Interface profiles extracted from (a) at early and late stages of the instability.

To further prove the previous findings and confirm theoretical predictions of the Stokes model quantitatively, multiple independent simulations (five for each case) are performed to gather statistics of the dominant modes with a statistical methodology proposed by Zhao, Sprittles & Lockerby (Reference Zhao, Sprittles and Lockerby2019). For each realisation, a discrete Fourier transform of the interface position $h(x)$ is applied to get the power spectral density (PSD) of the perturbations. The square root of the ensemble-averaged PSD ($H_{rms}$) at each time is plotted in figures 10(a,b) with a modal distribution (spectrum) fitted by the Gaussian function, where the peak represents the dominant wavenumber $k_{max}$ (see black dash-dotted lines). Extracting $k_{max}$ from the fitted spectrum at each time instant yields the insets in figures 10(a,b). Promisingly, $k_{max}$ converges rapidly to a constant in both cases, further supporting the findings in figure 9(c). Here, we compare these constants with the $k_{max}$ predicted by the Stokes model. Similar statistical analysis is applied for other cases to generate the symbols ($\lambda _{max} = 2 {\rm \pi}/k_{max}$) in figure 10(c), giving good agreement with theoretical predictions. Therefore, we can conclude that the dominant modes depend strongly on slip lengths with longer perturbation waves formed on a more slippery surface.

Figure 10. (a,b) The root mean square (r.m.s.) of non-dimensional perturbation amplitude versus non-dimensional wavenumber on fibres with different slip lengths, at four time instants: (a) $2\times 10^3$ (green), $4.5\times 10^3$ (black), $6\times 10^3$ (red) and $7\times 10^3$ (blue); (b) $1.2\times 10^2$ (green), $3.3\times 10^2$ (black), $4.2\times 10^2$ (red) and $4.8\times 10^2$ (blue). The inset shows the time history of the dominant wavenumber. (c) Variations of the dominant wavelengths with slip lengths: a comparison between the theoretical predictions of (3.16) and numerical solutions for two fibres of different radii.

Figure 11 shows the influence of slip on the droplet size. In figure 11(a), larger drops are formed on slippery fibres with obviously longer wavelengths compared to ones achieved on the no-slip boundary, consistent with the theoretical results in figure 10(c) qualitatively. Furthermore, we define the area between two bottom points of $h(x)$ as one droplet. Performing statistical analysis for all realisations gives the droplet-size distribution of slip and no-slip cases, illustrated in figure 11(b), whose mean and standard deviation are achieved by fitting the Gaussian function (red dash-dotted lines). The slip case has a larger mean value of the distribution, further proving the finding in figure 11(a). Besides the mean, slip is also found to lead to a wider distribution, with a larger standard deviation, possibly because the slip wall imposes fewer restrictions (shear stresses on the fibre, i.e. $\tilde {\tau }_{w} = \mu \,\partial _{\tilde {r}} \tilde {u}(\tilde {x},a)$) on the liquid compared to the no-slip wall. The perturbation wavelengths have more ‘freedom’ to expand or shrink, leading to a wider wavelength distribution. Note that figure 11 might not be the final equilibrium state. So droplet coalescence could happen with complicated nonlinear dynamics. However, these nonlinear behaviours are beyond the scope of this work and should be the subject of future investigation.

Figure 11. (a) Interface profiles of liquid films on fibres of radius $\alpha =0.7$. The results for two slip lengths are presented by a solid blue line for $l_{s}=0$, and a red dash-dotted line for $l_{s}=3$. (b) Distributions of the droplet size on the no-slip (upper plot) and slip (lower plot).

5. Comparisons with experimental results

Besides the numerical validations, the Stokes model is compared further with the experimental results in the previous work done by Haefner et al. (Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015). The experiments were done by measuring the evolutions of entangled polystyrene films with homogeneous initial thicknesses on no-slip and slip fibres of radius $a=9.6\ \mathrm {\mu }\textrm {m}$.

Figure 12 displays the dimensional dominant wavelengths $\tilde {\lambda }_{max}$ ($\lambda _{max} h_0$) of liquid films with different initial thicknesses $h_0$, where the blue triangles and red dots represent experimental data of the films on no-slip and slip fibres, respectively. Since $\lambda _{max}$ predicted by the lubrication model (3.18) is the constant ($\lambda _{max}=2\sqrt {2}{\rm \pi}$), independent of both $l_{s}$ and $\alpha$, $\lambda _{max} h_0$ grows linearly with increasing $h_0$ (see the black dashed line in figure 12). However, the lubrication model seemingly underestimates experimentally obtained $\lambda _{max}$. The possible reason is that most $\alpha$ values investigated in the experiment are smaller than $0.5$, outside the range of validity of the lubrication approximation. A series of dominant wavelengths with different $\alpha$ values are extracted from the Stokes model (3.16) to generate a blue solid line for no-slip cases, and a red dash-dotted line for slip cases with a fitting $\tilde {l}_{s}$ in figure 12. Note that the dimensionless slip length varies along the red dash-dotted curve with different $h_0$ values. Though both the blue and red curves still look linear due to small changes of $\lambda _{max}$, obvious differences between the Stokes predictions and the lubrication ones are found for the thick films ($h_0>40\ \mathrm {\mu }\textrm {m}$). Promisingly, despite dispersion of experimental data, the Stokes model can provide better predictions for dominant wavelengths of both no-slip cases and slip cases compared to the lubrication model.

Figure 12. Influence of the geometry on the dominant wavelength of the instability: a comparison between the theoretical predictions and experimental data. The experimental data come from the work of Haefner et al. (Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015).

In addition to the dominant wavelengths, we then make a further comparison between the experimental growth rates and theoretical ones, illustrated in figure 13, where the growth rate is non-dimensionalised by the ratio of the capillary velocity to the fibre radius, i.e. $\omega =({a}\tilde {\omega }/({\gamma /\mu }))$. However, using the parameters from Haefner et al. (Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015) ($\gamma /\mu = 294\ \mathrm {\mu }\mathrm {m}\ \min ^{-1}$ and $\tilde {l}_{s} = 0.3\,a$), the Stokes model is found to provide worse agreements with the experimental data than the lubrication model. Noticeably, predictions of the lubrication model match the experimental data even in the case of thick films ($h_0/a > 3$), which is unreasonable due to the thin-film assumption in the lubrication theory. The reason for the disagreement is currently unclear and should be the subject of future investigation.

Figure 13. Influence of the geometry on the growth rate of the dominant mode on no-slip and slip fibres. The solid lines and dashed lines are the predictions of the Stokes model and the lubrication model, respectively. Blue lines represent no-slip results, and red lines are the slip ones. The experimental data come from the work of Haefner et al. (Reference Haefner, Benzaquen, Bäumchen, Salez, Peters, McGraw, Jacobs, Raphaël and Dalnoki-Veress2015).