2.1. Problem definition and metric terms in general coordinates

In this section, we introduce the essential differential geometry tools to solve the problem of the coating on a generic substrate. For a complete description of differential geometry and general coordinates, we refer to Deserno (Reference Deserno2004). The derivation of the lubrication equation for generic curved substrates can be found in Roy et al. (Reference Roy, Roberts and Simpson2002), Thiffeault & Kamhawi (Reference Thiffeault and Kamhawi2006) and Wray, Papageorgiou & Matar (Reference Wray, Papageorgiou and Matar2017). The geometry is sketched in figure 2. We consider a generic substrate ${h}^{0}$, on which lies a fluid film of thickness ${h}$, and introduce a Cartesian reference frame $({x},{y},{z})$. The substrate is identified by the position vector $\boldsymbol {X}(x^{\{1\}}, x^{\{2\}})$, where ($x^{\{1\}},x^{\{2\}}$) denote the local coordinates used to parameterize the surface (e.g. the zenith and the azimuth for spherical coordinates, the radial coordinate and the azimuth for a cone). The flow equations are solved in the local and natural reference frame of the substrate. We introduce the local coordinate vectors parallel to the substrate ${\boldsymbol {e}}_{i}=\partial _i \boldsymbol {X}$, $i=1,2$ (not necessarily orthonormal), and the normal coordinate vector ${\boldsymbol {e}}_{3}={\boldsymbol {e}_{1}\times \boldsymbol {e}_{2}}/{|\boldsymbol {e}_{1}\times \boldsymbol {e}_{2}|}$. From the knowledge of the local coordinate vectors, we introduce the $2\times 2$ symmetric metric tensor components $\mathsf{G}_{i j}$ and the square root of the determinant of the metric on the substrate $w$, which is related to the area element on the surface $\mathrm {d}A$ through ${\rm d}A=w\,\mathrm {d}\kern0.06em x^{\{1\}}\,\mathrm {d}\kern0.06em x^{\{2\}}$

(2.1a,b)\begin{equation} \mathsf{G}_{i j}=\boldsymbol{e}_{i} \boldsymbol{\cdot} \boldsymbol{e}_{j}=\mathsf{G}_{j i}, \quad w=(\operatorname{det} \mathsf{G}_{i j})^{1 / 2}. \end{equation}

The metric tensor defines the generic line element $\mathrm {d} s$ as $\mathrm {d} s^{2}=\mathsf{G}_{11} (\mathrm {d}\kern0.06em x^{\{1\}})^{2}+2\mathsf{G}_{12} \,\mathrm {d}\kern0.06em x^{\{1\}} \,\mathrm {d}\kern0.06em x^{\{2\}}+\mathsf{G}_{22} (\mathrm {d}\kern0.06em x^{\{2\}})^{2}$. Therefore, the dimensions of each component depend on the considered parameterization $(x^{\{1\}},x^{\{2\}})$, so that each part that composes $\mathrm {d} s^{2}$ has the dimensions of the square of a length. We also introduce the second fundamental form and the curvature tensor, which respectively read, following Einstein's notation for the summation

(2.2a,b)\begin{equation} {\mathsf{S}}_{i j}=\partial_i \boldsymbol{e}_{j} \boldsymbol{\cdot} {\boldsymbol{e}}_3, \quad \mathsf{K}_{i}^{\{j\}}={\mathsf{S}}_{i k}\mathsf{G}^{\{kj\}}, \end{equation}

where $\mathsf{G}^{\{ij\}}$ are the inverse metric tensor components, i.e. $\mathsf{G}^{\{ij\}}=\mathsf{G}_{ij}^{-1}$. The mean $\mathcal {K}$ and the Gaussian $\mathcal {G}$ curvatures read $\mathcal {K}=\operatorname {tr}\boldsymbol{\mathsf{K}}$ and $\mathcal {G}=\operatorname {det} \boldsymbol{\mathsf{K}}$, respectively. A generic vector $\boldsymbol {f}$ can be written in terms of its covariant and contravariant base, i.e. $\boldsymbol {f}=f^{\{i\}}\boldsymbol {e}_i=f_i\boldsymbol {e}^{\{i\}}$, where $\boldsymbol {e}^{\{i\}}$ is the covector defined as ${\boldsymbol {e}}^{\{i\}} \boldsymbol {\cdot } {\boldsymbol {e}}_{j}=\delta _{ij}$. The two contravariant components, parallel to the substrate, of the gravity vector read $g_t^{\{i\}}=\boldsymbol {g}\boldsymbol {\cdot } \boldsymbol {e}^{\{i\}}$, while the normal one reads $g_3=\boldsymbol {g}\boldsymbol {\cdot }{\boldsymbol {e}}_3$. The gradient of a scalar function $f$ and the divergence of a generic vector $\boldsymbol {f}=f^{\{i\}}{\boldsymbol {e}}_i$ respectively read (Irgens Reference Irgens2019)

(2.3a,b)\begin{equation} \boldsymbol{\nabla} f={\partial_i f} \,\mathsf{G}^{\{i j\}} \boldsymbol{e}_{j}=\partial^{\{i\}} f\, \boldsymbol{e}_{j}, \quad \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{f} =w^{{-}1} \partial_{i}(w f^{\{i\}}). \end{equation}

The above-defined quantities and differential operators are enough to describe the coating problem on a generic substrate, introduced in the following section.

Figure 2.

Sketch of the coordinate systems employed in this analysis. A global Cartesian reference frame $(x,y,z)$ is considered. At each point, the position of the substrate is identified by the vector $\boldsymbol {X}$, which depends on the chosen parameterization $(x^{\{1\}},x^{\{2\}})$ of the substrate. The derivatives of the position vector identify the local reference frame on the substrate, on which the lubrication equation is solved.

2.2. Lubrication equation and drainage solution

We consider a thin viscous film, flowing on a substrate ${h}^{0}$, of thickness ${h}$ measured along the direction perpendicular to the substrate itself. The constant fluid properties are the density $\rho$, viscosity $\mu$ and the surface tension coefficient $\gamma$. In the absence of inertia, the lubrication model for a generic curved substrate was first derived by Roy et al. (Reference Roy, Roberts and Simpson2002) via central manifold theory. We non-dimensionalize the thickness with $h_i$ and the tangential directions with $R$, i.e. a characteristic film thickness (e.g. the initial one, if uniform) and a relevant length of the substrate (e.g. its equatorial radius), respectively. We introduce the drainage time scale $\tau ={\mu R}/({\rho g h_i^{2}}).$ Upon non-dimensionalization, the equation in coordinate-free form reads (Roy et al. Reference Roy, Roberts and Simpson2002; Howell Reference Howell2003; Roberts & Li Reference Roberts and Li2006; Thiffeault & Kamhawi Reference Thiffeault and Kamhawi2006)

(2.4)\begin{align} &(1-\delta\mathcal{K}h+\delta^{2}\mathcal{G}h^{2})\frac{\partial h}{\partial t}+\frac{1}{3Bo} \boldsymbol{\nabla}\boldsymbol{\cdot}\left[h^{3} \left( \boldsymbol{\nabla} \tilde{\kappa}-\frac{\delta}{2} h(2\mathcal{K} \boldsymbol{\mathsf{I}}-\boldsymbol{\mathsf{K}}) \boldsymbol{\cdot}\boldsymbol{\nabla}\mathcal{K}\right)\right]\nonumber\\ &\quad +\frac{1}{3} \boldsymbol{\nabla}\boldsymbol{\cdot}\left[h^{3} \left({ \boldsymbol{g}}_{t}- \delta h\left(\mathcal{K} \boldsymbol{\mathsf{I}}+\frac{1}{2} \boldsymbol{\mathsf{K}}\right) \boldsymbol{\cdot} {\boldsymbol{g}}_{t}+ \delta {{g}}_{3} \boldsymbol{\nabla} h\right)\right]=0, \end{align}

where $\tilde {\kappa }=\mathcal {K}+\delta (\mathcal {K}^{2}-2\mathcal {G})h+\delta \nabla ^{2}h$ is the free-surface curvature, $Bo=(\rho g R^{2})/\gamma$ is the Bond number, $\delta =h_i/R$ is the aspect ratio of the thin film and ${\boldsymbol {g}}_{t}$ and ${{g}}_{3}$ identify the gravity vector components tangent and normal to the substrate, respectively. The terms in the first and second brackets represent the flux induced by capillary and gravity effects, respectively. Capillary flow is induced, at leading order, by variations of the mean curvature of the substrate $\mathcal {K}$ and leads to film thinning and thickening in the neighbourhood of local maximum and minimum values of the curvature (Roy et al. Reference Roy, Roberts and Simpson2002). Corrections at order ${O}(\delta )$ introduce (i) free-surface curvature variations and (ii) higher-order terms of the substrate curvature. Gravity-induced fluxes are instead related, at leading order, by the gravity components tangential to the substrate $\boldsymbol {g}_t$. Corrections at ${O}(\delta )$ introduce hydrostatic pressure gradients along the thin film. This equation can be written in compact form as follows:

(2.5)\begin{equation} (1-{\delta\mathcal{K}}{h}+{\delta^{2}\mathcal{G}}{h}^{2})\frac{\partial {h}}{\partial {t}}+\frac{1}{3} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{q}} =0,\end{equation}

so-called conservation form, where ${\boldsymbol {q}} ={q}^{\{1\}}\boldsymbol {e}_1+{q}^{\{2\}}\boldsymbol {e}_2$ is the flux. The so-called drainage problem relies on two assumptions. The Bond number is assumed to be very large, i.e. $R^{2}/\ell _c^{2} \gg 1$, where $\ell _c=\sqrt {\gamma /(\rho g)}$ is the capillary length. After this first assumption, the problem accounts only for gravity effects resulting from drainage and hydrostatic pressure gradients. The latter terms are important when the aspect ratio $\delta =h_i/R$ is not negligible, i.e. for a thick film on a large substrate (compared with the capillary length) with small radius of curvature (compared with the film thickness). A limit case occurs when the substrate is locally flat. The leading-order solution is given by the hydrostatic pressure terms, since drainage is absent. A case in which both capillary and hydrostatic effects cannot be neglected occurs instead when the radius of curvature of the substrate is comparable to the film thickness, i.e. regions of extremely large curvature such as the tip of a cone. In the following, we restrict ourselves to the situation in which the film is very thin and the substrate does not present regions of extremely large curvature. Therefore, also $\delta \ll 1$ is considered and the drainage problem reads

(2.6)\begin{equation} \frac{\partial h}{\partial t}+\frac{1}{3} \boldsymbol{\nabla}\boldsymbol{\cdot}[h^{3} { \boldsymbol{g}}_{t} ]=0. \end{equation}

In this case, the flux per unit length is defined as $\boldsymbol {q}=q^{\{1\}}\boldsymbol {e}_1+q^{\{2\}}\boldsymbol {e}_2=h^{3}g_t^{\{1\}}\boldsymbol {e}_1+h^{3}g_t^{\{2\}}\boldsymbol {e}_2$. The solution of the drainage problem requires only the knowledge of $w$ and the tangential gravity vector components ${g}_{t}^{\{i\}}$.

The numerical implementation of the lubrication equation (2.4) is performed in the finite-element solver COMSOL Multiphysics, in which the lubrication equation is implemented in its conservation form (2.5). Quadratic Lagrangian elements are exploited for the numerical discretization, while the time marching is performed with the built-in backward differentiation formula (BDF) solver. In the case of (2.4), we solve for the variables $(h,\tilde {\kappa })$. We refer to the corresponding sections for more detail about the boundary conditions for the different substrates.

The validation procedure consists of a first mesh size validation. We thus verify the faithfulness of the employed parameterization $\boldsymbol {X}(x^{\{1\}},x^{\{2\}})$ by a comparison with the parameterization $\boldsymbol {X}=(x,y,h^{0}(x,y))$, so-called Monge parameterization (Thiffeault & Kamhawi Reference Thiffeault and Kamhawi2006; Mayo et al. Reference Mayo, McCue, Moroney, Forster, Kempthorne, Belward and Turner2015), reported in the electronic supplementary material (ESM available at https://doi.org/10.1017/jfm.2022.758) together with the other parameterizations employed in this work. We also verify the non-dimensionalization by solving the dimensional version of (2.4) and comparing the solution for each substrate with the non-dimensional model. To illustrate and complement the theoretical results, we finally compare in § 6 the drainage problem results with experiments performed following the procedure outlined in Lee et al. (Reference Lee, Brun, Marthelot, Balestra, Gallaire and Reis2016) and Jones et al. (Reference Jones, Jambon-Puillet, Marthelot and Brun2021), for diverse substrates.

2.3. Asymptotic theory – general expression for the thickness at a local maximum

The employed substrate-free expression of the lubrication equation is suitable for analytical results. In this section, we develop a general expression for the thickness at a local maximum of the substrate. In the vicinity of the local maximum, the smooth substrate is described through a Monge parameterization of the substrate, i.e. $(x^{\{1\}}=x,x^{\{2\}}=y)$, see the ESM for further detail. The generic substrate position vector in the vicinity of the maximum reads $\boldsymbol {X}=(x,y,h^{0}(x,y))$. The square root of the determinant of the metric tensor $w$ and the tangential gravity components respectively read

(2.7a–c)\begin{equation} w=\sqrt{1+\left(\partial_x h^{0}\right)^{2}+\left(\partial_y h^{0}\right)^{2}}, \quad g_t^{\{1\}}={-}\frac{\partial_x h^{0}}{w^{2}}, \quad g_t^{\{2\}}={-}\frac{\partial_y h^{0}}{w^{2}}. \end{equation}

We expand the drainage solution in the vicinity of the maximum identified by the point $\boldsymbol {x}=(x,y)=\textbf {0}$ by employing an asymptotic expansion in the spatial variables, i.e.

(2.8)\begin{equation} h(x,y,t)=H_0(t)+H_{11}(t)x +H_{12}(t)y + \dots .\end{equation}

Upon substitution of the decomposition (2.8) in (2.6), the ${O}(1)$ problem for $H_0(t)$ reads

(2.9) \begin{align} \frac{H_0'}{H_0^{3}}&=\left.\frac{ (\partial_{xx}h^{0}(1+ \left(\partial_{y}h^{0}\right)^{2})-2 \partial_{x}h^{0} \partial_{xy}h^{0} \partial_{y}h^{0}+\partial_{yy}h^{0} (\left(\partial_{x}h^{0}\right)^{2}+1))}{3 ((\partial_{y} h^{0})^{2}+(\partial_{x}h^{0})^{2}+1)^{2}}\right\vert_{\boldsymbol{x}=\textbf{0}}\nonumber\\ &={-}\left.\left(\frac{1}{3}\left(\boldsymbol{g}\boldsymbol{\cdot}{\boldsymbol{e}}_3\right) \mathcal{K}\right)\right\vert_{\boldsymbol{x}=\textbf{0}}, \end{align}

with the initial condition $H_0(0)=1$ in the case of a unitary initial thickness. At the maximum location, $\partial _{x}h^{0}=\partial _{y}h^{0}=0$, i.e. the normal vector and gravity are aligned. Therefore, the quantity $(\frac {1}{3}(\boldsymbol {g}\boldsymbol {\cdot }{\boldsymbol {e}}_3)\mathcal {K})\vert _{\boldsymbol {x}=\textbf {0}}$ simplifies to $\partial _{xx}h^{0}+\partial _{yy}h^{0}=-\mathcal {K}_p$, where $\mathcal {K}_p$ is the opposite of the mean curvature at the maximum. The resulting problem and associated solution read

(2.10)\begin{equation} \frac{H_0'(t)}{H_0(t)^{3}}={-}\frac{1}{3}\mathcal{K}_p \rightarrow H_0(t)=\frac{1}{\sqrt{\dfrac{2 \mathcal{K}_p t}{3}+1}}. \end{equation}

A general expression for the drainage in the vicinity of a local maximum is obtained. Note that the mean curvature at the local maximum is negative, and thus $\mathcal {K}_p>0$. The thickness at a local maximum depends on the mean curvature. From a geometrical point of view, the mean curvature represents variations of the tangential vectors along the surface. Since the normal to the surface and the gravity vector are aligned, the mean curvature determines the evolution of the tangential components of the gravity field, in the vicinity of the local maximum. In particular, an increase of $\mathcal {K}_p$ implies larger values of gravity in the tangent plane of the substrate moving away from the pole and thus a faster drainage and a lower thickness, and vice versa. This solution allows one to identify the limits of the considered drainage model. If $\mathcal {K}_p=0$, i.e. a locally flat substrate, there is no drainage and, therefore, the thickness is constant and equal to $H_0=1$. In this case, hydrostatic effects cannot be neglected since they are the leading-order effect and lead to a time-dependent drainage (Huppert & Simpson Reference Huppert and Simpson1980). Therefore, the drainage model is not suitable to describe locally flat substrates. A second limiting case occurs when $\mathcal {K}_p \rightarrow \infty$, i.e. the radius of curvature in the vicinity of the local maximum tends to zero. A classical example is the tip of a cone, parameterized with the radius $x^{\{1\}}=r$ and the azimuth $x^{\{2\}}=\varphi$. In this case, an exact solution $h \propto \sqrt {r}$ is obtained (see the ESM), which presents a zero thickness at the pole, in accordance with solution (2.10), which tends to zero as $\mathcal {K}_p \rightarrow \infty$. In that case, hydrostatic and capillary effects are crucial to define the thickness distribution in the vicinity of the tip. Another important limitation comes from the considered geometry. When the fluid is located below the substrate, the solution is formally analogous to (2.10), and predicts a progressive thinning. However, it is well known, in these situations, that hydrostatic pressure gradients and capillary effects play a key role since the Rayleigh–Taylor instability can occur (Balestra, Nguyen & Gallaire Reference Balestra, Nguyen and Gallaire2018b). In the case of converging geometries with a local minimum, the solution is formally analogous, but with $\mathcal {K}_p <0$. In this case, the thickness progressively increases and tends to infinity for $t =-3/(2\mathcal {K}_p) >0$, long after hydrostatic and capillary effects should have been considered. These combined effects contribute indeed to the Rayleigh–Taylor instability when the fluid lies below the substrate and for a levelling and flattening of the interface when it lies above.

Turning back to the situation where $\mathcal {K}_p > 0$ and remains finite, one can take the limit for $t \gg 1$ of solution (2.10), leading to

(2.11)\begin{equation} H_0(t)=\frac{1}{\sqrt{\dfrac{2 \mathcal{K}_p t}{3}}}+{O}\left(\frac{1}{t^{3/2}}\right),\end{equation}

i.e. the solution is independent of the initial condition. This result is in agreement with the analysis in Lee et al. (Reference Lee, Brun, Marthelot, Balestra, Gallaire and Reis2016), where the authors theoretically and experimentally showed an insensitivity of the film thickness with respect to the initial conditions.

In the following, we investigate the spatial evolution of the thin film thickness when moving away from the pole. We initially consider the case of an axisymmetric substrate, the spheroid.

3.1. Drainage problem

In this section, we consider the drainage of a thin film flowing on an spheroidal substrate of equatorial radius $R$ (i.e. $a=b=1$) and height $cR$. We non-dimensionalize the in-plane directions and substrate variables with the equatorial radius $R$. We parameterize the spheroidal surface via the zenith (or colatitude) ${x^{\{1\}}}=\vartheta$ and the azimuth ${x^{\{2\}}}=\varphi$

(3.1)\begin{equation} \boldsymbol{X}(\vartheta,\varphi)=(\sin\vartheta\cos \varphi, \sin \vartheta \sin \varphi,c \cos \vartheta). \end{equation}

A complete description of the metric and curvature tensors is reported in the ESM. The gravity term $g_t^{\{1\}}$ and $w$ are

(3.2a,b)\begin{equation} g_t^{\{1\}}(\vartheta)=\frac{c \sin (\vartheta)}{c^{2} \sin ^{2}(\vartheta)+\cos ^{2}(\vartheta)}, \quad w(\vartheta)=\frac{1}{\sqrt{2}}\sin (\vartheta){\sqrt{(1-c^{2})\cos (2 \vartheta)+1+c^{2}}}, \end{equation}

while $g_t^{\{2\}}=0$. In the case $c=1$, we recover the evolution equation for the spherical case, reported in the ESM. Following the previous section, we consider as initial condition a constant thickness on the substrate, i.e. $h(\vartheta,0)=1$. The problem is solved through an asymptotic expansion in the vicinity of the pole. We expand the solution at different orders in $\vartheta$

(3.3)\begin{equation} h(\vartheta,t)=H_0(t)+\vartheta^{2}H_2(t)+\vartheta^{4}H_4(t)+\vartheta^{6}H_6(t)+\dots . \end{equation}

We introduce this ansatz and expand in powers of $\vartheta$. At each order ${O}(\vartheta ^{n})$, one obtains an ordinary differential equation for $H_n$. The problem at orders ${O}(1)$ and ${O}(\vartheta ^{2})$ together with their solution read

(3.4a,b)\begin{gather} \frac{2}{3} cH_0^{3}+H_0'=0, \quad H_0(0)=1 \rightarrow H_0=\left(\frac{4 c t}{3}+1\right)^{{-}1/2},\end{gather}

(3.5)\begin{gather} \left.\begin{gathered} \left(\frac{2 c}{3}-c^{3}\right)H_0^{3}+4 cH_0^{2}H_2+H_2'=0, \quad H_2(0)=0,\\ \rightarrow H_2=\frac{(3 c^{2}-2) (64 \sqrt{3} c^{3} t^{3}+144 \sqrt{3} c^{2} t^{2}+108 \sqrt{3} c t+27 (\sqrt{3}-\sqrt{4 c t+3}))}{10 (4 c t+3)^{7/2}}. \end{gathered}\right\} \end{gather}

Applying the same procedure at orders ${O}(\vartheta ^{4})$ and ${O}(\vartheta ^{6})$, the solution up to ${O}(\vartheta ^{6})$ and at leading order for $t \gg 1$, reads (see Appendix A for further detail)

(3.6)\begin{align} h(\vartheta,t)&=\sqrt{\frac{3}{4{t}}}\frac{1}{\sqrt{c}} \left(1+\frac{1}{10} \left(3 c^{2}-2\right) \vartheta^{2}-\frac{\left(336 c^{4}-408 c^{2}+31\right) \vartheta^{4}}{4800}\right. \nonumber\\ &\quad \left.+\frac{\left(58\,464 c^{6}-115\,368 c^{4}+62\,667 c^{2}-4576\right)\vartheta^{6}}{1\,584\,000}\right)+{O}(\vartheta^{8}) +{O}\left(\frac{1}{t^{3/2}}\right) . \end{align}

Note that the ${O}(1)$ large-time solution is formally analogous to (2.10) with $\mathcal {K}_p=2 c$, i.e. the opposite of the mean curvature at the pole.

We perform numerical simulations of (2.6), in the region $0 <\vartheta <{\rm \pi} /2$. Owing to the hyperbolic nature of the equation, no boundary conditions are necessary at $\vartheta =0$ and $\vartheta ={\rm \pi} /2$, and thus we impose only the initial condition $h(\vartheta,0)=1$. Numerical convergence is achieved with the characteristic element size $\Delta \vartheta =1^{\circ }$. Figure 3 shows a comparison of the numerical solution of (2.6) at $t=100$ with the analytical ones at orders ${O}(\vartheta ^{2})$ (solid lines) and ${O}(\vartheta ^{6})$ (dashed lines), which shows an overall agreement. The solution at order ${O}(\vartheta ^{6})$ gives a better agreement with the numerics in a larger range of $\vartheta$. For $c>1.2$, the numerical and analytical solutions start to deviate for $\vartheta > 60^{\circ }$. The agreement with the solution at second order is good in most cases for $\vartheta <60^{\circ }$. The second-order term in (3.6) vanishes when $c^{*}=\sqrt {2/3}\approx 0.81$. Under these conditions, the solution at ${O}(\vartheta ^{2})$ is constant along the zenith. The approximation at order ${O}(\vartheta ^{6})$ does not admit a constant solution. However, the minimum variation of its integral in the region $0<\vartheta <{\rm \pi} /2$, with respect to the constant value given by employing $H_0(t)$, is obtained for $c \approx 0.74$. Independently of the considered order of the solution, for very small (respectively very large) values of $c$ the thickness decreases (respectively increases) when moving away from the pole. Moreover, for $c< c^{*}$, the numerical solution and the analytical one at order ${O}(\vartheta ^{6})$ present a non-monotonic behaviour, as shown in figure 3 for $c=0.4,0.6$, with an initial decrease followed by a slight increase for $\vartheta > 70^{\circ }$. The solutions for $c>c^{*}$ monotonically increase.

Figure 3.

(a) Sketch of the spheroidal substrate with varying $c$. (b) Comparison of the numerical solution at $t=100$ of (2.6) (coloured dots) with the analytical ones at order ${O}(\vartheta ^{2})$ (solid lines) and ${O}(\vartheta ^{6})$ (dashed lines). Different colours identify different values of $c$: $c=0.4$ (blue), $c=0.6$ (orange), $c=0.8$ (yellow), $c=1$ (purple), $c=1.2$ (green), $c=1.4$ (cyan).

The leading-order large-time analytical solution presents a temporal decay $h \sim {t}^{-1/2}$. The spherical case is recovered by imposing $c=1$ (Couder et al. Reference Couder, Fort, Gautier and Boudaoud2005; Lee et al. Reference Lee, Brun, Marthelot, Balestra, Gallaire and Reis2016; Qin, Xia & Gao Reference Qin, Xia and Gao2021). The large-time solution is independent of the initial thickness $h_i$. It is interesting to note the good agreement between the analytical and numerical solutions for $c<1.2$ and $\vartheta >1$, which is out of the expected range of validity of the asymptotic expansion. The relative size of the terms in the asymptotic expansion decreases as higher orders are considered, thus suggesting that the power series expansion may converge to the exact solution in the considered range of $\vartheta$.

A decrease of $c$ implies a reduction of the gravity component parallel to the substrate and thus a reduction of the pole thickness, for a given time horizon. The film thinning or thickening moving downstream of the pole is also related to the considered geometry, as the consequence of two competing effects, already shown close to the pole (§ 2.3), where the thickness evolution depends on the normal gravity component multiplied by the local mean curvature $(\boldsymbol {g}\boldsymbol {\cdot } \boldsymbol {e}_3)\mathcal {K}$. While at the pole gravity is aligned with the substrate normal, i.e. $\boldsymbol {g}\boldsymbol {\cdot } \boldsymbol {e}_3=1$, as we move away the normal gravity component $\boldsymbol {g}\boldsymbol {\cdot }\boldsymbol {e}_3$ decreases with the zenith owing to the slope increase of the substrate, leading to a first inhomogeneity mechanism. Moving away from the pole, this slope increase leads to a slower decrease of the thickness with time. This explains why, in the case of constant curvature, e.g. spheres or cylinders, the thickness increases moving toward the equator, in agreement with the results of Lee et al. (Reference Lee, Brun, Marthelot, Balestra, Gallaire and Reis2016) and Balestra et al. (Reference Balestra, Kofman, Brun, Scheid and Gallaire2018a). The second mechanism at hand is the evolution of the curvature $\mathcal {K}$ along the zenith direction. Curvature variations induce an accumulation of fluid in regions of lower curvature, characterized by a slower decrease of the thickness with time. Spheroids with small height are characterized by a mean curvature that increases away from the pole, and vice versa for spheroids of large height. The former are thus likely to present a decreasing thickness moving downstream, and vice versa, as observed in figure 3. As a result, the thickness distribution is a result of the competition between variations of slope and mean curvature, which may induce thinning or thickening of the fluid layer.

Therefore, the transition does not occur when the mean curvature is constant (i.e. $c=1$), but when there is a balance between the thickness variations due to the change in mean curvature and those induced by slope variations. This value can be obtained by considering the quantity on the right-hand side of (2.9), i.e. $(\boldsymbol {g}\boldsymbol {\cdot }{\boldsymbol {e}}_3) \mathcal {K}$, which, in the vicinity of the pole, reads

(3.7)\begin{equation} ({\boldsymbol{g}}\boldsymbol{\cdot} {\boldsymbol{e}}_3) \mathcal{K} \approx{-}(2 c -c (3 c^{2}-2)\vartheta^{2}+{O}(\vartheta^{4})), \end{equation}

which is constant for $c=c^{*}$, i.e. the value that causes the ${O}(\vartheta ^{2})$ contribution to vanish. Note that the same transition value can be obtained by evaluating how the quantity $({1}/{w})\partial _1 (wg_t^{\{1\}})$ perturbs the ${O}(1)$ solution. The non-monotonic behaviours at large $\vartheta$ for $c< c^{*}$ are related to higher-order terms.

3.2. Spreading problem

Typical coating applications involve the spreading of an initial volume of fluid located close to the top of the considered geometry (Takagi & Huppert Reference Takagi and Huppert2010). Here, following previous works (Huppert Reference Huppert1982a), we recover some typical relevant quantities such as the position and thickness of the spreading front. An initial volume of fluid $V$ is released on the substrate. We impose the conservation of mass in general coordinates, under the assumption $\delta =0$ (Roy et al. Reference Roy, Roberts and Simpson2002):

(3.8)\begin{equation} \int\int_{S} {h}({x}^{\{1\}},{x}^{\{2\}},{t}) {w} \,\mathrm {d}\kern0.06em x^{\{1\}}\,\mathrm {d}\kern0.06em x^{\{2\}}=V, \end{equation}

where $V$ is the initial volume released on the substrate and $S$ is the region of the substrate, parameterized with $({x}^{\{1\}},{x}^{\{2\}})$, which contains the fluid and varies with time because of the moving front. Far from the contact line, the drainage solution ${h}({x}^{\{1\}},{x}^{\{2\}},t)$ can be employed, while capillary effects are relevant only in the close vicinity of the front (Huppert Reference Huppert1982a). For a fixed substrate geometry, ${w}$ is known, and thus relation (3.8) is an implicit equation with the front position as an unknown. A typical assumption to simplify the analysis is the employment of the large-time drainage solution.

We consider an initial volume of fluid of constant height $h_i=1$ released at $t=0$ in the region $0<\varphi < 2{\rm \pi}$, $0<\vartheta < \vartheta _0$. Owing to the invariance along the azimuthal direction, the conservation of the initial fluid volume (3.8) reads

(3.9)\begin{equation} \int_0^{\vartheta_F(t)} h(\vartheta,t)w(\vartheta)\,\mathrm{d}\vartheta= \int_0^{\vartheta_0} 1\,w(\vartheta)\,\mathrm{d}\vartheta, \end{equation}

where $\vartheta _F(t)$ is the front angle; the analytical expression (3.6) for $h(\vartheta,t)$ is employed. Equation (3.9) is numerically solved in Matlab via the built-in function ‘fsolve’. Figure 4(a) shows the evolution of the front angle $\vartheta _F$ with time, for different values of $\vartheta _0$ and $c$. An increase in $\vartheta _0$ leads to larger values of $\vartheta _F$, for fixed time. At small times, an increase in $c$ leads to larger $\vartheta _F$; however, at large times, the opposite behaviour is observed. In figure 4(b) we report the thickness at the front $h_F=h(\vartheta _F(t),t)$, which presents slight variations with $c$.

Figure 4.

Spreading of an initial volume of fluid on a spheroid. (a) Variation of the front angle $\vartheta _F$ with time and (b) of the thickness at the front $h_F$ with $\vartheta _F$, for different values of $c$ (coloured lines) and $\vartheta _0$ (different clusters of curves). The black dashed lines correspond to the analytical approximation of the relation $\vartheta _F(t)$ and $h_F(\vartheta _F)$, while the stars are the values recovered by a numerical simulation of the complete model with $c=0.6$, $Bo=500$, $\delta =10^{-3}$. (c) Numerical thickness distribution obtained from the complete model with $c=0.6$, $Bo=500$, $\delta =10^{-3}$ as a function of $\vartheta$ at different times: $t=20$ (blue), $t=40$ (orange), $t=60$ (yellow), $t=80$ (purple), $t=100$ (green), $t=120$ (cyan), $t=140$ (maroon), $t=160$ (blue). The black dashed lines denote the corresponding leading-order large-time drainage solution.

We approximate these results by considering an expansion for $\vartheta \ll 1$, by employing equation (3.4a,b) for $t \gg 1$, and $w(\vartheta )=\vartheta +{O}(\vartheta ^{2})$. In this case, both the right-hand side and left-hand side of (3.9) can be analytically integrated and an explicit relation for $\vartheta _F$ is found, together with an expression of the thickness at the front $h_F$

(3.10a,b)\begin{equation} \sqrt{\frac{3}{4{t}}}\frac{1}{\sqrt{c}}\frac{\vartheta_F^{2}}{2}=\frac{\vartheta_0^{2}}{2} \rightarrow \vartheta_F=\vartheta_0 \left(\frac{4ct}{3}\right)^{1/4}, \quad h_F=\left(\frac{\vartheta_0}{\vartheta_F}\right)^{2}. \end{equation}

Note that this expression with $c=1$ coincides with the solution on a sphere (Takagi & Huppert Reference Takagi and Huppert2010). These results, reported in the black dashed line in figure 4(a,b), agree well with the implicit equation for small values of $\vartheta$. The velocity of the front $U_F=\mathrm {d}\vartheta _F/\mathrm {d}t=({c}/{192})^{1/4}t^{-3/4}$ decreases with time. Therefore, the front slows down as moving downstream toward the equator, for all values of $c$.

We verify the faithfulness of this approach by comparing it with the numerical results of the complete model (2.4) with parameters $c=0.6$, $Bo=500$ and $\delta =10^{-3}$ (figure 4c). To simulate the spreading on the substrate, we consider a precursor film of size $h_{pr}=0.005$ (Troian, Wu & Safran Reference Troian, Wu and Safran1989b; Kondic & Diez Reference Kondic and Diez2002) with the following initial condition (Balestra et al. Reference Balestra, Badaoui, Ducimetière and Gallaire2019):

(3.11)\begin{equation} h(\vartheta,0)=\frac{h_i-h_{pr}}{2}\left(1-\tanh\left({100\left(\vartheta-\vartheta_0\right)}\right)\right)+h_{pr}.\end{equation}

Figure 4(c) shows the evolution with time of the film thickness, with $\vartheta _0=20^{\circ }$. In the vicinity of the front, a capillary ridge connects the film to the precursor one. Far from the front, the drainage solution approximates well the thin film evolution. In figure 4(a,b), we report also the position and the values of the maximum thickness at the ridge, with a good agreement with the analytical approach.

The spreading velocity decreases with time and is proportional to $c^{1/4}$, in the vicinity of the pole. As $c$ increases, for fixed $\vartheta _F$, the tangential gravity component increases while the area invaded by the fluid does not vary, at leading order ($w \approx \vartheta$), close to the pole. The propagation velocity therefore increases since a faster drainage is observed with increasing $c$. Nevertheless, at large times, spheroids with smaller $c$ present larger values of $\vartheta _F$. Close to the equator, the tangential gravity component is almost vertical and thus the film velocity is not strongly affected by $c$. Nevertheless, for fixed equatorial radius, the distance covered for a small increment $\mathrm {d} \vartheta _F$ increases with $c$, at large $\vartheta _F$, therefore implying a reduction of the spreading velocity $\mathrm {d}\vartheta _F/\mathrm {d}t$.

In this section, we described the competition between the substrate's slope and curvature in defining the drainage and spreading patterns on an axisymmetric substrate, the spheroid. In the ESM, we report also the case of a paraboloid, which instead always shows a decreasing mean curvature and thus an increasing thickness moving away from the pole. The spheroid analysis was simplified thanks to the absence of odd terms in the asymptotic expansion in $\vartheta$. To better understand the role of the curvature in modifying the drainage, we now consider the torus, a substrate in which the symmetry with respect to $\vartheta$ is broken.

4.1. Drainage problem

In this section, we consider the drainage of a thin film flowing on a toroidal substrate of tube radius $R$ and distance $dR$ between the axis of revolution and the centre of the tube (see figure 5a). The torus is thus generated by the rotation along the azimuthal direction of a circular cross-section whose centre is located at a distance $d$ from the axis of rotation. Non-dimensionalizing the in-plane directions and substrate variables with $R$, the following parameterization based on the zenith $\vartheta$ and the azimuth $\varphi$ is employed:

(4.1)\begin{equation} \boldsymbol{X}(\vartheta,\varphi)=((d+\sin\vartheta)\cos \varphi, (d+\sin \vartheta) \sin \varphi, \cos \vartheta). \end{equation}

The position along the cylinder, at each azimuthal circular cross-section, is defined through the zenith $\vartheta$. Two limiting cases are identified; the first one occurs for $d\rightarrow \infty$, which leads to the cylindrical case, reported in the ESM. The second case occurs for $d=1$, in which the points at $\vartheta =-90^{\circ }$ are in contact, leading to the so-called horn torus. The gravity term $g_t^{\{1\}}$ and $w$ read

(4.2a,b)\begin{equation} g_t^{\{1\}}(\vartheta)= \sin (\vartheta), \quad w(\vartheta)={d+\sin (\vartheta)}.\end{equation}

The same procedure employed for the drainage solution of the spheroidal case is adopted. However, in this case we cannot a priori neglect the odd terms in the asymptotic expansion, i.e. $h(\vartheta,t)=H_0(t)+\vartheta H_1(t)+\vartheta ^{2}H_2(t)\dots \,$. The resulting problems, at different orders in $\vartheta$, are reported in Appendix B. For the sake of brevity, the large-time solution at ${O}(\vartheta ^{4})$ reads

(4.3)\begin{align} h&=\sqrt{\frac{3}{2t}} \left(\frac{19\,377 \vartheta^{4}}{176\,000 d^{4}}-\frac{1409 \vartheta^{3}}{11\,000 d^{3}}-\frac{7477 \vartheta^{4}}{147\,840 d^{2}}+\frac{31 \vartheta^{2}}{200 d^{2}}\right. \nonumber\\ &\quad\left.+\frac{91 \vartheta^{3}}{2640 d}-\frac{\vartheta}{5 d}+\frac{43 \vartheta^{4}}{10\,752}+\frac{\vartheta^{2}}{16}+1\right)+{O}(\vartheta^{5})+{O}\left(\frac{1}{t^{3/2}}\right) . \end{align}

The cylinder thickness distribution is recovered for $d \rightarrow \infty$ (Balestra et al. Reference Balestra, Kofman, Brun, Scheid and Gallaire2018a). The ${O}(1)$ solution is analogous to the cylinder case for any value of $d$. The drainage problem is numerically solved in the domain $-{\rm \pi} /2<\vartheta <{\rm \pi} /2$. Numerical convergence is achieved with $\Delta \vartheta =0.5^{\circ }$. Figure 5 shows a comparison between the numerical and large-time analytical solutions of the drainage problem, for different values of $d$ in the range $-{\rm \pi} /2<\vartheta <{\rm \pi} /2$. The distribution is not symmetric with respect to $\vartheta =0$. In particular, the thickness is larger for negative values of $\vartheta$, i.e. on the inner side of the torus, while for $\vartheta >0$ the thickness is almost constant. These differences are enhanced as $d$ decreases. The numerical solution compares well with the analytical one at ${O}(\vartheta ^{4})$ while, at ${O}(\vartheta ^{2})$, the agreement is good only in the vicinity of the top.

Figure 5.

(a) Sketch of the axisymmetric flow configuration for the coating of a torus. (b) Drainage solution on a torus at $t=300$, numerical (coloured dots) and analytical solutions at order ${O}(\vartheta ^{2})$ (solid lines) and ${O}(\vartheta ^{4})$ (dashed lines), for $d=1.1$ (blue), $d=1.25$ (orange), $d=1.5$ (yellow), $d=2.5$ (purple), $d=5$ (green).

At the top ($\vartheta =0$), $\mathcal {K}=-1$ and therefore the film drains as in the cylinder case, locally. The different thickness distributions in the two sides of the circular cross-section of the torus result from the non-symmetric drainage with respect to $\vartheta$. While the slope is symmetric with respect to $\vartheta$, the mean curvature decreases on the inner part and decreases on the outer part. Following the discussion of § 3.1, a decreasing (respectively increasing) curvature induces an increasing (respectively decreasing) thickness. Therefore, much larger thicknesses are attained on the inner part than on the outer one, where the thickness slightly decreases, in the vicinity of the top. The slight increase on the outer part observed at large $\vartheta$ is due to the saturation of the mean curvature value, which remains almost constant, while the substrate's slope increases. From a quantitative point of view, we consider the product between the normal component of gravity and the mean curvature

(4.4)\begin{equation} \left(\boldsymbol{g}\,{\cdot}\, \boldsymbol{e}_3\right) \mathcal{K} \approx{-} \left(\left(-\frac{1}{a^{2}}-\frac{1}{2}\right) \vartheta^{2}+\frac{\vartheta}{a}+1\right), \end{equation}

which shows a decrease on the inner part, thus highlighting an accumulation of fluid downstream, and vice versa. The same result could be obtained by considering how the flux perturbs the ${O}(1)$ solution.

4.2. Spreading problem

We now present the results for the spreading of a volume of fluid on a torus. We consider an initial volume of fluid of thickness $h=1$ in the region $-\vartheta _0<\vartheta <\vartheta _0$. The breaking of symmetry with respect to $\vartheta =0$ results in two different spreading fronts for $\vartheta <0$ (inner side) and $\vartheta >0$ (outer side). However, at $\vartheta =0$, the drainage gravity component is exactly zero, i.e. $q^{\{1\}}=h^{3}g_t^{\{1\}}h^{3}=0$. Therefore, the total volume on each side of the torus is conserved since there is no flux at $\vartheta =0$. Note that, when hydrostatic or capillary effects are considered, the flux is not exactly zero at the top. A preliminary analysis showed that appreciable variations of the mass on the two sides (of the order of $2\,\%$) are observed for $Bo=250$ and $\delta =0.1$, when either pure capillary or pure hydrostatic effects are considered, in addition to drainage. For larger values of $Bo$ or smaller values of $\delta$, these differences rapidly decrease. In the limit $Bo\rightarrow \infty$ and $\delta =0$ (i.e. the considered drainage problem), a zero flux at the top of the torus is numerically observed.

The conservation of mass for the two regions reads

(4.5a,b)\begin{align} \int_0^{\vartheta^{O}_F(t)} h(\vartheta,t)w(\vartheta)\,\mathrm{d}\vartheta= \int_0^{\vartheta_0} w(\vartheta)\,\mathrm{d}\vartheta, \quad \int_{-\vartheta^{I}_F(t)}^{0} h(\vartheta,t)w(\vartheta)\,\mathrm{d}\vartheta= \int_{-\vartheta_0}^{0} w(\vartheta)\,\mathrm{d}\vartheta, \end{align}

where $\vartheta ^{O}_F(t)$ and $\vartheta ^{I}_F(t)$ are the front angle on the outer and inner part, respectively, $w(\vartheta )=d+\sin (\vartheta )$ and $h(\vartheta,t)$ is given by (4.3). Note that the two integrals on the right-hand side do not assume the same value, since $w(\vartheta )$ is not symmetric with respect to $\vartheta =0$. Equations (4.5a,b) are implicit integrals that are solved in Matlab through the built-in function ‘fsolve’. A first analytical approximation is found by taking the ${O}(1)$ approximation, leading to

(4.6)\begin{equation} \frac{\vartheta^{O}_F}{\vartheta_0}=\frac{\vartheta^{I}_F}{\vartheta_0}=\sqrt{\frac{2t}{3}}, \end{equation}

i.e. the solution of ${O}(1)$ does not depend on $d$ and is analogous to the spreading on a cylinder (Balestra et al. Reference Balestra, Kofman, Brun, Scheid and Gallaire2018a). The thickness at the front thus reads $h_F={\vartheta _0}/{\vartheta _F}$. A better approximation that includes the curvature of the torus can be obtained by considering the ${O}(\vartheta )$ approximation of the integrand

(4.7)\begin{gather} \left.\begin{gathered} \int_0^{\vartheta^{O}_F(t)} \sqrt{\frac{3}{2t}}\left(\frac{4}{5}\vartheta+d\right)\mathrm{d}\vartheta= \int_0^{\vartheta_0} (d+\vartheta)\,\mathrm{d}\vartheta\\ \rightarrow \left(\vartheta_F^O\right)^2 +\frac{5d}{2}\vartheta^{O}_F-\left(\frac{5}{2}({\rm d}\vartheta_0+\vartheta_0^{2}/2)\right)\sqrt{\frac{2t}{3}}=0, \\ \rightarrow \vartheta^{O}_F(t)=\frac{1}{2}\left(-\frac{5d}{2}+\sqrt{\frac{25d^{2}}{4}+4({\rm d}\vartheta_0+\vartheta_0^{2}/2)\frac{5}{2}\sqrt{\frac{2t}{3}}}\right), \end{gathered}\right\} \end{gather}

(4.8)\begin{gather} \left.\begin{gathered} \int_{-\vartheta^{I}_F(t)}^{0} \sqrt{\frac{3}{2t}}\left(\frac{4}{5}\vartheta+d\right)\mathrm{d}\vartheta= \int_{-\vartheta_0}^{0} (d+\vartheta)\,\mathrm{d}\vartheta\\ \rightarrow \left(\vartheta_F^I\right)^2-\frac{5d}{2}\vartheta^{I}_F +\left(\frac{5}{2}(d\vartheta_0-\vartheta_0^{2}/2)\right)\sqrt{\frac{2t}{3}}=0, \\ \rightarrow \vartheta^{I}_F(t)=\frac{1}{2}\left(\frac{5d}{2}-\sqrt{\frac{25d^{2}}{4}-4(d\vartheta_0-\vartheta_0^{2}/2)\frac{5}{2}\sqrt{\frac{2t}{3}}}\right). \end{gathered}\right\} \end{gather}

Figure 6(a,b) shows the behaviours of $\vartheta ^{O}_F$, $\vartheta ^{I}_F$ and the front thicknesses $h^{O}_F$ and $h^{I}_F$ on the inner and outer sides of the torus, respectively, for different values of $d$ and $\vartheta _0$. As concerns panel (a), for a fixed time, the front angle on the inner side is always larger than the one on the outer side. An increase of $d$ leads to a decrease (respectively increase) of $\vartheta _F$ on the inner (respectively outer) side. The front thickness does not strongly depend on $d$, even if some differences can be appreciated on the inner side, for large values of $\vartheta ^{O}_F$. The ${O}(1)$ approximation gives a reasonable agreement in the prediction of the front angle and thickness. In particular, it appears to be the lower (respectively upper) limit for the inner (respectively outer) sides, as $d$ increases. The order ${O}(\vartheta )$ approximations follow well the implicit relations (4.5a,b). We compare these analytical results with a numerical simulation of the complete model (2.4) with parameters $d=1.25$, $Bo=500$, $\delta =10^{-3}$, $h_{pr}=0.005$, initial condition

(4.9)\begin{gather} h(\vartheta,0)=\frac{h_i-h_{pr}}{2}\left(1-\tanh(100(\vartheta-\vartheta_0)\right))+h_{pr}, \quad \mathrm{for} \ \vartheta>0, \end{gather}

(4.10)\begin{gather} h(\vartheta,0)=\frac{h_i-h_{pr}}{2}\left(1-\tanh(100(-\vartheta-\vartheta_0)\right))+h_{pr}, \quad \mathrm{for} \ \vartheta<0, \end{gather}

and $\vartheta _0=10^{\circ }$ (see figure 6c). The agreement between the numerical front angle, given by the maximum thickness location, and the theoretical one is very good, and also the maximum thickness well follows the front thickness predicted by the theory.

Figure 6.

Spreading of an initial volume of fluid on a torus. (a) Variation of the front angle $\vartheta _F$ with time and (b) of the thickness at the front $h_F$ with $\vartheta _F$, for different values of the initial angle $\vartheta _0$ and $d$. The solid and dot-dashed lines denote the values of $\vartheta _F$ and $h_F$ on the outer and inner sides, respectively. The black and red dashed lines correspond to the ${O}(1)$ and ${O}(\vartheta )$ analytical approximations of the relation $\vartheta _F(t)$ and $h_F(\vartheta _F)$, respectively, while the stars are the values recovered by a numerical simulation of the complete model with $d=1.25$, $Bo=500$, $\delta =10^{-3}$, precursor film thickness $h_{pr}=0.005$. (c) Numerical thickness distribution obtained from the complete model with $d=1.25$, $Bo=500$, $\delta =10^{-3}$ as a function of $\vartheta$ at different times: $t=10$ (blue), $t=20$ (orange), $t=30$ (yellow), $t=40$ (purple). The black dashed lines denote the corresponding large-time drainage solutions.

In analogy with the drainage solution, the faster spreading attained on the inner region is related to the substrate geometry. For a fixed angular distance from the top, the area covered by the spreading fluid is larger on the outer region than on the inner one. Therefore, for a fixed time, the fluid spreads faster on the inner region, reaching larger values of $\vartheta _F$ than on the outer region. Interestingly, the solution at ${O}(1)$ does not capture the symmetry breaking, since, at the top, a torus locally coincides with a cylinder. Nevertheless, the ${O}(\vartheta )$ approximation already captures the asymmetry of the substrate.

The torus case shows non-symmetric drainage and spreading along the zenith direction. In the following, we present how these analyses can be extended to non-axisymmetric substrates which are characterized by a three-dimensional, non-uniform along the azimuthal direction, drainage. We chose as a testing ground an ellipsoid with three different axes.

5.1. Numerical drainage solution

In this section, we study the coating of an ellipsoidal substrate of horizontal semiaxes $aR$, $bR$ and vertical semiaxis $R$ (see figure 1); gravity is pointing downward. In non-dimensional form, the following parameterization holds:

(5.1)\begin{equation} \boldsymbol{X}(\vartheta,\varphi)=(a\sin\vartheta\cos \varphi,b \sin \vartheta \sin \varphi,\cos \vartheta). \end{equation}

We identify different limiting cases, depending on the values of $a$ and $b$. If $a=b=1$, we recover the spherical case; if $a=b \neq 1$ the resulting substrate is an axisymmetric ellipsoid of unitary height and equatorial radius $a=b$, whose results can be recovered from those of § 3.1. Note that the time scale is different since the in-plane directions and substrate variables are non-dimensionalized with the height and not with the equatorial radius, in the present section. When $a \neq b$ the axisymmetry is broken since the two axes at the equator are different. In the following, we assume that $b \geqslant a$ and consider the range $0.4< a,b<2$. Note that the solutions for $a>b$ can be recovered by simply translating of $\varphi =90^{\circ }$ the solution for $b \geqslant a$ (obtained by swapping the desired values of $a$ and $b$).

The metric components vary along the $\varphi$ direction; in particular, the metric tensor is not diagonal (see the ESM). The local coordinates system defined by the parameterization is thus non-orthogonal and a second gravity component $g_t^{\{2\}}(\vartheta,\varphi )$, appears. The square root of the determinant of the metric and the gravity terms now read

(5.2)$$\begin{gather} w(\vartheta,\varphi)=\sqrt{\sin ^{4}(\vartheta) \left(a^{2} \sin ^{2}(\varphi)+b^{2} \cos ^{2}(\varphi)\right)+a^{2} b^{2} \sin ^{2}(\vartheta) \cos ^{2}(\vartheta)}, \end{gather}$$

(5.3)$$\begin{gather}g_t^{\{1\}}(\vartheta,\varphi)=\frac{\sin (\vartheta) \left(a^{2} \sin ^{2}(\varphi)+b^{2} \cos ^{2}(\varphi)\right)}{\sin ^{2}(\vartheta) \left(a^{2} \sin ^{2}(\varphi)+b^{2} \cos ^{2}(\varphi)\right)+a^{2} b^{2} \cos ^{2}(\vartheta)}, \end{gather}$$

(5.4) $$\begin{gather}g_t^{\{2\}}(\vartheta,\varphi)= \frac{\sin (\varphi) \cos (\varphi) \left(a^{2} \cos (\vartheta)-b^{2} \cos (\vartheta)\right)}{\begin{array}{l}a^{2} b^{2} \cos ^{2}(\vartheta) \cos ^{2}(\varphi)+a^{2} b^{2} \cos ^{2}(\vartheta) \sin ^{2}(\varphi)+a^{2} \sin ^{2}(\vartheta) \sin ^{2}(\varphi)+b^{2} \sin ^{2}(\vartheta) \cos ^{2}(\varphi)\end{array}}. \end{gather}$$

We solve (2.6) by imposing periodic boundary conditions in $0<\varphi <2{\rm \pi}$ and the initial condition $h(\vartheta,\varphi,0)=1$. Numerical convergence is achieved with a characteristic mesh size of $0.9^{\circ }$. Figures 7 and 8 respectively show the resulting film distributions and a section at $\vartheta ={\rm \pi} /4$ for different values of $a$ and $b$, at $t=100$. We first increase the value of $b$, with $a=1$. For $b=1.2$ (panel (a)), the thickness presents modulations along the azimuthal direction, with a maximum thickness localized at $\varphi =k{\rm \pi}$ ($k=0,1,2$), i.e. along the direction of the smaller axis $a$. These modulations are enhanced as $b$ increases (panel (b)), with larger values of the attained thickness. Two regions of low thickness are localized at $\varphi ={\rm \pi} /2+k{\rm \pi}$, along the larger axis $b$. The same trends are observed further increasing $b$ (panel (c)). When $b=1$ and $a$ decreases (panel (d)), the thickness also presents modulations along the azimuthal direction, but the thickness always increases moving downstream. The thickness decreases as $a$ decreases. Similar patterns are also obtained when small values of $a$ and large values of $b$ are considered. The numerical solution of (2.6) shows the presence of modulations of the thickness along the azimuthal direction. According to § 3.1, spheroids with small (respectively large) height were characterized by a decrease (respectively increase) of the thickness. We can extend these considerations to an ellipsoid by considering the drainage along the principal directions defined by $(x,y)$, see figure 1(c). Since the drainage component along the azimuthal direction is identically zero along the two principal semiaxes, the flow locally behaves like the spheroidal case of § 3.1. Therefore, we expect to follow these trends along the two semiaxes, depending on $a$ and $b$. In the axisymmetric case, the thickness increases downstream for height–radius ratios larger than $0.74$, which corresponds to $a,b \lessapprox 1.35$. Therefore, when $a,b\lessapprox 1.35$ we always observe an increase of the thickness with $\vartheta$, as observed in figure 9(a–c) (see also figure 3 for the cases with $c>c^{*}$). However, the thickness presents clear modulations owing to the non-uniform drainage when $a \neq b$. Similarly, when $a,b\gtrapprox 1.35$ one expects a decrease of the thickness followed by a slight increase at large $\vartheta$, with modulations if $a \neq b$, as shown in figure 9(d–g). The intermediate situation occurs when $a \lessapprox 1.35$ and $b \gtrapprox 1.35$, characterized by an increase of the thickness along the $x$ direction and a decrease along the $y$ direction, as observed in figure 7(b,c,g).

Figure 7.

Numerical solution of (2.6) at $t=100$ as a function of $(\vartheta,\varphi )$, for different values of the semiaxes $a$ and $b$, with $a\leqslant 1$ and $b \geqslant 1$.

Figure 8.

Numerical solution of (2.6) at $t=100$ and $\vartheta ={\rm \pi} /4$ as a function of $(\varphi )$: (a) $a=1$ and $b=1$ (blue), $b=1.2$ (orange), $b=1.4$ (yellow), $b=1.6$ (purple), $b=1.8$ (green), $b=2$ (cyan); (b) $b=1$ and $a=0.4$ (blue), $a=0.6$ (orange), $a=0.8$ (yellow), $a=1$ (purple).

Figure 9.

Numerical solution of (2.6) at $t=100$ as a function of $(\vartheta,\varphi )$, for different values of the semiaxes $a$ and $b$, with $a,b<1$ and $a,b>1$.

The modulations of the thickness distribution are related to the variation of drainage with the azimuth. In the vicinity of the minor semiaxis, the tangential gravity component along the zenith is larger than close to the major semiaxis. Higher velocities are thus attained along the minor semiaxis, displacing more fluid downstream than along the major semiaxis. This process induces transport of fluid from progressively farther and farther regions and thus a secondary flow from the major semiaxis (associated with low velocities) to the minor semiaxis (associated with large velocities). In the light of this discussion, one may wonder if these patterns persist with time or merely represent a snapshot of a more intricate evolution. In the following, we also aim at clarifying this aspect by deriving an analytical solution for the drainage problem.

5.2. Analytical drainage solution

In this section, we derive an analytical drainage solution and compare it with the numerical results of the previous section. In analogy with § 3.1, we perform an asymptotic expansion in powers of $\vartheta$, with $\vartheta \ll 1$. The solution at order ${O}(1)$ does not depend on $\varphi$ since the solution at the pole has to be unique. We thus consider the following expansion, in which the odd terms have been removed because of symmetry:

(5.5)\begin{equation} h(\vartheta,\varphi,t)=H_0(t)+\vartheta^{2}H_2(\varphi,t)+\vartheta^{4}H_4(\varphi,t)+\vartheta^{6}H_6(\varphi,t)+\dots . \end{equation}

We expand (2.6) at various orders in $\vartheta$. At order ${O}(1)$, one obtains the following ordinary differential equation (ODE):

(5.6)\begin{align} \frac{1}{3} \left(\frac{1}{a^{2}}+\frac{1}{b^{2}}\right) H_0(t)^{3}+H_0'(t)=0 \rightarrow H_0(t)=\frac{1}{\sqrt{\dfrac{2}{3} t \left(\dfrac{1}{a^{2}}+\dfrac{1}{b^{2}}\right)+1}}=\frac{1}{\sqrt{\alpha t +1}}, \end{align}

where $\alpha =\frac {2}{3}(1/a^{2}+1/b^{2})$. Also in this case, the ${O}(1)$ solution reduces to $({3}/(2\mathcal {K}_pt))^{1/2}$ at late time, with $\mathcal {K}_p=(1/a^{2}+1/b^{2})$. The equation at order ${O}(\vartheta ^{2})$ reads

(5.7) \begin{align} &\frac{\partial H_2(\varphi,t)}{\partial t}\nonumber\\ &\quad =\frac{H_0(t)^{2} ((b^{2}-a^{2}) \sin (2 \varphi) \dfrac{\partial H_2}{\partial \varphi}+2 H_2 ((a^{2}-b^{2}) \cos (2 \varphi)-2 (a^{2}+b^{2})))}{2 a^{2} b^{2}}\nonumber\\ &\qquad+\frac{\begin{array}{l}H_0(t)^{3} ((a^{4} (b^{2}-2)-a^{2} b^{4}+2 b^{4}) \cos (2 \varphi) +2 (a^{4} (-(b^{2}-1))\\\quad+a^{2} (b^{2}-b^{4})+b^{4}))\end{array}}{6 a^{4} b^{4}}, \end{align}

which is a parabolic partial differential equation (PDE) in $H_2(\varphi,t)$. We numerically solve (5.7) with initial condition $H_2(\varphi,0)=0$. The periodic boundary conditions at $\varphi = [0,2{\rm \pi} ]$ are automatically imposed thanks to a Fourier spectral collocation method implemented in Matlab. The time stepping is performed by employing the built-in function ‘ode23t’, with a tolerance of $10^{-6}$. Numerical convergence is achieved with $100$ collocation points.

Figure 10(a) shows the spatio-temporal evolution of the second-order solution $H_2(\varphi,t)$, for $a=0.5$ and $b=1.5$. An initial growth in absolute value until $t \approx 0.3$ is followed by a slow decay at large times. In figure 10(b) we report the $H_2$ profiles rescaled with $H_0$, at different times in the slow-decay regime. The second-order solution $H_2$ is ${\rm \pi}$-periodic and the maximum is attained at $\varphi =k{\rm \pi}$, i.e. along the smaller axis of the ellipsoid. At $\varphi =k{\rm \pi} /2$, i.e. along the larger axis of the ellipsoid, the correction reaches much smaller values. As time increases, the profiles collapse on a single curve, suggesting that a large-time solution characterized by a separation of variables is possible, i.e. $H_2(\varphi,t)=H_0(t)H_2^{*}(\varphi )$. We introduce this decomposition in (5.7). Exploiting (5.6), the temporal dependence disappears and the following ODE for $H_2^{*}(\varphi )$ is obtained:

(5.8) \begin{align} &-a^{4} b^{2} \cos (2 \varphi)+2 a^{4} b^{2}+2 a^{4} \cos (2 \varphi)-2 a^{4}+a^{2} b^{4} \cos (2 \varphi)+2 a^{2} b^{4}\nonumber\\ &\quad +3 a^{2} b^{2} (a^{2}-b^{2}) \sin (2 \varphi) H_2^{*\prime}(\varphi) -2 a^{2} b^{2} H_2^{*}(\varphi) (3 (a^{2}-b^{2}) \cos (2 \varphi)-5 (a^{2}+b^{2}))\nonumber\\ &\quad -2 a^{2} b^{2}-2 b^{4} \cos (2 \varphi)-2 b^{4}=0, \end{align}

whose solution reads

(5.9) \begin{align} H_2^{*}(\varphi)&=C_1 \sin (2 \varphi) \sin ^{-{5 (a^{2}+b^{2})}/{3 (a^{2}-b^{2})}}(\varphi) \cos ^{{5 (a^{2}+b^{2})}/{3 (a^{2}-b^{2})}}(\varphi)\nonumber\\ &\quad -\frac{1}{8 a^{2} b^{2} (4 a^{2}+b^{2}) (a^{2}+4 b^{2})}(a^{6} (7 b^{2}-4) +26 a^{4} (b^{4}-b^{2})\nonumber\\ &\quad +a^{2} b^{4} (7 b^{2}-26)+(a^{2}-b^{2}) (a^{4} (b^{2}+4)+a^{2} b^{2} (b^{2}+14)+4 b^{4})\nonumber\\ &\quad \times \vphantom{(a^{4} (b^{2}+4)+a^{2} b^{2} (b^{2}+14)+4 b^{4})}\cos (2 \varphi)-4 b^{6}), \end{align}

where $C_1$ is a constant to be determined. However, it is observed that $C_1 \neq 0$ implies an unbounded behaviour. Therefore, we impose $C_1=0$ to prevent non-physical solutions. The analytical result for $H_2^{*}$ is reported in figure 10(b), with an excellent agreement with the numerical solution. We then investigate the effect of $a$ and $b$ by considering different cases at $t=100$, reported in figure 11. An increase of $b$ for fixed $a=0.5$ (panel (a)) leads to a decrease of $H_2^{*}$ in the region $\varphi =k{\rm \pi} /2$, while an increase in $a$ for fixed $b=1.5$ (panel (b)) shows an overall decrease of $H_2^{*}$. Also for these cases, an excellent agreement with the analytical solution is observed. At ${O}(\vartheta ^{2})$, the large-time analytical solution can be written in compact form as

(5.10)\begin{equation} h \approx H_0(t) (1+ (f(a,b) +g(a,b)\cos(2\varphi))\vartheta^{2}).\end{equation}

The modulations observed in the numerical simulations of the previous section are captured by the ${O}(\vartheta ^{2})$ term, which is a ${\rm \pi}$-periodic function of the azimuth. These modulations are present as long as $g(a,b)=(a^{2}-b^{2}) (a^{4} (b^{2}+4)+a^{2} b^{2} (b^{2}+14)+4 b^{4} ) \neq 0$. The only case in which modulations are absent occurs when $g(a,b)=0$ and thus $a=b$, i.e. the spheroidal case. In the latter case, the solution reads $H_2^{*}=\frac {1}{10}({3}/{b^{2}}-2)$ and is formally analogous to the second-order solution of the spheroid with $c=1/b$ (see § 3.1).

Figure 10.

Drainage on an ellipsoid with $a=0.5$ and $b=1.5$. (a) Spatio-temporal evolution of $H_2$: iso-contours of $H_2$ in the $(\varphi,t)$ plane. (b) Second-order correction $H_2^{*}=H_2/H_0\approx H_2\sqrt {\alpha t}$ as a function of $\varphi$ at different times: $t=0.4$ (blue), $t=1$ (orange), $t=5$ (yellow), $t=10$ (purple) $t=30$ (green), $t=50$ (cyan), $t=70$ (maroon), $t=90$ (black), $t=100$ (red). The black stars denote the late-time analytical solution for $H_2^{*}$ from (5.9).

Figure 11.

(a) Second-order correction $H_2^{*}=H_2/H_0\approx H_2\sqrt {\alpha t}$ as a function of $\varphi$ at $t=100$, for $a=0.5$ and increasing b: $b=0.6$ (blue), $b=0.8$ (orange), $b=1$ (yellow), $b=1.2$ (purple), $b=1.4$ (green), $b=1.6$ (cyan), $b=1.8$ (maroon), $b=2$ (black); $(b)$ $H^{*}_2$ as a function of $\varphi$ at $t=100$, for $b=1.5$ and increasing a: $a=0.4$ (blue), $a=0.6$ (orange), $a=0.8$ (yellow), $a=1$ (purple), $a=1.2$ (green), $a=1.4$ (cyan). The black stars denote the late-time analytical solution for $H_2^{*}$ from (5.9).

The faithfulness of the analytical solution is verified against the numerical simulations of § 5.1 in figure 12. For the comparison, we consider the solution at orders ${O}(\vartheta ^{2})$ and ${O}(\vartheta ^{6})$. The higher-order problems, together with their solutions $H_4$ and $H_6$, are reported in Appendix C. The same large-time behaviour is observed. In general, the analytical solutions at ${O}(\vartheta ^{6})$ compare well with the numerical ones, while those at ${O}(\vartheta ^{2})$ are accurate only in the vicinity of the pole. The analytical solution at ${O}(\vartheta ^{6})$ deviates from the numerical one for $a<0.8$. The agreement for $a>0.8$ is satisfactory for any value of $b$.

Figure 12.

Comparison at three different azimuthal sections between the numerical (coloured dots) and the quasi-analytical solution for an ellipsoid at ${O}(\vartheta ^{2})$ (solid lines) and ${O}(\vartheta ^{6})$ (dashed lines) at $t=100$, for different values of $a$ and $b$.

In this section, we derived an analytical approximation for the drainage problem. The problem was solved by employing an asymptotic expansion in a first stage, followed by a separation of variables at each order of the expansion. The final structure of the analytical approximation (5.10) is characterized by a time dependence separated by the spatial one, similarly to the previous cases. Nevertheless, the power-series expansion presents terms that depend on the azimuth. The simple form of (5.10) captures well the ${\rm \pi}$-periodicity of the drainage solution, induced by the differences in drainage along the minor and major axes. In the spreading problem, these modulations may play a crucial role.

5.3. Spreading problem

In this section, we consider the spreading of an initial volume of fluid of height $h_i=1$ contained in the region $0<\vartheta <\vartheta _0$, $0<\varphi <2{\rm \pi}$. Figure 13 shows the evolution of the film thickness with time, for different values of $a$ and $b$, obtained employing the complete model (2.4) with initial condition formally analogous to (3.11), i.e. invariant along the azimuthal direction. For $t=10$, the maximum thickness position, at which the front is located, is modulated along the azimuthal direction. This modulation accentuates with time and a region of large thickness forms at $\varphi =k{\rm \pi}$ (along the shorter axis) while the thickness is much lower at $\varphi =k{\rm \pi} /2$. Therefore, the front presents two peaks of large thickness aligned along the shorter axis. This effect is enhanced when larger (respectively lower) values of $b$ (respectively $a$) are considered.

Figure 13.

Iso-contours of the numerical solution for the spreading of (2.4), at different times, for an ellipsoid and $Bo=100$, $\delta =10^{-2}$, precursor film $h_{pr}=0.02$ and $\vartheta _0=20^{\circ }$; (a) $a=0.8$, $b=1$, (b) $a=0.8$, $b=1.6$, (c) $a=1$, $b=1.6$, (d) $a=1.4$, $b=1.8$.

In figure 14, we report a zoom in the region $0<\varphi <{\rm \pi}$ for one simulation of the complete model (2.4) with $Bo=1000$ and $\delta =10^{-2}$, together with a three-dimensional rendering of the spreading on the ellipsoid, viewed from the top. The black lines denote the streamlines of the flux $\boldsymbol {q}=q^{\{1\}}\boldsymbol {e}_1+q^{\{2\}}\boldsymbol {e}_2$. The flux streamlines are almost parallel to the azimuthal direction at low values of $\vartheta$, then bend and align along the zenith direction as $\vartheta$ increases. An exception to this behaviour is observed at $\varphi =0,{\rm \pi} /2$, in which the flow streamlines are always parallel to the zenith direction. The three-dimensional rendering highlights the formation of two, finger-like, front peaks of large thickness along the shorter axis, while the fluid slowly spreads along the larger axis.

Figure 14.

Iso-contours of the numerical spreading solution at different times of (2.4) for an ellipsoid with $a=1$, $b=1.4$, $Bo=1000$, $\delta =10^{-2}$, precursor film $h_{pr}=0.02$ and $\vartheta _0=10^{\circ }$. The black solid lines denote the streamlines of the volume flux per unit length $\boldsymbol {q}$ and the white line the average front $\bar {\vartheta }_F$.

A scaling law for the spreading front and thickness is obtained by neglecting the modulations of the front, and assuming a constant average value along the azimuth, i.e. $\bar {\vartheta }_F=1/(2{\rm \pi} )\int _0^{2{\rm \pi} }\vartheta _F(\varphi,t)\,\mathrm {d}\varphi$. The conservation of volume reads

(5.11)\begin{equation} \int_0^{2{\rm \pi}} \int_0^{\bar{\vartheta}_F(t)} h(\vartheta,\varphi,t)w(\vartheta)\,\mathrm{d}\vartheta\,\mathrm{d}\varphi= \int_0^{2{\rm \pi}} \int_0^{\vartheta_0} w(\vartheta)\,\mathrm{d}\vartheta\,\mathrm{d}\varphi, \end{equation}

where $h=H_0(t)+\vartheta ^{2}H_2(\varphi,t)+\vartheta ^{4}H_4(\varphi,t)+\vartheta ^{6}H_6(\varphi,t)$ is the asymptotic solution obtained in the previous section, and $w$ is given by (5.2). Also in this case, an analytical approximation is found by employing the large-time ${O}(\vartheta )$ approximation (5.6), with $w=ab\vartheta +{O}(\vartheta ^{2})$, leading to the following expressions for the average front position and thickness:

(5.12a,b)\begin{equation} \bar{\vartheta}_F=\vartheta_0\left(\frac{2}{3}\left(\frac{1}{a^{2}}+\frac{1}{b^{2}}\right)t\right)^{1/4}, \quad \bar{h}_F=\left(\frac{\vartheta_0}{\bar{\vartheta}_F}\right)^{2}. \end{equation}

The azimuth-averaged numerical solution of (5.11) and the analytical approximation (5.12a,b) are reported in figure 15, displaying a good agreement for low values of $\vartheta _0$ and large values of $a$, while the results start to diverge for large $\vartheta _0$ and small $a$.

Figure 15.

Spreading of an initial volume of fluid on an ellipsoid. $(a)$ Variation of the average front angle $\bar {\vartheta }_F$ with time and $(b)$ of the average thickness at the front $\bar {h}_F$ with $\bar {\vartheta }_F$, for $a=0.6$ (dash-dotted lines), $a=1$ (solid lines) and different values of the initial angle $\vartheta _0$ and $b$. The black dashed lines correspond to the analytical approximation of the relation $\bar {\vartheta }_F(t)$ and $\bar {h}_F(\vartheta _F)$, while the stars are the values recovered by a numerical simulation of the complete model with $a=1$, $b=1.4$, $Bo=1000$, $\delta =10^{-2}$, precursor film $h_{pr}=0.02$ and $\vartheta _0=10^{\circ }$.

Figure 16 shows the evolution of the front position and thickness with time, obtained from a numerical simulation of the complete model with $a=1$, $b=1.4$, $Bo=1000$, $\delta =10^{-2}$, precursor film $h_{pr}=0.02$ and $\vartheta _0=10^{\circ }$. The values of front position and thickness are averaged and compared with the analytical prediction. The analytical and numerical simulation results show similar trends. However, at large times, the modulations of the front are very large and the front travels much faster along the shorter axis than along the longer one.

Figure 16.

Maximum thickness (a) position and (b) value recovered from the numerical spreading simulation with $a=1$, $b=1.4$, $Bo=1000$, $\delta =10^{-2}$, precursor film $h_{pr}=0.02$ and $\vartheta _0=10^{\circ }$. Different colours denote different times $10 \leqslant t \leqslant 120$, with step size $\Delta t=10$. In the insets, we report a comparison between theoretical (black dashed line) and numerical (red dots) (a) average front position and (b) average front thickness.

The spreading problem on an ellipsoid is characterized by a different front speed along the azimuthal direction, which leads to an accumulation of fluid and a faster spreading along the smaller axis. Peaks of large thickness form together with a modulation of the front, prior to any fingering instability. These modulations are similar to those observed in the previous section for the drainage solution. As already explained, larger velocities induce transport of fluid from regions of lower velocity to regions of larger velocity. As a result, fluid accumulates and forms the observed peaks of large thickness. These velocity differences lead to a progressively more pronounced bending of the front. Therefore, a fingering instability analysis necessarily needs to consider the non-uniform spreading of the fluid along the ellipsoid, which may lead to the preferential formation of fingers. While this analysis focused on the spreading in the absence of surface tension, further studies may involve the formation of fingers resulting from the driven contact-line instability.

In this section, we described the drainage and spreading solution for the coating on an ellipsoid. We obtained an analytical solution that compares well with the numerical one. We showed the potential of general coordinates and asymptotic expansions to obtain a two-dimensional analytical solution suitable for a physical interpretation of the drainage and spreading process, in complement to the previous results for axisymmetric geometries. A large-time solution characterized by the separation of temporal and the two spatial dependencies was obtained. Modulations of the drainage solution and spreading front were explained in terms of the different slopes along the principal semiaxes, which induce an accumulation of fluid along the minor axis.