1. Introduction
Consider the behaviour of a thin viscous fluid layer coating a solid substrate. For a substrate with constant curvature, such as a sphere or cylinder of infinite extent, the effect of surface tension is to smooth out the fluid surface over time, resulting in a fluid layer of uniform depth over the substrate. However, when the substrate curvature varies, such as for an ellipsoid or an oblique cylinder, surface tension induces fluid movement so as to produce a non-uniform fluid depth – the precise final configuration then depends on the specified substrate. Further, the effect of shear, induced from e.g. a free-stream air flow, is directional. Shear induces movement of fluid across the substrate, which results in a shallower fluid layer at the ‘front’ of the substrate, and a deeper fluid layer at the ‘rear’. The addition and removal of fluids via e.g. mass injection can also impact the onset of instabilities. The main aim of this work is to consider these competing effects of shear, surface tension and mass injection for the prototypical scenario of thin-film flow on an ellipse.
Shear and surface tension are two crucial effects in the free-surface flow of viscous fluids over topography. For instance, to eliminate defects in industrial coating processes, where fluid is transferred to a moving substrate, the surface tension of the fluid must be carefully controlled (Weinstein & Ruschak Reference Weinstein and Ruschak2004; Ding et al. Reference Ding, Liu and Harris2016). Surface tension is also responsible for preventing the rupture and de-wetting of thin fluid layers in biological applications, such as in pulmonary airways (Grotberg Reference Grotberg1994) and on the surface of the eye (Braun Reference Braun2012). The effect of wind shear is rarely studied in isolation, and is typically paired with body forces, such as those induced by gravity or electric fields. For instance, for a fluid layer coating an inclined substrate, it is possible to obtain stationary surface profiles due to the competing effects of wind shear and gravity (King & Tuck Reference King and Tuck1993; Meng et al. Reference Meng, Papageorgiou and Vanden-Broeck2024). Analogous results have also been obtained in this manner by Tseluiko & Kalliadasis (Reference Tseluiko and Kalliadasis2011) for the flow of liquid on the interior of a cylindrical surface. Further, in the flow of a fluid partially coating an inclined substrate, rivulet patterns are known to emerge from a fingering instability at the fluid–solid contact line (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011); the resultant rivulet patterns can then be stabilised by wind shear (Saber & El-Genk Reference Saber and El-Genk2004), and have also been observed on the exterior of curved substrates by Takagi & Huppert (Reference Takagi and Huppert2010) and Mayo et al. (Reference Mayo, McCue and Moroney2013), and on the interior by Hammoud (Reference Hammoud2016).
When the depth of the fluid is small with respect to a horizontal length scale, the lubrication approximation can be employed to derive the thin-film equation (Oron et al. Reference Oron, Davis and Bankoff1997; Myers Reference Myers1998; Craster & Matar Reference Craster and Matar2009) which governs the fluid surface, and indeed a significant number of works have previously investigated such thin-film flows over curved substrates. For instance, since the pioneering analytical and experimental study on fingering instabilities by Moffatt (Reference Moffatt1977), the flow surrounding rotating cylinders has been extensively studied. Time-dependent solutions were obtained for this formulation through the method of multiple scales by Hinch & Kelmanson (Reference Hinch and Kelmanson2003). Karabut (Reference Karabut2007), Thiele (Reference Thiele2011) and Lopes et al. (Reference Lopes, Thiele and Hazel2018) considered thin-film flows on rotating cylinders under the effects of surface tension and gravity. It was shown that bifurcations occurred at specific ranges of rotational speeds. The effect of shear has also been included in this problem by Newell & Viljoen (Reference Newell and Viljoen2019) and Mitchell et al. (Reference Mitchell, Duffy and Wilson2022), who considered the inclusion of air with laminar flow in the far field, and mass flux was incorporated into this formulation by Noakes et al. (Reference Noakes, King and Riley2011). Further, the behaviour of thin fluid films on rotating substrates with varying curvature has also been studied; a rotating ellipse was considered by Li et al. (Reference Li, Carvalho and Kumar2017), more general curved substrates by Parrish et al. (Reference Parrish, Carvalho and Kumar2022) and Shepherd et al. (Reference Shepherd, Boujo and Sellier2024) and an ellipsoid in three dimensions was studied by Duruk et al. (Reference Duruk, Boujo and Sellier2021). Substrate curvature is also known to suppress instabilities for the flow on the interior of a cylinder (Trinh et al. Reference Trinh, Kim, Hammoud, Howell, Chapman and Stone2014). The effect of varying substrate curvature has also been studied for non-rotating substrates by Kalliadasis et al. (Reference Kalliadasis, Bielarz and Homsy2000) and Kalliadasis & Homsy (Reference Kalliadasis and Homsy2001), who obtained steady solution profiles with capillary ridges induced by the substrate and studied the resultant stability of these. Additionally, Schwartz & Weidner (Reference Schwartz and Weidner1995) studied a thin-film flow coating a corner, in which the fluid was seen to eventually de-wet close to the location of highest curvature.
The effect of shear has widely been considered in thin-film flows with solidification, a phase change from water to ice, for aeroplanes flying at high altitude (Myers et al. Reference Myers, Charpin and Chapman2002). The shear strength is large due to the high air velocity involved. Here, ice can form from the thin-film flow on cold substrates, such as aeroplane wings. More recently, ice formation within aeroplane engines has received extensive numerical and experimental study, due to the subsequent power loss that can occur when the ice detaches and damages critical engine components. For this latter case in aeroplane turbine engines, known as ice crystal icing (Beaugendre et al. Reference Beaugendre, Morency and Habashi2003), the substrate is heated. The thin layer of viscous fluid exists as a stage prior to the formation of ice, which occurs due to the accretion of small airborne ice particles. This has many similarities to our present study. For instance, the solid substrate, such as an engine stator or rotor blade, has substantial curvature variation with two regions of high curvature and two regions of low curvature. The ellipse studied in this work has very similar properties to this. The aim is to understand the water movement, temperature and freezing time (Villedieu et al. Reference Villedieu, Trontin and Chauvin2014; Wright et al. Reference Wright, Struk, Bartkus and Addy2015). More recent models (Currie Reference Currie2020; Connolly Reference Connolly2021) consider only the effect of shear when determining the thin film; some recent work has applied numerical simulations to study films under a wider range of operating conditions (Lallement et al. Reference Lallement, Trontin, Laurent and Villedieu2018; Chen et al. Reference Chen, Yi, Wang and Liu2023) but more theoretical analysis is required. This present work is strongly motivated by the authors’ recent research on the development of thin-film models for ice crystal icing (Peters et al. Reference Peters, Shelton, Tang and Trinh2024) where it was recognised that fundamental analysis of shear and surface-tension-driven thin films on curved surfaces was required.
In this work, we consider the formulation for a viscous fluid surrounding a two-dimensional elliptical substrate of varying curvature, with the effects of both shear and surface tension included. Motivated by the prior application in aeroengines (Peters et al. Reference Peters, Shelton, Tang and Trinh2024), the effects of gravity are not included as they are often negligible considering the strong rotating flows in engines. Additionally, a source and sink of mass will be included which form a simple model for both mass impingement due to water carried by the incoming airstream and droplet pinchoff, respectively. The thin-film equation for a fluid coating an arbitrary substrate has previously been derived by e.g. Myers et al. (Reference Myers, Charpin and Chapman2002), Howell (Reference Howell2003) and Roy et al. (Reference Roy, Roberts and Simpson2002). This formulation was subsequently solved analytically (in an axisymmetric case without shear and gravity) for large substrate curvature by Howell (Reference Howell2003), where matched asymptotic expansions were used to resolve the finite-time shock formation of the leading-order outer solution. When shear and surface tension are included in our formulation, an instability will be observed to emerge between these two competing effects. We can define a line of symmetry parallel to the external shear flow, such that both symmetric and asymmetric solutions are obtained to the steady formulation. It will be necessary for us to include the effect of mass impingement, in a similar manner to that used by Noakes et al. (Reference Noakes, King and Riley2011), in the steady formulation to obtain steady solutions at the same parameter values used to discover the instability in time-evolution simulations. We will demonstrate that, through a formal linear stability analysis, the steady symmetric solutions are unstable and the asymmetric solutions are stable. The instability then emerges from a symmetric profile, and counter-intuitively results in the movement of fluid perpendicular to the direction of shear. This instability will be shown to exist only an interval of shear strengths, and for sufficiently high substrate curvature variation. An analogous symmetry-breaking instability was discovered by Thiele et al. (Reference Thiele, Brusch, Bestehorn and Bär2003) for compactly supported thin-film droplets on a substrate with periodic wettability.
We begin in § 2 by detailing the mathematical formulation for our thin-film problem surrounding an ellipse. In § 3, analytical and numerical methods are detailed that are used to study steady solutions (§ 3.1), linear stability of these (§ 3.2) and the time-evolution problem (§ 3.3). The solution space of steady solutions, obtained with numerical continuation, and the linear stability of these with respect to symmetric and asymmetric temporal perturbations are studied in § 4. In § 5, we perform time-evolution simulations with different values of the shear strength to demonstrate how a symmetry-breaking instability emerges. Lastly, we conclude and discuss our results in § 6.
2. Mathematical formulation
In this work, we consider the flow of a thin viscous fluid surrounding an elliptical substrate, with the inclusion of the competing effects of shear, surface tension, substrate curvature and mass impingement. We first introduce a general curvilinear coordinate system
$(s_1,s_2)$
which parameterises the surface of a three-dimensional substrate,
$\boldsymbol{R} (s_1(x,y,z),s_2(x,y,z) )$
, and then show how this formulation simplifies when the substrate has an elliptical cross-section. Following this general introduction, the remainder of this work will be two-dimensional.
In order to parameterise a coordinate in three-dimensional space, we introduce the normal vector to the surface by
$\boldsymbol{n}(s_1,s_2)$
. Then, a three-dimensional coordinate is defined close to the substrate in Cartesian coordinates
$\boldsymbol{r}=(x,y,z)$
by
\begin{equation} \boldsymbol{r}=\boldsymbol{R}(s_1,s_2)+ \epsilon \eta \boldsymbol{n}(s_1,s_2), \end{equation}
provided that
$\epsilon \ll 1$
. Here,
$\epsilon = L_h/L$
is the ratio between the tangential and normal length scales from the substrate, and
$\eta$
parameterises the normal direction from the substrate for which the fluid surface is given by
$\eta=h(s_1,s_2,t)$
. The surface of the thin viscous fluid surrounding the three-dimensional substrate is defined by evaluating (2.1) at
$\eta =h(s_1,s_2,t)$
. Similarly, locations within the fluid are obtained by taking
$0 \lt \eta \lt h(s_1, s_2, t)$
, and the substrate is specified by
$\eta =0$
.
The demand that the fluid is thin, with
$\epsilon \ll 1$
, both allows for the parameterisation (2.1) to exist, and for the Navier–Stokes equations governing this fluid flow to be simplified under the lubrication approximation. As shown by Myers et al. (Reference Myers, Charpin and Chapman2002), the thin-film equations for the above formulations satisfy (see their equation (31))
\begin{equation} \begin{aligned} \frac {\partial h}{\partial t} + \frac {1}{E^{1/2}}\frac {\partial }{\partial s_1} \bigg [ A_1\frac {h^2}{2}+ \bigg (B \boldsymbol{g} \boldsymbol{\cdot }\boldsymbol{e_1} -\frac {1}{E^{1/2}}\frac {\partial p}{\partial s_1}\bigg )\frac {h^3}{3} \bigg ]&\\ +\frac {1}{G^{1/2}}\frac {\partial }{\partial s_2} \bigg [ A_2\frac {h^2}{2}+ \bigg (B \boldsymbol{g} \boldsymbol{\cdot }\boldsymbol{e_2} -\frac {1}{G^{1/2}}\frac {\partial p}{\partial s_2}\bigg )\frac {h^3}{3} \bigg ]&=f. \end{aligned} \end{equation}
The solution,
$h(s_1,s_2,t)$
, of the nonlinear partial differential equation (PDE) (2.2) represents the distance of the thin film above the substrate (in the normal direction), as seen by evaluating relation (2.1) at
$\eta =h(s_1,s_2,t)$
. The scalar quantities
$E$
and
$G$
, defined by
$E= \partial \boldsymbol{R} /\partial s_1 \boldsymbol{\cdot }\partial \boldsymbol{R} / \partial s_1$
and
$G=\partial \boldsymbol{R} /\partial s_2 \boldsymbol{\cdot }\partial \boldsymbol{R} / \partial s_2$
, are coefficients of the first fundamental form of the substrate,
$\boldsymbol{R}$
. The functions
$A_1(s_1,s_2)$
and
$A_2(s_1,s_2)$
represent the interfacial shear stress in the
$s_1$
and
$s_2$
directions, respectively. The unit vectors,
$\boldsymbol{e}_1$
and
$\boldsymbol{e}_2$
are given by
$\boldsymbol{e}_1= (\partial \boldsymbol{R} / \partial s_1)/E^{1/2}$
and
$\boldsymbol{e}_2= (\partial \boldsymbol{R} / \partial s_2)/G^{1/2}$
, and gravity enters the formulation through the unit vector
$\boldsymbol{g}=(0,0,-1)$
. Lastly,
$f(s_1,s_2,t)$
is a source function, and
$p(s_1,s_2,t)$
is the surface pressure distribution, both of which we will specify later in (2.6c
) and (2.6d
).
All quantities in (2.1) and (2.2) are non-dimensional. We have non-dimensionalised by introducing
$L$
and
$U$
as length and velocity scales tangential to the substrate, and
$L_h$
and
$\epsilon U$
as length and velocity scales normal to the substrate. The fluid pressure has been non-dimensionalised with the lubrication scale
$\mu U L / L_h^2$
, and time with respect to
$L/U$
. The key non-dimensional constants are the aspect ratio,
$\epsilon$
, and the Bond number,
$B$
, which are defined by
\begin{equation} \epsilon = \frac {L_h}{L} \quad \text{and} \quad B= \epsilon ^2 \frac {g L^2}{\nu U}. \end{equation}
Here,
$g$
is the gravitational constant, and
$\nu$
is the kinematic viscosity of the fluid. Note that the gravitation term is included in (2.2) for completeness, but it is not considered in this study (see § 2.2 for more details). Additional non-dimensional parameters will appear in our specifications of the source,
$f$
, pressure,
$p$
, and shear function,
$A_1$
.
2.1. The elliptical substrate
$(a)$
The cross-section of the elliptical substrate,
$\boldsymbol{R}(\theta )$
is shown in the (
$x,y$
)-plane. This is parameterised by
$0 \leqslant \theta \leqslant 2 \pi$
, with the point
$\theta =0$
corresponding to the ‘front’ of the ellipse, and
$\theta =\pi$
to the ‘rear’. The location of the fluid surface is then measured from this substrate with the normal vector
$\boldsymbol{n}$
, as in formula (2.1). A source and sink are located around
$\theta =0$
and
$\theta =\pi$
, respectively, and panel
$(b)$
shows the form of the corresponding mass impingement function
$f(\theta )$
from (2.6c
).
Although the general equation (2.2) is three-dimensional, we shall consider in this work only the two-dimensional dynamics.
We specify the substrate to be a cylinder with elliptical cross-section given in terms of the curvilinear coordinates
$s_1=\theta$
and
$s_2=z$
by the parametrisation
\begin{equation} \boldsymbol{R}(\theta ,z)=({-}\cos {(\theta )}, \delta \sin {(\theta )},z). \end{equation}
Due to the two-dimensional nature of the problem, dependency on
$s_2=z$
is dropped in the resultant analysis. The non-dimensional quantity
$\delta$
, denoted as the ellipse axis ratio, is the ratio between the minor and major axis lengths of the ellipse,
$b$
and
$a$
, hence
\begin{equation} \delta = \frac {b}{a}. \end{equation}
The outward normal vector to this substrate is
$\boldsymbol{n}=({-}\delta \cos {(\theta )},\sin {(\theta )},0)/E^{1/2}$
, where
$E=\sin ^2 {(\theta )}+\delta ^2 \cos ^2{(\theta )}$
. When
$\delta =1$
, the substrate cross-section in the
$(x,y)$
plane is circular. When
$\delta \gt 1$
, the ellipse extends further into the
$y$
direction, and vice versa for
$\delta \lt 1$
. The location
$\theta =0$
corresponds to the ‘front’ of the ellipse, from which increasing the value of
$\theta$
parameterises the substrate in a clockwise manner. This is demonstrated in figure 1
$(a)$
.
Note that the rotation of the ellipse by an angle of
$\pi /2$
corresponds to the change
$\delta \mapsto 1/\delta$
. However, since our non-dimensionalisation fixes the horizontal axis to have unit length, as seen in figure 1, the perimeter changes under this scaling and so to recover the same physical problem we would require also that the total fluid mass
$\mathcal{M}$
undergoes the transformation
$\mathcal{M} \mapsto \mathcal{M}/\delta$
.
2.2. The thin-film equation for an ellipse
Since the elliptical substrate (2.4) remains unchanged in the
$s_2=z$
direction, we make the assumption that the film height
$h(\theta ,t)$
also has no dependency on
$z$
. Then, the governing thin-film equation (2.2) simplifies to
\begin{equation} E^{1/2}\frac {\partial h}{\partial t} + \frac {\partial }{\partial \theta } \bigg ( \frac {h^2}{2}A(\theta ) -\frac {h^3}{3 E^{1/2}}\frac {\partial p}{\partial \theta } \bigg )= f(\theta ), \end{equation}
where
$E=\sin ^2{(\theta )}+\delta ^2 \cos ^2{(\theta )}$
. The shear, source, and pressure functions are defined by
\begin{align} A(\theta )&=S \frac {\sin {(\theta )}}{E^{1/2}}, \end{align}
\begin{align} f(\theta )&= \lambda \Big ( \mathrm{e}^{-4 \theta ^2} + \mathrm{e}^{-4(\theta -2 \pi )^2} - \mathrm{e}^{-4 (\theta -\pi )^2}\Big ), \end{align}
\begin{align} p(\theta ,t) &= \frac {C\delta }{E^{3/2}} - \epsilon C \bigg ( \frac {1}{E}\frac {\partial ^2 h}{\partial \theta ^2}+\frac {\sin {(2\theta )}(\delta ^2-1)}{2E^2}\frac {\partial h}{\partial \theta }+\frac {\delta ^2 h}{E^3}\bigg ). \end{align}
Each of the three non-dimensional constants
$S$
,
$\lambda$
and
$C$
, in ((2.6b
)–(2.6d
)) represent the strength of the effects of shear, mass impingement and surface tension, respectively, and are defined by
\begin{equation} S=\epsilon \frac {\tau _1 L}{ \rho \nu U}, \qquad C= \epsilon ^2\frac {\sigma }{\rho \nu U}, \qquad \lambda = \epsilon ^{-1}\frac {{\dot {m}}}{\rho U}. \end{equation}
Above,
$\rho$
is the fluid density,
$\nu$
is the kinematic viscosity,
$\tau _1$
is the shear stress coefficient,
$\sigma$
is the coefficient of surface tension and
$\dot {m}$
is the source strength.
The shear function,
$A(\theta )$
, in (2.6b
) has been obtained by assuming that the shear is induced by an external flow parallel to the
$x$
-axis, and thus is given by
$A(\theta )=S\, \boldsymbol{e}_{x} \boldsymbol{\cdot }\boldsymbol{e}_1$
.
The source term
$f(\theta )$
, shown in figure 1
$(b)$
, consists of Gaussian functions which result in mass gain near the front of the ellipse,
$\theta =0$
, and mass loss near the rear of the ellipse,
$\theta =\pi$
. A similar source term involving Heaviside functions was used by Noakes et al. (Reference Noakes, King and Riley2011) to study the thin-film flow on the interior of a rotating cylinder. The thin-film reduction of the free-surface pressure, as given in (2.6d
), is more subtle (Howell Reference Howell2003). In general,
$p=-C(x_{\theta }y_{\theta \theta }-y_{\theta }x_{\theta \theta })/{(x_{\theta }^2+y_{\theta }^2)^{3/2}}$
, where
$x(\theta )=-\cos {(\theta )}-\epsilon \delta h(\theta ,t) \cos {(\theta )}/E^{1/2}$
and
$y(\theta )=\delta \sin {(\theta )}+\epsilon h{(\theta ,t)} \sin {(\theta )}/E^{1/2}$
from (2.1). In (2.6d
), following Howell (Reference Howell2003) and taking the thin-film limit, we have retained the
$O(\epsilon )$
terms in the expansion of pressure which contain the second-order partial derivative of the surface
$h$
. This term is retained to regularise finite-time blowup that occurs when only the leading-order component of curvature is used. It is also retained so as to capture the dynamics for both circular and elliptical substrates: notice that in the limit of a circular cylinder,
$\delta \to 1$
and
$E \to 1$
, the leading-order pressure is constant, and therefore the dynamics of (2.6a
) is driven by the
$O(\epsilon )$
terms of (2.6d
).
The total mass of the system will be an important condition to consider for both steady and time-dependent simulations. In Cartesian coordinates, this is given by
$\mathcal{M}=\int \int _{\varOmega } \mathrm{d}x\mathrm{d} y$
, where
$\varOmega$
is the domain of the thin film. By transforming this integral from
$(x,y)$
to
$(\theta ,\eta )$
, where
$0 \leqslant \theta \lt 2 \pi$
and
$0 \leqslant \eta \leqslant h(\theta ,t)$
, with the Jacobian determinant
$\epsilon E^{1/2}+\epsilon ^2 \delta \eta / E$
, we obtain
$\mathcal{M}=\epsilon \int _{0}^{2 \pi } [h E^{1/2} + \epsilon \delta h^2 /(2E) ]\mathrm{d}\theta$
. In considering the leading-order component of this as
$\epsilon \to 0$
, and writing
$\mathcal{M} = \epsilon M$
, we obtain an expression for the mass of the system
\begin{equation} M=\int _{0}^{2 \pi } h(\theta ,t) E(\theta )^{1/2} \mathrm{d}\theta . \end{equation}
For time-dependent simulations, this quantity is used to track conservation of mass. For steady solutions, we fix the value of
$M$
as an additional constraint on the system. Note that, for a fixed value of
$M$
, there exists an exact solution to the thin-film equation (2.6a
) with
$\delta =1$
(for which
$E(\theta )=1$
),
$S=0$
and
$\lambda =0$
, which is given by
$h(\theta ,t)=M/(2\pi )$
. This exact solution will be used as a starting point when obtaining steady solutions to the thin-film equation for other values of
$(\delta ,S,\lambda )$
by numerical continuation in § 3.1.
2.3. Summary of each physical effect
Our thin-film equation (2.6a
) contains the effects of shear (2.6b
), mass impingement (2.6c
) and surface tension (2.6d
). In addition to the three non-dimensional parameters
$S$
,
$\lambda$
and
$C$
characterising the strength of each of these effects, we are also interested in the effect of changing the ellipse axis ratio
$\delta$
. The qualitative effects of including each of these parameters in a time-evolution simulation is displayed in figure 2.
The different time-dependent behaviours that emerge in solutions of the thin-film equation (2.6a
) with the initial condition
$h(\theta ,0)=1$
are shown at
$t=1$
. These solutions were evolved using the time-evolution method described in § 3.3 with
$\epsilon =0.01$
. Arrows depict the movement of fluid. In
$(a\textrm{i})$
, we show a solution profile with
$C=\lambda =0$
and non-zero shear,
$S\gt 0$
. The effect of shear is to transport fluid from the front of the ellipse to the rear. In
$(b\textrm{i})$
, a solution is shown with
$C=S=0$
, and non-zero mass source/sink,
$\lambda \gt 0$
. Here, fluid is gained at the front of the ellipse, and lost at the rear. In
$(c\textrm{i})$
, we show a solution with non-zero surface tension,
$C=1$
, and varying substrate curvature
$\delta \lt 1$
. The effect of surface tension is to drive fluid towards regions of the ellipse with low curvature. Note that in these plots we have enhanced the film height by a factor of
$10$
for visibility. The function
$h(\theta )$
for each of these solutions is shown in panels
$(a\textrm{ii}$
–
$c\textrm{ii})$
.
-
i. The shear strength,
$S$
. This parameter appears in our idealised model for the shear, defined by
$A(\theta )=S \boldsymbol{e}_{x} \boldsymbol{\cdot }\boldsymbol{e}_1$
, from (2.6b). It induces a flux of the fluid towards the rear of the ellipse, an example of which is shown in figure 2
$(a)$
. -
ii. The source strength,
$\lambda$
. This appears in the definition of the mass impingement function
$f(\theta )$
from (2.6c), which contains a source of fluid near the front of the ellipse, and a sink near the rear of the ellipse. The result is a build up of fluid near the front of the ellipse, and net movement of fluid induced towards the rear. An example solution profile containing this effect is shown in figure 2
$(b)$
. For the stationary solutions considered in § 3.1, solutions with
$S\gt 0$
are only attainable for
$\lambda \gt 0$
. -
iii. The capillary number,
$C$
. This appears in the definition of the pressure (2.6d). Without loss of generality, we specify either
$C=0$
, corresponding to zero surface tension, or
$C=1$
, corresponding to non-zero surface tension. This is because if there exists a solution
$h(\theta ,t)$
of (2.6a) with the parameter values
$(S,\lambda ,C)$
, where
$C\gt 0$
, then
$h(\theta ,Ct)$
is a solution of (2.6a) with parameter values
$(S/C, \lambda /C,1)$
. Essentially, by restricting positive
$C$
to take the value of unity, the time scale and definition of the non-dimensional parameters are changed. When
$C=1$
, the effect of surface tension is to minimise the curvature of the fluid surface, and thus fluid moves towards low-curvature regions of the elliptical substrate. This is depicted for
$\delta \lt 1$
in figure 2
$(c)$
. -
iv. The ellipse axis ratio,
$\delta$
. This parameter contributes to each of the above three effects, since the function
$E(\theta )=\sin ^2(\theta )+\delta ^2 \cos ^2(\theta )$
, which appears in ((2.6b
)–(2.6d
)), involves
$\delta$
. Further,
$\delta$
also appears directly in the definition of
$p(\theta )$
from (2.6d
). We consider the values of
$0.5 \leqslant \delta \leqslant 1$
in order to study the effect of substrate curvature on this system.
3. Methodology
We now detail the analytical and numerical methods used to obtain steady solutions, analyse the linear stability of these and perform general time-evolution simulations.
3.1. The steady boundary-value problem
Seeking steady solutions
$h=h(\theta )$
, we set
$\partial h/\partial t=0$
in (2.6a
) and integrate once to obtain the third-order ordinary differential equation (ODE)
\begin{align} &\frac {h^2}{2}A(\theta ) -\frac {h^3}{3 E^{1/2}}\frac { {d\!}{}}{\operatorname {d\!}{}\theta }\left [\frac {C\delta }{E^{3/2}} - \epsilon C \bigg ( \frac {1}{E} \frac { {d}^2 h}{\operatorname {d\!}{}\theta ^2}+\frac {\sin {(2\theta )}(\delta ^2-1)}{2E^2}\frac {\operatorname {d\!}{}h}{\operatorname {d\!}{}\theta }+\frac {\delta ^2 h}{E^3}\bigg )\right ] \nonumber\\ &\quad = \varLambda + \int _{0}^{\theta }f(\theta ') \mathrm{d}\theta ', \end{align}
once the pressure from (2.6d
) is substituted. Since
$A(0)=0$
,
$E(0)=\delta ^2$
and
$E'(0)=0$
, the constant of integration
$\varLambda$
satisfies
\begin{equation} \varLambda = -\frac {\epsilon C}{\delta ^2} (h'''(0)+h'(0)). \end{equation}
We consider periodic solutions that satisfy
$h(\theta )=h(\theta +2 \pi )$
.
Examining (3.1), we see that there are six associated constants:
$\{C, \epsilon , S, \delta , \lambda , \varLambda \}$
. Throughout the paper, we take
$C=1$
without loss of generality, and fix
$\epsilon =0.01$
. In our computations, there is the potential to enforce two further constraints
\begin{align} \text{[$\mathscr{M}$, fixed mass]} & \qquad & \int _{0}^{2 \pi } [h(\theta )-1] E(\theta )^{1/2 }\mathrm{d}\theta &=0, \end{align}
\begin{align} \text{[$\mathscr{A}$, fixed amplitude]} & \qquad & h(\pi ) &= H_1, \end{align}
with
$H_1$
specified and the mass expression follows from (2.8). Note that the mass constraint
$\mathscr{M}$
forces the mean fluid level to be equal to one. The total mass consequently depends on
$\delta$
(since
$E(\theta ) = \sin ^2 \theta + \delta ^2 \cos ^2\theta$
), and this value is chosen to facilitate comparison with the time-dependent simulations from § 5 which use the initial condition
$h(\theta ,0)=1$
.
For instance, specifying both constraints
$\mathscr{M}$
and
$\mathscr{A}$
above, we can then regard the constants
$\lambda$
and
$\varLambda$
of the original six-parameter set as being unknowns of the system. Choices of parameters to set/solve and constraints to impose are clarified below and in Appendix A; different choices are required in order to explore each particular ‘slice’ of the high-dimensional bifurcation space.
In § 4 we study branches of steady solutions in the
$(H_1,\lambda )$
plane by numerical continuation. These include solutions symmetric about
$\theta =0$
and
$\theta =\pi$
, as well as asymmetric solutions. Since symmetric solutions satisfy
$h'(0)=0$
and
$h'''(0)=0$
, they will have
$\varLambda =0$
from (3.2). We only consider solutions for which the film height remains positive across the entire domain, which eliminates solutions such as localised droplets that can form in experiments; this latter situation of dewetting is discussed further at the end of § 3.3.
3.1.1. Spectral solver
The thin-film equation (3.1) is solved with a pseudospectral numerical method that uses Newton iteration.
Let us first discuss the algorithm when
$\mathscr{M}$
and
$\mathscr{A}$
are imposed via (3.3). As noted, we select
$C = 1$
and
$\epsilon = 0.01$
(typically), and then fix values of
$S$
and
$\delta$
. A value of
$H_1$
is also fixed as part of constraint
$\mathscr{A}$
. The unknowns consist of values of
$h$
across the domain,
$\lambda$
, and
$\varLambda$
.
Specifically, we introduce
$N$
collocation points across the domain
$0 \leqslant \theta \lt 2 \pi$
, defined by
$\theta _j=2 \pi j /N$
for
$j=0,1,2 \ldots , N-1$
. Then, an initial guess is specified for the
$N+2$
constants
$(h(\theta _0),\ldots ,h(\theta _{N-1}),\lambda ,\varLambda )$
. Derivatives in the thin-film equation (3.1) are computed spectrally with the relation
$h^{(n)}(\theta )=\mathscr{F}^{-1}[(\mathrm{i} \boldsymbol{\cdot }\text{sgn} (k))^n\mathscr{F}\,[h(\theta )]]$
. Here,
$\mathscr{F}$
denotes the Fourier transform which is evaluated using the fast Fourier transform algorithm. By evaluating (3.1) at each collocation point, and enforcing
$\mathscr{M}$
and
$\mathscr{A}$
, we obtain
$N+2$
conditions to be satisfied. These are minimised by Newton iteration, in which the Jacobian is constructed numerically, such that the residual (the
$L^2$
-norm of the
$N+2$
conditions) lies below
$10^{-8}$
.
The spectral method detailed above produces one solution for the specified values of
$\{\delta , S, H_1\}$
. We now describe how additional solutions are obtained by numerical continuation.
For example, for solutions with
$S\gt 0$
(which consequently have
$\lambda \gt 0$
) we implement continuation within the
$(H_1,\lambda )$
-plane. Here, given one known solution we attempt to converge onto another nearby solution. Our algorithm steps so that the distance between neighbouring solutions is inversely proportional to the local curvature along the solution branch (evaluated from three prior data points). This allows solution branches with high-curvature variations to be computed. The initial guess taken is either a previously computed numerical solution, or the exact solution
$h(\theta )=1$
,
$\lambda =0$
,
$\varLambda =0$
, obtained with
$\delta =1$
,
$S=0$
,
$H_1=1$
and
$M - 2 \pi = 0$
(cf. (3.3a
)). Bifurcations are detected along solution branches by a change in sign of the determinant of the Jacobian obtained in Newton iteration, and these bifurcations are examined in detail in § 4.
We present more specific details of the continuation process, as applied to generate many of the results in this work, in Appendix A.
3.2. Linear stability of steady solutions
By analysing the temporal stability of small-amplitude perturbations to steady solutions, we will be able to explain symmetry-breaking behaviour that emerges in time-evolution simulations of § 5.
We first consider a small time-dependent perturbation,
$h_1(\theta ,t)$
, about a steady solution,
$h_0(\theta )$
, by writing
\begin{equation} h(\theta ,t) \sim h_0(\theta )+ \varepsilon h_1(\theta ,t), \end{equation}
where
$\varepsilon \ll 1$
. In substituting this expansion into the thin-film equation (2.6a
), we obtain at
$O(\varepsilon ^0)$
the steady equation satisfied by
$h_0(\theta )$
. Next, at
$O(\varepsilon )$
we obtain a linear fourth-order ODE for
$h_1(\theta ,t)$
. This is given by
\begin{equation} \frac {\partial h_1}{\partial t}= \alpha _4 \frac {\partial ^4 h_1}{\partial \theta ^4} + \alpha _3 \frac {\partial ^3 h_1}{\partial \theta ^3} + \alpha _2 \frac {\partial ^2 h_1}{\partial \theta ^2} + \alpha _1 \frac {\partial h_1}{\partial \theta } +\alpha _0 h_1, \end{equation}
where the coefficients
$\alpha _i(\theta )$
are known in terms of the steady solution by the relations
\begin{align} \alpha _0= & \frac {Ah_0'-A^{\prime }h_0}{E^{1/2}}+\bigg [\frac {2 h_0 h_0' }{E}+\frac {h_0^2 E' }{2E^2}\bigg ]p_0^{\prime }-\frac {h_0^2 p_0''}{E} \nonumber \\ &+\frac {\epsilon C \delta ^2 h_0^2 }{E^6}\bigg [ E^{\prime \prime }h_0 -\frac {9 {E^{\prime }}^2h_0}{2} -3 E E^{\prime } h_0^{\prime }\bigg ], \end{align}
\begin{align} \alpha _1 = & \frac {Ah_0}{E^{1/2}}+\frac {h_0^2 p_0'}{E}+\frac {\epsilon C h_0^2}{6 E^5}\bigg [(\delta ^2-1)\cos {(2 \theta )}\Big (9EE^{\prime }h_0+6 E^2 h_0^{\prime }\Big )+ 13 \delta ^2 h_0 E^{\prime }\nonumber\\ &+6\delta ^2 E h_0^{\prime }+(\delta ^2-1)\sin {(2 \theta )}\Big ( 2 E E^{\prime \prime }h_0-7{E^{\prime }}^2h_0+4E^2h_0-6EE^{\prime }h_0^{\prime }\Big )\bigg ], \end{align}
\begin{align} \alpha _2&=-\frac {\epsilon C h_0^2}{6 E^4} \bigg [2h_0 E E^{\prime \prime }+\frac {9(\delta ^2-1)(3 E' h_0+ 2 E h_0^{\prime })\sin {(2 \theta )}}{6}\nonumber\\ &\hphantom {=}\,-4 h_0 E (\delta ^2-1)\cos {(2 \theta )}-5{E^{\prime }}^2 h_0 -6 E E^{\prime }h_0^{\prime }-2 \delta ^2 h_0\bigg ], \end{align}
\begin{equation} \alpha _3= -\epsilon C \frac {h_0^2}{6 E^3}\bigg [ h_0 \sin {(2 \theta )}(\delta ^2-1)-5 h_0 E^{\prime }-6 E h_0^{\prime } \bigg ], \qquad \alpha _4=-\epsilon C \frac {h_0^3}{3 E^2}. \end{equation}
Here,
$A(\theta )$
is the shear specified in (2.6b
), and
$p_0(\theta )$
is defined as the
$O(\varepsilon ^{0})$
component of the pressure, obtained by evaluating (2.6d
) at
$h=h_0(\theta )$
. All parameter values
$S, \delta , C, \epsilon , \lambda$
are fixed in the linear stability analysis. Note that while mass (3.3a
) and amplitude (3.3b
) conditions were enforced when calculating steady solutions in § 3.1, these are not enforced on perturbations in the linear stability analysis.
Next, we assume that the perturbation takes the form
\begin{equation} h_1(\theta ,t)= \mathrm{e}^{\sigma t} \bar {h}(\theta ) +c.c., \end{equation}
where
$c.c.$
denotes complex conjugate. In (3.7), the constant
$\sigma$
and function
$\bar {h}(\theta )$
are in general complex valued;
$\text{Re}[\sigma ]$
corresponds to the growth/decay of this mode in time, and
$\text{Im}[\sigma ]$
corresponds to oscillations of this mode in time. Since the PDE (3.5) is linear, it will be possible to solve for
$\bar {h}$
in isolation without regard for the complex-conjugate component. Substitution of ansatz (3.7) into the PDE (3.5) then yields a fourth-order ODE for
$\bar {h}(\theta )$
which contains the unknown constant
$\sigma$
. This ODE is given by
\begin{equation} \sigma h_1= \alpha _4 \frac { {d\!}{} ^4 h_1}{\operatorname {d\!}{} \theta ^4} + \alpha _3 \frac { {d\!}{}^3 h_1}{\operatorname {d\!}{} \theta ^3} + \alpha _2 \frac { {d\!}{}^2 h_1}{\operatorname {d\!}{} \theta ^2} + \alpha _1 \frac {\operatorname {d\!}{} h_1}{\operatorname {d\!}{} \theta } +\alpha _0 h_1, \end{equation}
which we will solve numerically by applying a finite-difference discretisation. In discretising the domain by
$\theta _j=2 \pi j / N$
for
$j=0,\ldots ,N-1$
, we define
$\bar {h}_j=\bar {h}(\theta _j)$
and thus obtain
$N+1$
unknown constants:
$N$
from
$\bar {h}_j$
, and one from
$\sigma$
.
The solution vector is defined by
$\boldsymbol{x}=[\bar {h}_0, \ldots , \bar {h}_{N-1}]^{T}$
. We then discretise each of the derivatives in (3.8) with a centred fourth-order finite-difference approximation. Evaluating the discretised ODE (3.8) at each point in the domain,
$\theta _j$
, then yields
$N$
algebraic equations for
$N+1$
unknown constants. We represent these as the linear system
\begin{equation} \boldsymbol{L} \boldsymbol{x} = \sigma \boldsymbol{x}, \end{equation}
where
$\boldsymbol{L}$
is an
$N \times N$
matrix. Since the ODE (3.8) is local,
$\boldsymbol{L}$
is sparse and only contains non-zero entries on the five main diagonals, as well as the uppermost three and lowermost three diagonals (induced by periodicity). In general, for the
$j$
th row, we have entries
\begin{equation} L_{j,k}=\left \{ \begin{aligned} -\frac {\alpha _4(\theta _j)}{6(\Delta \theta )^4} +\frac {\alpha _3(\theta _j)}{8(\Delta \theta )^3} &\quad \text{if $k=j-3$},\\ \frac {2 \alpha _4(\theta _j)}{(\Delta \theta )^4} -\frac {\alpha _3(\theta _j)}{(\Delta \theta )^3}-\frac {\alpha _2(\theta _j)}{12(\Delta \theta )^2}+\frac {\alpha _1(\theta _j)}{12(\Delta \theta )} &\quad \text{if $k=j-2$},\\ -\frac {13\alpha _4(\theta _j)}{2(\Delta \theta )^4} +\frac {13\alpha _3(\theta _j)}{8(\Delta \theta )^3}+\frac {4\alpha _2(\theta _j)}{3(\Delta \theta )^2}-\frac {2\alpha _1(\theta _j)}{3(\Delta \theta )} &\quad \text{if $k=j-1$},\\ \frac {28\alpha _4(\theta _j)}{3(\Delta \theta )^4}-\frac {5\alpha _2(\theta _j)}{2(\Delta \theta )^2}+\alpha _0(\theta _j) &\quad \text{if $k=j$},\\ -\frac {13\alpha _4(\theta _j)}{2(\Delta \theta )^4} -\frac {13\alpha _3(\theta _j)}{8(\Delta \theta )^3}+\frac {4\alpha _2(\theta _j)}{3(\Delta \theta )^2}+\frac {2\alpha _1(\theta _j)}{3(\Delta \theta )} &\quad \text{if $k=j+1$},\\ \frac {2 \alpha _4(\theta _j)}{(\Delta \theta )^4} +\frac {\alpha _3(\theta _j)}{(\Delta \theta )^3}-\frac {\alpha _2(\theta _j)}{12(\Delta \theta )^2}-\frac {\alpha _1(\theta _j)}{12(\Delta \theta )} &\quad \text{if $k=j+2$},\\ -\frac {\alpha _4(\theta _j)}{6(\Delta \theta )^4} -\frac {\alpha _3(\theta _j)}{8(\Delta \theta )^3} &\quad \text{if $k=j+3$}, \end{aligned}\quad \right . \end{equation}
and zero otherwise, where
$\Delta \theta = 2 \pi /N$
. The exception to this is for
$j=1,2,3$
for which when
$k \leqslant 0$
in formula (3.10), we take
$k \mapsto k+N$
, and for
$j=N-2, N-1, N$
for which when
$k \geqslant N+1$
in (3.10) we take
$k \mapsto k-N$
.
We have thus reduced the study of time-dependent perturbations to the eigenvalue problem (3.9), which we solve for in MATLAB with the inbuilt function eig which implements the QZ algorithm. The most important quantity is the growth rate, Re
$[\sigma ]$
. Steady solutions will be unstable if there exists an eigenvalue with Re
$[\sigma ]\gt 0$
, and stable otherwise.
3.3. The time-evolution problem
We now describe our numerical implementation of the time-evolution problem. The simulations shown later in § 5 will all use the initial condition
$h(\theta ,0)=1$
. The numerical scheme used is an implicit–explicit method where the diffusive terms are treated implicitly and the nonlinear terms explicitly (see e.g. Hammoud Reference Hammoud2016 Appendix B). We first write the thin-film equation (2.6a
) as
\begin{equation} E^{1/2 }\frac {\partial h}{\partial t} + \frac {\partial }{\partial \theta }\bigg ( B_3 \frac {\partial ^3 h}{\partial \theta ^3} + B_2 \frac {\partial ^2 h}{\partial \theta ^2} + B_1 \frac {\partial h}{\partial \theta } + B_0 \bigg ) = f, \end{equation}
where the functions
$B_{i}(\theta ,t)$
depend on both
$\theta$
and
$h(\theta ,t)$
. These functions are given by
\begin{equation} \left .\begin{aligned} B_0(\theta ,t)&=\frac {A h^2}{2} + \frac { C \delta E'h^3}{2 E^3}- \frac { \epsilon C \delta ^2 E' h^4}{E^{9/2}},\\ B_1(\theta ,t)&= \frac {\epsilon Ch^3 }{3 E^{7/2}}( (\delta ^2-1) [\cos {(2 \theta )} E-\sin {(2 \theta )}E']+\delta ^2 ),\\ B_2(\theta ,t)&=\frac {\epsilon C h^3}{6 E^{5/2}}( (\delta ^2-1)\sin {(2 \theta )}-2E^{\prime } ), \qquad B_3(\theta ,t)=\frac {\epsilon C h^3}{3 E^{3/2}}.\\ \end{aligned} \right \} \end{equation}
The problem is treated semi implicitly, such that to obtain values for
$h(\theta ,t)$
at the subsequent time step we require the solution of
$N$
algebraic equations (3.14), which is performed through sparse matrix inversion. The functions
$B_i(\theta ,t)$
will be treated explicitly in time, whereas the derivatives of
$h(\theta ,t)$
in (3.11) will be treated implicitly.
We discretise the domain by introducing the points
$(\theta _i,t_n)$
, given by
\begin{equation} \theta _{j}=2 \pi j / N \quad \text{and} \quad t_n = n \Delta t, \end{equation}
for
$j=0,\ldots ,N-1$
and
$n=0,1,2,\ldots$
, and where the spatial step size is
$\Delta \theta = 2 \pi /N$
and the temporal step size is
$\Delta t$
. The solution values at each gridpoint are then defined by
$h_{j}^{n}=h(\theta _{j},t_n)$
. Derivatives are approximated with finite differences: a second-order centred scheme is used for spatial derivatives, and a first-order forward Euler scheme for the time derivative. The first-order time derivative is discretised with a forward scheme applied at
$(\theta _{j},t_n)$
, yielding
$\partial h /\partial t = (h^{n+1}_{j}-h^{n}_{j})/(\Delta t)$
. This approximation has first-order accuracy in time.
For the spatial derivatives, we use centred second-order finite-difference approximations, in which the temporal behaviour is treated implicitly to resolve stiffness associated with nonlinearity and the fourth-order derivative of the problem. These are first applied to
$\partial /\partial \theta [B_3 h_{\theta \theta \theta } + B_2 h_{\theta \theta } + B_1 h_{\theta } + B_0]$
at
$\theta _{j}$
, yielding
$\partial [B_1 h_{\theta }]/\partial \theta = [(B_1^n h_{\theta }^{n+1})_{j+1/2}-(B_1^n h_{ \theta }^{n+1})_{j-1/2}]/(\Delta \theta )$
for instance. Note here the implicit behaviour in time: we have evaluated coefficients at
$t_n$
, and derivatives of the solution at
$t_{n+1}$
. Further, we have evaluated this first-order derivative at the midpoints between spatial grid points, defined by
$\theta _{j+1/2}=(j+1/2)\Delta \theta$
. Next, when evaluating the derivatives of
$h$
at
$\theta _{j+1/2}$
and
$\theta _{j-1/2}$
, we again use a centred second-order scheme to evaluate these at the gridpoints
$\boldsymbol{H}_{j+1/2} = (h_{j-1}, h_{j}, h_{j+1}, h_{j+2})$
for
$\theta _{j+1/2}$
, and
$\boldsymbol{H}_{j-1/2} = (h_{j-2}, h_{j-1}, h_{j}, h_{j+1})$
for
$\theta _{j-1/2}$
. The formulas for these centred finite differences are given by
$\partial h/\partial \theta \rvert _{j+1/2} = (0,-1,1,0)\boldsymbol{\cdot }\boldsymbol{H}_{j+1/2}/\Delta \theta$
,
$\partial ^2 h/\partial \theta ^2 \rvert _{j+1/2} = (1,-1,-1,1)\boldsymbol{\cdot }\boldsymbol{H}_{j+1/2}/[2(\Delta \theta )^2]$
and
$\partial ^3 h/\partial \theta ^3 \rvert _{j+1/2} = ({-}1,3,-3,1)\boldsymbol{\cdot }\boldsymbol{H}_{j+1/2}/(\Delta \theta )^3$
(Fornberg Reference Fornberg1988).
In summary, we obtain an equation evaluated at
$\theta _j$
for the five unknown quantities
$h^{n+1}_{j-2}$
,
$h^{n+1}_{j-1}$
,
$h^{n+1}_{j}$
,
$h^{n+1}_{j+1}$
and
$h_{j+2}^{n+1}$
, given by
\begin{align} & (\Delta t)\bigg [\frac {E^{1/2}(\theta _j)}{\Delta t}-\frac {{B_1}^n_{j-1}+2{B_1}^n_{j}+{B_1}^n_{j+1}}{2(\Delta \theta )}-\frac {{B_2}^n_{j-1}-{B_2}^n_{j+1}}{2(\Delta \theta )^2}+\frac {3{B_3}^n_{j-1}+6{B_3}^n_j+3{B_3}^n_{j+1}}{(\Delta \theta )^3}\bigg ] h_{j}^{n+1} \nonumber\\ &\qquad -(\Delta t)\bigg [\frac {{B_2}^{n}_{j-1}+{B_2}^{n}_{j}}{2(\Delta \theta )^2}-\frac {{B_3}^{n}_{j-1}+{B_3}^{n}_{j}}{(\Delta \theta )^3}\bigg ] h_{j-2}^{n+1}+(\Delta t)\bigg [\frac {{B_2}^{n}_{j}+{B_2}^{n}_{j+1}}{2(\Delta \theta )^2}+\frac {{B_3}_{j}^{n}+{B_3}_{j+1}^{n}}{(\Delta \theta )^3}\bigg ] h_{j+2}^{n+1} \nonumber\\ &\qquad +(\Delta t)\bigg [\frac {{B_1}^n_{j-1}+{B_1}^n_{j}}{2(\Delta \theta )}-\frac {{B_2}^n_{j-1}+{B_2}^n_j+{B_2}^n_{j+1}}{2(\Delta \theta )^2}-\frac {3{B_3}^n_{j-1}+4{B_3}^n_j+{B_3}^n_{j+1}}{(\Delta \theta )^3}\bigg ] h_{j-1}^{n+1} \nonumber\\ &\qquad +(\Delta t)\bigg [\frac {{B_1}^n_{j}+{B_1}^n_{j+1}}{2(\Delta \theta )}-\frac {{B_2}^n_{j-1}+2{B_2}^n_j+{B_2}^n_{j+1}}{2(\Delta \theta )^2}-\frac {{B_3}^n_{j-1}+4{B_3}^n_{j}+3{B_3}^n_{j+1}}{(\Delta \theta )^3}\bigg ] h_{j+1}^{n+1} \nonumber\\ &\quad = f(\theta _j) \Delta t + E^{1/2}(\theta _j) h_{j}^n - \frac {\Delta t}{E^{1/2}(\theta _j) \Delta \theta }({B_0}_{j+1}^n - {B_0}_{j}^n). \end{align}
Evaluating this equation at
$j=0,\ldots ,N-1$
gives
$N$
algebraic equations for
$N$
unknown constants. Note that, when evaluating (3.14) at
$j=0$
,
$j=1$
,
$j=N-2$
and
$j=N$
, we use the periodic boundary conditions to evaluate the unknown constants. For instance, when
$j=0$
, by periodicity of the solution
$h^{n+1}_{-1}$
is equal to
$h^{n+1}_{N-1}$
and
$h^{n+1}_{-2}$
is equal to
$h^{n+1}_{N-2}$
.
We formulate this set of algebraic equations as the linear system
$\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$
, where the matrix
$\boldsymbol{M}$
is of size
$N \times N$
and the vectors
$\boldsymbol{x}$
and
$\boldsymbol{b}$
are both of size
$N \times 1$
. The
$j$
th entry of the solution vector
$\boldsymbol{x}$
is the unknown
$h^{n+1}_{j}$
. The
$j$
th entry of the vector
$\boldsymbol{b}$
is the right-hand side of (3.14), which is a known quantity. The matrix
$\boldsymbol{A}$
is sparse, with most elements being equal to zero. The non-zero entries then correspond to the coefficients on the left-hand side of (3.14). In general,
$\boldsymbol{A}$
has non-zero entries on the five main diagonals and also non-zero entries induced by periodicity in the top right and bottom left corners. In denoting the element
$[\boldsymbol{A}]_{i,j}$
by
$A_{i,j}$
, these additional non-zero entries are
$A_{1,N-1}$
,
$A_{1,N}$
,
$A_{2,N}$
,
$A_{N-1,1}$
,
$A_{N,1}$
and
$A_{N,2}$
. Inverting this linear system by QR factorisation then yields the solution values at the next time step,
$h_j^{n+1}$
, and the residual associated with this inversion is verified to be at machine precision.
While we only consider solutions that are fully wetted with
$h(\theta ,t)\gt 0$
, some numerical simulations will display behaviour where at some critical value of time the film height reaches zero at a certain value of
$\theta$
. This is associated with dewetting, which in experiments for instance would result in a partially wetted solution that coats only a portion of the substrate. There are multiple methods that can resolve this behaviour by adjusting the governing formulation, for instance by introducing a moving contact line (Myers Reference Myers1998), or a disjoining pressure function (Craster & Matar Reference Craster and Matar2009). We implement a precursor layer thickness (Moriarty et al. Reference Moriarty, Schwartz and Tuck1991), in which the minimum fluid thickness is fixed to a small value,
$h_p=0.005$
. In our semi-implicit numerical scheme this involves setting, at each time step, the film height to
$h_p$
at any values of
$\theta$
where
$h(\theta ,t)\lt h_p$
. Note that this results in the addition of mass into the system when the film height falls below the threshold
$h_p=0.005$
, which can be measured by (2.8). For our simulations in §§ 5.2 and 5.3, this threshold value will only be reached for that with
$S=100$
in figure 10
$(d)$
. In this case, there is a very slight increase in the total mass across the simulation, which is not noticeable in the plotted profiles. All other simulations will have
$h(\theta ,t)\gt h_p$
.
4. Steady numerical solutions
We now obtain numerical solutions to the steady boundary-value problem (3.1), and also calculate the linear temporal stability of these solutions. Firstly, in § 4.1, solutions with
$S=0$
(zero shear strength) are considered. We then consider solutions with
$S\gt 0$
in § 4.2. These consequently have non-zero mass injection,
$\lambda$
, and the resultant solution space consists of both symmetric and asymmetric solutions, as well as non-trivial symmetry-breaking instabilities.
4.1. Solutions with
$S=0$
(zero shear strength)
We first consider solutions with
$S=0$
, corresponding to zero shear strength via the function
$A(\theta )$
defined in (2.6b
). We fix the total mass to the value obtained by evaluating (2.8) with a fluid of uniform height
$h=1$
, and specify both
$C=1$
and
$\epsilon =0.01$
. This results in the free parameters
$\delta$
and
$\lambda$
. As detailed in § 3.1, we typically fix the value of
$\delta$
and
$h(\pi )$
for which
$\lambda$
(the source strength from (2.6c
)) and
$\varLambda$
(the constant of integration in (3.1)) are obtained as part of the solution.
The location of fully wetted steady numerical solutions, with
$S=0$
,
$C=1$
and
$\epsilon =0.01$
are shown in the
$(\delta ,h(\pi ))$
plane in panel
$(a)$
. The mass (2.8) is fixed as an additional constraint to take the value
$M=\int _{0}^{2 \pi } E^{1/2}\mathrm{d}\theta$
. These solutions were obtained using the spectral numerical scheme detailed in § 3.1. Numerical continuation was used to calculate the branch bounding the shaded region, along which the minimum value taken by
$h(\theta )$
is equal to
$0.01$
. For solutions within this shaded region, the minimum value of
$h(\theta )$
is larger than
$0.01$
, and solutions are then characterised by the values
$H_1=h(\pi )$
and
$H_2=h(\pi /2)-h(3 \pi /2)$
, the latter of which measures asymmetry. The solutions shown in panel
$(b)$
have
$\delta =1$
,
$H_2=0$
, and varying
$H_1$
. Those shown in panel
$(c)$
have
$\delta =0.99$
,
$H_1=0.5$
and
$H_2=\{0,0.5,1,1.5,2,2.5 \}$
. Three solutions with min
$[h(\theta )]=0.01$
are then shown in panel
$(d)$
.
Figure 3. Long description
Panel A: A line graph shows the relationship between h(π) and δ. The x-axis is labeled δ ranging from 0.98 to 1.00, and the y-axis is labeled h(π) ranging from 0 to 2. The graph includes points labeled b1 to b5 and d1 to d3, with a shaded region indicating specific solution characteristics. Panel B: A line graph shows the relationship between h and θ. The x-axis is labeled θ ranging from 0 to 2π, and the y-axis is labeled h ranging from 0 to 2. Multiple curves labeled b1 to b5 are displayed. Panel C: A line graph shows the relationship between h and θ for different values of H2. The x-axis is labeled θ ranging from 0 to 2π, and the y-axis is labeled h ranging from 0 to 2. Curves are shown for H2 values from 0 to 2.5. Panel D: A line graph shows the relationship between h and θ for different solutions. The x-axis is labeled θ ranging from 0 to 2π, and the y-axis is labeled h ranging from 0 to 2. Curves are labeled d1, d2, and d3.
An exact solution is known in this parameter space with
$\delta =1$
,
$\lambda =0$
and
$\varLambda =0$
, which has the surface profile
$h(\theta )= 1$
. Further solutions are then obtained from
$h(\theta )=1$
by numerical continuation. When
$\delta =1$
, the solutions exhibit translational invariance and so given one solution
$h(\theta )$
, the profile
$h(\theta - c)$
forms another solution for all
$0 \leqslant c \leqslant 2 \pi$
. This invariance is broken for
$\delta \lt 1$
, for which the lines of symmetry of the ellipse (about
$x=0$
and
$y=0$
) means that given one solution,
$h(\theta )$
, there exist three additional solutions. These are given by
$h(2 \pi -\theta )$
,
$h(\pi -\theta )$
and
$h(\theta -\pi )$
, corresponding to the symmetries about the
$x$
axis, the
$y$
axis and both the
$x$
and
$y$
axes, respectively.
The solutions studied in § 4.2 with
$S\gt 0$
will exhibit symmetry about the line
$\theta =\pi$
as well as bifurcations that break this symmetry. The degree of asymmetry can be measured via the quantity
$H_2 = h(\pi /2) - h(3 \pi /2)$
, which is zero for symmetric solutions that satisfy
$h(\theta ) = h(\pi -\theta )$
. In figure 3
$(b)$
we show solutions with
$\delta =1$
and
$H_2=0$
. Solutions with
$H_2\gt 0$
can be constructed from these by invoking translational symmetry, which results in solutions existing for
$\delta =1$
throughout a region of the
$(H_1, H_2)$
plane. As the value of
$\delta$
decreases, solutions smoothly transition from those profiles with
$\delta = 1$
. For instance, those with
$\delta =0.99$
and
$H_1=0.5$
are shown in figure 3
$(c)$
for
$H_2 = \{0,0.5,1,1.5,2,2.5\}$
. Solutions with
$H_2\lt 0$
can then be obtained from these by symmetry about
$x=0$
and
$y=0$
. Note that when the shear strength
$S$
is slowly increased from zero, we only obtain solutions that bifurcate from those with
$S=0$
and
$H_2=0$
.
The solution space, shown in the
$(\delta ,h(\pi ))$
plane in figure 3
$(a)$
, is dense with solutions existing within the shaded region. As the solutions in figure 3
$(a)$
approach the bounding curve of the shaded region, the minimum value of
$h(\theta )$
tends to zero. This curve was computed by numerical continuation along with the additional constraints that
$\text{min}[h(\theta )]=0.01$
and
$H_2=0$
. Due to this additional condition imposed, we compute solution branches in the
$(\delta ,h(\pi ))$
plane by treating one of
$\delta$
or
$h(\pi )$
as an unknown of the formulation. Further details on this continuation procedure are given in Appendix A. For
$\delta \lt 1$
, fully wetted solutions stop existing due to the film height reaching zero at either
$\theta =0$
or
$\theta =\pi$
. For instance, solutions
$d1$
and
$d2$
in figure 3 have a film thickness equal to
$0.01$
at
$\theta =0$
and
$\theta =\pi$
, respectively. Solution
$d3$
, which is that with the smallest value of
$\delta = 0.9798$
, then has both
$h(0)=h(\pi )=0.01$
.
4.2. Solutions with
$S\gt 0$
(non-zero shear strength)
Solutions obtained with
$S=0$
in the previous section were found to exist only for a small range of values near
$\delta =1$
within
$0.9798 \lt \delta \lt 1.0200$
. Recall that
$\delta$
is the elliptical axis ratio, and that the substrate has larger variations in curvature as
$\delta$
decreases below unity. We now calculate numerical solutions with
$S\gt 0$
and
$\delta \lt 1$
in order to demonstrate the instabilities that emerge with larger variations in the substrate curvature. The following discussion mainly relates to results illustrated in figures 4 and 5.
4.2.1. Solutions with
$S=10$
and
$\delta =0.5$
Numerical solutions of the steady problem (3.1) are shown for the parameter values
$S=10$
,
$C=1$
,
$\epsilon =0.01$
and
$\delta =0.5$
. Solutions exist which are characterised by different values of the injection constant
$\lambda$
, and branches of these are shown in the
$(h(\pi ),\lambda )$
plane in panel
$(a)$
. The branch consisting of points
$b1$
–
$b5$
consists of solutions symmetric about
$\theta =\pi$
, which are shown in panel
$(b)$
. The branch consisting of points
$c1$
–
$c5$
consists of asymmetric solutions, which are shown in panel
$(c)$
. The bifurcation was identified by a change in sign of the determinant of the Jacobian matrix associated with the numerical scheme. The two values of
$(h(\pi ),\lambda )$
indicated by crosses at
$(0.3812,0)$
and
$(2.564,0)$
are taken from asymmetric and symmetric time-dependent solutions obtained later in figure 4
$(c1)$
.
We first consider solutions that exist with
$S=10$
and
$\delta =0.5$
. We will also fix
$C=1$
,
$\epsilon =0.01$
and the value of
$h(\pi )$
. Both of
$\lambda$
and
$\varLambda$
are obtained as unknowns of the problem. We were able to converge directly onto a solution with
$h(\pi )=1$
using Newton iteration with an initial condition
$h(\theta )=1$
,
$\lambda =0$
,
$\varLambda =0$
, and this resulted in a symmetric solution with
$\lambda \gt 0$
. By performing numerical continuation in the
$(h(\pi ),\lambda )$
plane, the branch of symmetric solutions shown in figure 4
$(a)$
was obtained. Our numerical continuation scheme terminated at the points labelled
$b1$
and
$b5$
, at which the decay of Fourier coefficients in our spectral representation of the solution became too poor for Newton iteration to converge. At the end of the branch denoted by
$b1$
, the solution shown in figure 4
$(b)$
has a global minimum at
$\theta =\pi$
where the film height is near zero. The other end of the solution branch, labelled
$b5$
, is associated with solutions where both
$\lambda \to 0$
and the film thickness over a portion of the domain is approaching zero.
A bifurcation from the symmetric branch of solutions shown in figure 4
$(a)$
was identified by tracking the determinant of the Jacobian matrix associated with each solution. This determinant changed sign either side of the point labelled
$b3,c5$
. The constant
$\varLambda$
, which by definition in (3.2) is zero for symmetric solutions, was fixed to a small non-zero value to latch onto the new solution branch associated with the bifurcation. This additional branch consists of asymmetric solutions, and as shown in 4
$(c)$
these solutions can display significant asymmetry. For instance, solution
$c1$
has one pronounced maxima across the entire domain. The asymmetric branch terminated numerically close to
$\lambda =0$
at solution
$(c1)$
due to the Fourier coefficient decay becoming too poor to be captured with a fixed number,
$N=1001$
, of Fourier coefficients. This was associated with the film thickness over a portion of the domain approaching zero.
The temporal stability of steady solutions from figure 4 is shown. In panel
$(a)$
, the maximum growth rate,
$\text{Re}[\sigma ]$
is shown for both symmetric perturbations (dashed line) and asymmetric perturbations (solid line) to symmetric and asymmetric solutions from the branches in figure 4
$(a)$
. Solutions on the asymmetric branch are stable with respect to all linear perturbations. Solutions on the symmetric branch are stable to symmetric perturbations, but unstable to asymmetric perturbations. The full eigenvalue spectrum is shown for the two labelled solutions in panels
$(b)$
and
$(c)$
: here, filled and hollow circles represent the eigenvalues for symmetric and asymmetric perturbations, respectively.
Figure 5. Long description
Panel A: A line graph shows the maximum growth rate as a function of h(π). The x-axis ranges from 0.6 to 1.2, and the y-axis ranges from 0 to 1.0. The graph includes two branches: the symmetric branch (dashed line) and the asymmetric branch (solid line). The inset zooms in on the region between 0.70 and 0.72 on the x-axis, showing detailed behavior around points labeled b and c. Panel B: A scatter plot for symmetric perturbations shows the imaginary part of σ (Im(σ)) on the y-axis ranging from -50 to 50 and the real part of σ (Re(σ)) on the x-axis ranging from -200 to 0. Filled circles represent eigenvalues. Panel C: A scatter plot for asymmetric perturbations shows the same axes and ranges as Panel B, with hollow circles representing eigenvalues.
The stability of solutions along both the symmetric and asymmetric branches from figure 4
$(a)$
is shown in figure 5. Note that in this linear stability analysis we fix the parameter values
$\lambda$
and
$\delta$
, and then the control parameter
$h(\pi )$
is implicitly perturbed via the expansion
$h \sim h_0(\theta ) + \epsilon h_1(\theta ,t)$
from (3.4). Here, we show in 5
$(a)$
the maximum growth rate for both symmetric and asymmetric perturbations along each branch. Solutions along the symmetric branch have both asymmetric and symmetric eigenfunctions in the linear stability analysis, whereas those along the asymmetric branch have only asymmetric eigenfunctions. We see that solutions along the symmetric branch are always unstable. When
$h(\pi )\lt 0.72$
, the dominant growth rate is associated with symmetric perturbations, whereas for
$h(\pi )\gt 0.72$
is is associated with asymmetric perturbations. Solutions along the asymmetric branch in figure 5
$(a)$
are seen in the inset to be stable. The full eigenvalue spectrum for a symmetric and an asymmetric solution is shown in 5
$(b,c)$
. Here,
$\text{Re}[\sigma ]$
corresponds to the growth rate of the temporal perturbation, and
$\text{Im}[\sigma ]$
to the temporal frequency. Eigenvalues corresponding to symmetric modes are shown shaded, whereas those corresponding to asymmetric modes are shown hollow.
Note that the symmetric branch in figure 5
$(a)$
is stable with respect to symmetric perturbations for
$h(\pi )\gt 1.1$
. In the time-evolution simulations described later in § 5, we begin with a symmetric initial condition that initially evolves onto one of these symmetric steady states, which is stable with respect to symmetric perturbations. However, because the symmetric steady state is unstable, the time-evolution simulations eventually destabilise and switch to an asymmetric steady state that is stable with respect to asymmetric perturbations.
4.2.2. Solutions with
$S\gt 0$
and
$0.5\lt \delta \lt 1$
The steady solutions shown in figures 4 and 5 contained symmetric solutions that were unstable to asymmetric perturbations, and asymmetric solutions that were fully stable. These results used
$S=10$
and
$\delta =0.5$
. We now demonstrate how the instability changes as both the shear strength
$S$
and elliptical axis ratio
$\delta$
vary.
Firstly, we consider the stability of solutions with
$\delta =0.5$
and varying shear strength,
$S$
. Specifically, we also fix
$\lambda =2$
for which
$h(\pi )$
becomes an unknown. When
$S=10$
, there are two solutions with
$\lambda =2$
, one symmetric and the other asymmetric. The branches associated with either of these are shown in the
$(S,h(\pi ))$
plane in figure 6
$(a)$
, and these are calculated by numerical continuation. Solutions along both of these branches are shown in figure 6
$(b)$
. The asymmetric branch begins at the point
$b1$
with
$S=8.206$
,
$h(\pi )=0.1596$
, develops significant asymmetry as seen in solution
$b5$
, and then ends at point
$b3$
where
$S=26.70, h(\pi )=6.512$
.
The maximum growth rate is shown for symmetric and asymmetric perturbations for both the symmetric branch in figure 6
$(c)$
, and for the asymmetric branch in 6
$(d)$
. At the start of this branch at
$S=9.3$
symmetric perturbations are unstable with
$\text{Re}[\sigma ]=70.87$
and asymmetric perturbations are stable with
$\text{Re}[\sigma ]=-0.008693$
. The symmetric perturbations stabilise where the branch has infinite gradient in the
$(S,h(\pi ))$
plane at
$S=7.689$
, and the asymmetric perturbations destabilise at
$S=8.196$
where the branch has negative gradient in the
$(S,h(\pi ))$
plane. Between the turning point at
$S=7.689$
and the location
$S=26.68$
, solutions on this branch are unstable to asymmetric perturbations. Within the same range, solutions along the asymmetric branch are stable, as shown in figure 6
$(d)$
.
The behaviour and stability of steady solutions is shown as the shear strength,
$S$
, varies. The fixed parameters are
$C=1$
,
$\epsilon =0.01$
,
$\delta =0.5$
, and
$\lambda =2$
. In panel
$(a)$
, branches of symmetric and asymmetric solutions are shown in the
$S,h(\pi )$
plane. Three solutions along each of the branches are shown in
$(b)$
. The largest growth rate,
$\text{Re}[\sigma ]$
, corresponding to symmetric perturbations (dashed line) and asymmetric perturbations (solid line) is shown in
$(c)$
and
$(d)$
for each of the two solution branches. For the range
$7.689 \lt S \lt 26.68$
, solutions on the symmetric branch are stable with respect to symmetric perturbations, and unstable for asymmetric perturbations. Solutions on the asymmetric branch are stable within the same range of
$S$
.
Stability results for steady symmetric solutions across the
$0.5 \leqslant \delta \leqslant 1$
and
$0 \leqslant S \leqslant 30$
parameter space are shown for
$\lambda =2$
. The value of the largest growth rate,
$\text{Re}[\sigma ]$
, is shown for the linear stability of each solution. Negative growth rates are displayed in dark blue, and positive with colours between light blue and red. The lower black curve corresponds to the lower bound of shear values,
$S$
, at which symmetric solutions exist. Solutions are unstable between the two black curves, where there exists an eigenvalue with
$\text{Re}[\sigma ]\gt 0$
associated with an asymmetric eigenvector, and the largest value of
$\delta$
at which the instability exists is given by
$\delta =0.8150$
. The stability of solutions for
$\delta =0.5$
was previously shown in figure 6
$(c)$
.
This instability associated with asymmetric perturbations to symmetric solutions, shown for
$\delta =0.5$
in figure 6
$(c)$
, is extended to the range
$0.5\lt \delta \lt 1$
in figure 7. To calculate the stability of symmetric solutions within this range, we repeated the calculation shown for
$\delta =0.5$
in 6
$(c)$
to
$100$
values of
$\delta$
within
$0.5\lt \delta \lt 1$
, each spaced apart by
$0.005$
. The maximum growth rate of any perturbation is then reported along the resultant parameter values in figure 7. We see that the symmetry-breaking instability, associated with the range
$8.186\lt S\lt 26.65$
at
$\delta =0.5$
extends to a wider range of values of
$\delta$
, and exists up to the value
$\delta =0.8150$
. Above this value, for
$\delta \gt 0.8150$
, all symmetric solutions are stable.
5. Numerical solutions of the time-evolution problem
The simulations reported in this section use a time step of
$\Delta t = 0.001$
and
$N=400$
spatial gridpoints, corresponding to
$\Delta \theta = 2 \pi / 400$
, unless specified otherwise. The solutions are evolved until a certain time is reached, which is typically
$t = 10$
. Simulations on smaller mesh sizes, up to
$N = 6000$
and
$\Delta t = 1 \times 10^{-5}$
are used when required to verify mesh independence (cf. § 5.2). Comparisons with the steady-state solutions for large
$t\gg 1$
further serve to verify numerical convergence.
5.1. The method of characteristics from Howell (Reference Howell2003)
We begin by studying an exact time-dependent solution that exists with
$S=0$
,
$\lambda =0$
, and when only the leading order of the pressure function
$p(\theta ,t)$
from (2.6d
) is retained. The resultant PDE contains at most first-order derivatives in space and time and is integrable, however, shocks will be seen to appear at finite time. When the
$O(\epsilon )$
terms in the pressure (2.6d
) are retained in the thin-film equation (2.6a
), the PDE contains spatial derivatives up to fourth order that regularise the formation of these shocks. The analysis follows the methodology of Howell (Reference Howell2003), who derived analytical solutions over a general substrate for purely curvature-driven thin-film flows. In setting
$S=0$
, and
$\lambda =0$
, we obtain the first-order nonlinear PDE
\begin{equation} \frac {\partial h}{\partial t} +\frac {1}{E^{1/2}}\frac {\partial }{\partial \theta }\bigg (q \frac {h^3}{3} \bigg )=0, \end{equation}
where we have denoted
$q(\theta )=3C\delta E'(\theta )/(2 E^{3})$
. This quantity is closely related to the curvature,
$\kappa =-C\delta /E^{3/2}$
, of the elliptical substrate by the relation
$q(\theta )=\kappa '(\theta )/E^{1/2}$
.
We will solve (5.1) with the method of characteristics. The idea behind this method is to introduce a coordinate transformation, from
$(\theta ,t)$
to
$(s,\tau )$
, such that the PDE (5.1) will involve differentiation with respect to
$\tau$
only. We first transform the solution by writing
$p(\theta ,t)=q^{1/3}h(\theta ,t)$
, which results in the PDE
$\partial _{t} p+ (q^{1/3}p^2/E^{1/2}) \partial _{\theta }p=0$
. Next, we write the independent variables as functions of
$s$
and
$\tau$
by
$\theta (s,\tau )$
and
$t(s,\tau )$
, and define the solution
\begin{equation} g(s,\tau ) \equiv p(\theta (s,\tau ),t(s, \tau )). \end{equation}
Then, we may use the chain rule to obtain the relation
$\partial g / \partial \tau = (\partial t / \partial \tau )(\partial p/\partial t)+(\partial \theta / \partial \tau )(\partial p / \partial \theta )$
. With an appropriate specification of
$\partial t / \partial \tau$
and
$\partial \theta / \partial \tau$
the left-hand side of PDE (5.1) can be reduced to only involve partial differentiation with respect to
$\tau$
. We thereby obtain three coupled first-order PDEs
\begin{equation} \frac {\partial t}{\partial \tau }=1 ,\qquad \frac {\partial \theta }{\partial \tau }=\frac {q^{1/3}}{E^{1/2}} g^2,\qquad \frac {\partial g}{\partial \tau }=0. \end{equation}
The boundary conditions considered at
$\tau =0$
are given by
\begin{equation} t(s,0)= 0, \qquad \theta (s,0)= s, \qquad g(s,0)= g_0(s). \end{equation}
Note that the original initial condition of
$h(\theta ,0)=h_0(\theta )$
implies that
$p(\theta ,0)=q^{1/3}h_0(\theta )$
, and so
$g_0(s)$
introduced in (5.4c
) above is given by
$g_0(s)=q(s)^{1/3}h_0(s)$
.
Equations (5.3a ) and (5.3c ) may be solved directly along with the boundary condition (5.4a ) and (5.4c ) to find
\begin{equation} t(s,\tau )=\tau \quad \text{and} \quad h(\theta ,t)=\frac {q(s)^{1/3}h_0(s)}{q(\theta )^{1/3}}. \end{equation}
To obtain the solution (5.5b
), we first solved for
$g(s,\tau )=g_0(s)=q(s)^{1/3}h_0(s)$
and then used the relation
$h(\theta ,t)=q^{1/3} p(\theta ,t)$
. To evaluate the right-hand side of (5.5b
), we require the inverse mapping
$s(\theta ,t)$
. For this, we first require a solution for
$\theta (s,\tau )$
, which can be obtained by solving PDE (5.3b
), yielding
\begin{equation} \int _{s}^{\theta }\frac {E(\vartheta )^{1/2}}{q(\vartheta )^{1/3}} \mathrm{d}\vartheta =q(s)^{2/3}h_0(s)^2 \tau . \end{equation}
Equation (5.5c
) specifies the dependence of
$\theta (s,\tau )$
on
$s$
and
$\tau$
. While this expression cannot be directly inverted, it may be solved numerically. This is performed later in § 5.2 in which we compare this exact solution with numerical time-evolution simulations in order to verify our implementation in the absence of shear (
$S=0$
), mass injection (
$\lambda =0$
), and higher-order capillary terms.
A shock will form in the above exact solution when the characteristic projections, which emanate from the initial condition at
$t=0$
, intersect with one another. It is possible to find an analytical expression for the time at which this occurs. To obtain this, we first note that by the chain rule,
$\partial _{s}g= \partial _{\theta }g \partial _{s}\theta$
, since
$\partial _{s}t=0$
by (5.3a
). Then, as
$\partial _{s}g(s,\tau ) = \partial _s{g_0}(s)$
is bounded, a shock forms in the solution when
$\partial _{s}\theta =0$
(corresponding to the blow up of
$\partial _{\theta }g$
). Thus, by taking the partial derivative of expression (5.5c
) with respect to
$s$
, and setting
$\partial _{s} \theta =0$
, we obtain the expression
\begin{equation} t_{s}(s)= -\frac {3E(s)^{1/2}}{6q(s)h_0(s)h_0^{\prime }(s)+2 h_0(s)^2q^{\prime }(s)}. \end{equation}
The minimum time,
$t_c$
, at which a shock forms is then obtained from the conditions
$t^{\prime }_s(s_c)=0$
and
$t_s(s_c)=t_c$
. The functions
$E$
,
$q$
and
$h_0$
are known from the substrate geometry and initial conditions, and so we may directly calculate
$t_c$
. The values for this, and their dependency on
$\delta$
, are shown in figure 8
$(a)$
. The physical location,
$\theta _c$
, at which the shock first forms is shown in figure 8
$(b)$
. Note that, due to symmetry about
$\theta =0$
and
$\theta =\pi$
, there are in general four values in the domain at which this occurs:
$\theta _c$
,
$-\theta _c$
,
$\pi -\theta _c$
and
$\theta _c-\pi$
. The value
$0 \leqslant \theta _c \leqslant \pi /2$
was obtained by first solving for
$s_c$
and
$t_c$
, and then evaluating the integral in (5.5c
) to find
$\theta _c$
.
The results in figure 8 show the time and location at which characteristics of the exact solution first intersect. This occurs in two different ways, corresponding to whether
$q(s_c)$
takes zero or non-zero values. These two different types of shock formation were studied in detail by Howell (Reference Howell2003). The first, a typical shock formation in which the spatial derivative of the solution
$h(x,t)$
becomes unbounded requires
$q(s_c) \neq 0$
, and corresponds to the approximate values
$0\lt \delta \lt 0.87$
and
$1.16 \lt \delta \lt \infty$
. The second, which requires
$q(s_c)=0$
(corresponding to a local minimum or maximum of the substrate curvature) results in the formation of a localised puddle of fluid, and occurs in our example when
$0.87 \leqslant \delta \leqslant 1.16$
. This second case requires more effort to solve analytically due to the turning point in the equation when
$q(s_c)=0$
, and requires the consideration of an inner solution in the work of Howell (Reference Howell2003). The boundaries between these two regimes are shown by dashed lines in figure 8.
The dependence of the exact solution to (5.1) with respect to the ellipse axis ratio
$\delta$
is shown. In
$(a)$
, the critical time for shock formation,
$t_c$
, is shown against
$\delta$
. These values for
$t_c$
are obtained from the minimum points of solution (5.6) for the general shock time across the domain. In
$(b)$
, we show the physical location,
$\theta _c$
, at which the shock emerges. The dashed lines in each plot separate the different types of shock formation. Outside the calculated interval (approximately
$0.87 \leqslant \delta \leqslant 1.16$
), the initial shock that forms results in the spatial derivative,
$\partial h / \partial \theta$
, becoming unbounded at
$\theta _c$
. Within this interval,
$q(s_c)=0$
, and so the initial shock results in a ‘puddle’ which has a bounded spatial derivative.
In the remainder of this work, we consider
$\epsilon \gt 0$
in expression (2.6d
) for the surface pressure, which inhibits the formation of shocks.
5.2. Comparison with the exact solution
The exact solution to the thin-film PDE (2.6a
) with
$S=\lambda =0$
,
$C=1$
, and higher-order terms neglected in the pressure function was detailed in § 5.1 using the method of characteristics. We now validate our numerical time-evolution solution by comparing with this exact solution.
Numerical simulations on the exterior of an ellipse with
$\delta = 0.5$
are shown with the initial condition
$h(\theta ,0)=1$
. Simulation
$(a)$
has
$S=\lambda =0$
,
$C=1$
and neglects
$O(\epsilon )$
terms in the pressure from (2.6d
), so an exact solution exists from the method of characteristics, and simulation
$(b)$
has
$S=\lambda =0$
,
$C=1$
and
$\epsilon =0.01$
. The surface profile for simulation
$(a)$
are shown at
$t=0.04$
in
$(a\textrm{i})$
and
$(a\textrm{ii})$
, and those for simulation
$(b)$
at
$t=6$
in
$(b\textrm{i})$
and
$(b\textrm{ii})$
. The first of these shows this solution with respect to
$\theta$
, and the second projects onto the ellipse in
$(x,y)$
space by formula (2.1). In
$(a\textrm{iii})$
and
$(b\textrm{iii})$
, the numerical solutions
$h(\theta ,t)$
are shown in a continuous heat map across the domain
$0\leqslant \theta \leqslant \pi$
and
$0 \leqslant t$
. The analytical solution for simulation
$(a)$
predicts a shock at
$t_c=0.04303$
and
$\theta _c=0.8069$
, and the characteristics of the exact solution, obtained from solving (5.5c
) for
$\theta (s,\tau )$
, are overlaid.
Figure 9. Long description
Panel A: (ai) A line graph shows the fluid surface profile with the vertical axis labeled h and the horizontal axis labeled θ ranging from 0 to 2π. (aii) A scatter plot projects the surface profile onto the ellipse in space using formula (2.1). (aiii) A heat map displays the numerical solution across the domain θ and t, with a continuous gradient indicating the surface profile. Panel B: (bi) A line graph shows the fluid surface profile with the vertical axis labeled h and the horizontal axis labeled θ ranging from 0 to 2π. (bii) A scatter plot projects the surface profile onto the ellipse in space using formula (2.1). (biii) A heat map displays the numerical solution across the domain θ and t, with a continuous gradient indicating the surface profile. The analytical solution predicts a shock at specific θ and t values, with characteristics of the exact solution overlaid.
The exact solution for
$h(\theta ,t)$
was given in (5.5b
). Note that this expression involves the characteristic variable
$s$
, which can be related to
$\theta$
and
$t$
by numerically evaluating the integral on the left-hand side of (5.5c
). Further, shocks form in the analytical solution when the characteristic projections self-intersect. The time,
$t_s(s)$
, at which this occurs was given in (5.6) as a function of the domain,
$s$
. The earliest time,
$t_c$
, at which a shock forms is therefore a minimum of
$t_s(s)$
, which by definition occurs at
$s_c$
. We thus have the conditions
$t_s'(s_c)=0$
and
$t_s(s_c)=t_c$
, from which the physical location of shock formation,
$\theta _c$
, can be obtained by evaluating (5.5c
). The dependency of
$t_c$
and
$\theta _c$
on
$\delta$
was previously shown in figure 8.
We now consider a numerical simulation with
$S=\lambda =0$
,
$C=1$
,
$\delta =0.5$
and
$h(\theta ,0)=1$
. The exact solution with these parameter values has
$t_c=0.04303$
and
$\theta _c=0.8069$
(to four significant digits). The results of this simulations are shown in figure 9(
$a\textrm{i}$
–
$a\textrm{iii}$
). This simulation used
$N=6000$
gridpoints in order to resolve the solution close to the formation of a shock (typically,
$N=400$
gridpoints would otherwise be appropriate). The numerical solution is shown at
$t=0.04$
in figures 9
$(a\textrm{i})$
and 9
$(a\textrm{ii})$
, where the high surface curvature associated with shock formation is observed. In figure 9
$(a\textrm{iii})$
, a heatmap of the numerical solution is shown for
$0 \leqslant t \leqslant 0.042$
. The black lines shown superimposed in this image are the characteristic projections of the exact solution, obtained by solving (5.5c
).
The shock formation observed in figure 9 is a consequence of the pressure term (2.6d
) only involving the leading-order component as
$\epsilon \to 0$
. The subdominant components of (2.6d
) (which involve higher-order spatial derivatives of
$h(\theta ,t)$
) suppress this shock formation. This is demonstrated by the numerical simulation shown in figure 9
$(b\textrm{i}$
–
$b\textrm{iii})$
, which has the same parameter values as figure 9, but with
$\epsilon =0.01$
used for lower-order components of the pressure function. This solution is able to be evolved beyond the value of
$t_c=0.04303$
, and the solutions shown in figures 9
$(b\textrm{i})$
and 9
$(b\textrm{ii})$
are from
$t=0.6$
. These contain only two maxima across the domain
$0 \leqslant \theta \lt 2 \pi$
. One peak is located at
$\theta =\pi /2$
, and the other at
$\theta =3 \pi /2$
, which correspond to minima of the substrate curvature. The heatmap in figure 9
$(b\textrm{iii})$
demonstrates that initially four maxima develop, similar to the simulation in figure 9, before transitioning to two maxima at later times.
The simulations in this section have been calculated in the absence of shear (
$S=0$
), for which there is symmetry in the solution about
$\theta =\pi /2$
and
$\theta =\pi$
(provided that the initial condition is also symmetric here). In the next section, we calculate time-dependent results with non-zero shear, and study the symmetry-breaking instability that emerges at certain values of
$S$
.
5.3. Competition between shear and surface tension
Four numerical simulations are shown, each with
$\lambda =0$
,
$C=1$
,
$\epsilon =0.01$
and
$\delta =0.5$
. Each grouping of three panels corresponds to shear strength with
$(a)$
$S=0.1$
,
$(b)$
$S=1$
,
$(c)$
$S=10$
and
$(d)$
$S=100$
. The top two insets of each grouping show the solution at
$t = 8$
(solid) and
$t = 2$
(dashed). The bottom heatmap shows the solution across
$0 \leqslant t \leqslant 8$
. The initial condition for each simulation was
$h(\theta ,0)=1$
. Note that the presence of positive shear (
$S\gt 0$
) breaks symmetry about
$\theta =\pi /2$
. Simulations in groups
$(a)$
,
$(b)$
and
$(d)$
result in fluid profiles that are symmetric about
$\theta =\pi$
. In simulation
$(c)$
with
$S=10$
, the profile is initially symmetric about
$\theta =\pi$
but rapidly transitions at
$t=4$
to an asymmetric profile. Note that in panels
$(a\textrm{ii}$
–
$d\textrm{ii})$
, we have amplified the film height by a factor of ten for visibility. The values of
$h(\pi )$
for each of the two solutions shown in panel
$(c\textrm{i})$
were previously shown alongside the steady bifurcation diagram in figure 4.
We now consider time-dependent solutions with non-zero shear,
$S\gt 0$
. The simulations in this section all have
$\lambda =0$
,
$\epsilon =0.01$
,
$C=1$
,
$\delta =0.5$
,
$h(\theta ,0)=1$
and differing values of
$S$
. Solutions with
$S=0$
were previously shown in § 5.2, where surface tension was observed to drive fluid towards areas of low substrate curvature (
$\theta =\pi /2$
and
$3\pi /2$
). The inclusion of a positive shear, which moves fluid from the ‘front’ of the ellipse at
$\theta =0$
to the ‘rear’ at
$\theta =\pi$
, will break front-to-back asymmetry of the solution about
$\theta =\pi /2$
. However, we also crucially demonstrate in this section that for a certain range of shear values the symmetry about
$\theta =\pi$
is also broken due to an instability that emerges.
Time-evolution simulations are shown in figure 10 for four different values of the shear strength,
$S$
. The main figure 10 is divided into four sets of grouped panels corresponding to
$(a)$
$S=0.1$
;
$(b)$
$S=1$
;
$(c)$
$S=10$
; and
$(d)$
$S=100$
. Each simulation begins from the same initial condition at
$t=0$
. For each of the four
$S$
values, we also show solution profiles at
$t = 2$
(dashed) and
$t = 6$
(solid) in the two top insets of each group. Heatmaps show the magnitude of the solution,
$h(\theta ,t)$
.
Those simulations, shown in groupings
$(a)$
and
$(b)$
of figure 10, and corresponding to the smaller values,
$S=0.1$
and
$S=1$
, are similar to those seen previously with
$S=0$
in figure 9
$(b)$
. In these, the solution evolves to a state with two maxima in the domain, located closer to the rear of the ellipse at
$\theta =\pi$
. Indeed, from these we may anticipate that the effect of larger values of
$S$
is to move the majority of the fluid to the rear of the ellipse at
$\theta =\pi$
. This is the behaviour observed in figure 10
$(d)$
, for which
$S=100$
results in a single maxima of the solution at
$\theta =\pi$
. Note that for this simulation with
$S=100$
in figure 10
$(d)$
there is a slight increase in mass over time due to our precursor layer thickness that forces the minimum fluid height to be at or above the value
$h_p=0.005$
.
From the simulation with
$S=10$
in figure 10
$(c)$
, we see that the same behaviour as above is initially observed by the
$t=2$
profile (dashed); this solution is symmetric about
$\theta =\pi$
. However, there is a rapid transition that occurs near to
$t=4$
during which the solution switches to an asymmetric profile. This is shown in the corresponding heatmap. In this example, the resultant solution, shown at
$t=6$
(solid) in the top two insets, has a single peak in the domain at
$\theta \approx 1.269 \pi$
. However, in general, the asymmetric solution can develop the asymmetry either such that the maximum develops at
$\theta \lt \pi$
or
$\theta \gt \pi$
.
The location,
$\theta _{{max}}$
, of the maximum fluid height above the substrate is shown for a range of shear values,
$0 \leqslant S \leqslant 40$
, when
$\delta =0.5$
. The steady solutions from figure 6, with
$\lambda =2$
, are shown with solid lines, and the time-dependent solutions from figure 10 at
$t=8$
, with
$\lambda =0$
, are shown with dotted lines. Each of these solutions have
$C=1$
and
$\epsilon =0.01$
. The circled locations correspond to the specific simulations shown in figure 10 with
$S=0.1$
,
$S=1$
and
$S=10$
. Within a range of shear values,
$8.4 \lesssim S \lesssim 26.8$
, the time-dependent simulations initially tended towards solutions associated with the symmetric branch, before destabilising and tending towards solutions associated with the stable asymmetric branch. Note that by symmetry considerations there is also an asymmetric branch for
$\pi /2 \leqslant \theta _{{max}} \lt \pi$
which we have omitted as only a single maximum will form in each asymmetric simulation.
The results of additional simulations performed with
$ \delta =0.5$
and values
$0 \leqslant S \leqslant 50$
are shown in figure 11. Here, we show the location,
$\theta _{{max}}$
, where the maximum fluid height occurs at
$t = 8$
. We see that for
$0 \leqslant S \lesssim 8.4$
, and
$26.8 \lesssim S \leqslant 50$
, there is either a single peak at
$\theta _{{max}}=\pi$
, or two at
$\theta _{{max}}$
and
$2\pi -\theta _{{max}}$
, since the solution profiles are symmetric. For the values
$8.4 \lesssim S \lesssim 26.8$
, the solution at
$t=8$
is asymmetric and so contains only one maximum within the domain. Note that this instability can develop on either side of
$\theta =\pi$
, which results in an asymmetric profile with a single maximum located either within
$\pi /2 \lt \theta _{{max}} \lt \pi$
or
$\pi \lt \theta _{{max}} \lt 3\pi /2$
. Only the case for
$\pi \lt \theta _{{max}} \lt 3\pi /2$
was shown in figure 11.
This symmetry-breaking instability directly follows from the results of § 4 that studied steady solutions and their linear stability. The results in figure 11 match closely with the range of shear values at which the symmetric solutions with
$\delta =0.5$
were unstable in figure 6
$(c)$
. It was seen in figure 6 that when
$\delta =0.5$
, the branch of symmetric solutions was linearly stable with respect to symmetric temporal perturbations. Further, the same solutions were linearly unstable with respect to asymmetric perturbations, and there also existed branches of asymmetric solutions that were stable with respect to all linear perturbations. The time-evolution simulations of this section have confirmed that the linearly stable solution is attractive with respect to the specified symmetric initial condition
$h(\theta ,0)=1$
. Further, the results show that the asymmetric modes which initially develop from computational error (at the level of machine precision) eventually destabilise the otherwise symmetric solution profile. This result in the time-dependent profile attracting the stable asymmetric solution.
6. Conclusions and discussion
We have studied the two-dimensional thin-film flow of a viscous fluid surrounding a elliptical substrate. The key parameters are the shear strength
$S$
and the ratio
$\delta$
between the semi-minor and semi-major axis of the ellipse. When the substrate curvature varies across the domain, as occurs with
$\delta \lt 1$
, the resultant effect of surface tension is to move fluid from regions of high substrate curvature (
$\theta =0$
,
$\pi$
) towards regions of low substrate curvature (
$\theta =\pi /2$
,
$3\pi /2$
). In contrast, the effect of the shear is to move fluid from the front of the ellipse (
$\theta =0$
) to the rear of the ellipse (
$\theta =\pi$
). These effects are therefore in competition with one another.
We have demonstrated that there can exist a symmetry-breaking instability induced by the inclusion of both these effects. The symmetric solutions obtained in time-dependent simulations are unstable, and display a rapid transition to an asymmetric solution. This was studied in detail for both linear stability of steady solutions in § 4 and time-evolution simulations in § 5. The instability was observed to persist for a wide range of values of the shear strength, and this range decreased as the substrate curvature became more uniform (corresponding to increasing values of
$\delta$
). The instability we have discovered therefore only exists for fluid flows on substrates with sufficiently high curvature variations and certain values of the shear strength.
As noted in § 3.3, we have strictly confined our analysis to the consideration of perfect wetting, and in the time-dependent numerical simulations, the film height is assumed to never decrease lower than a specified precursor layer thickness. In this vein, our approach has followed the motivating work by Myers (Reference Myers1998, Reference Myers2001) and Myers et al. (Reference Myers, Charpin and Chapman2002) and consideration of the original aeroengine application, where understanding the situation of perfectly wetted thin films with shear is important (cf. Peters et al. Reference Peters, Shelton, Tang and Trinh2024 and references therein). An important next step is to consider alternative models for thin films to allow for the situation of partial wetting. One particularly elegant model that does not require explicit treatment of the moving contact line singularity is the use of a disjoining pressure term in the thin-film equation (see e.g. works by Schwartz Reference Schwartz1998; Thiele et al. Reference Thiele, Brusch, Bestehorn and Bär2003; Mayo et al. Reference Mayo, McCue and Moroney2013; and the review in Thiele Reference Thiele2018). The disjoining pressure model would conserve mass in the extreme case (cf. figure 10
$d$
) in which our use of a precursor layer thickness does not. In particular, the work of Thiele (Reference Thiele2011) in understanding the behaviour of dewetting films on rotating cylinders (but without shear and injection) is of particular relevance in better understanding the situation of the dewetting dynamics.
This latter situation of flows on rotating cylinders contains a rich literature, as previously mentioned in our introduction, including Yih (Reference Yih1960), Ruschak & Scriven (Reference Ruschak and Scriven1976), Pukhnachev (Reference Pukhnachev1977) and Peterson et al. (Reference Peterson, Jimack and Kelmanson2001). For example, Duffy & Wilson (Reference Duffy and Wilson1999) study curtain flows on the outside of a rotating circular cylinder, including a sink and source term, but without surface tension or shear. We have not been able to find a direct comparison with the regime studied in our work here (including surface tension, shear and sink/source on differing elliptical geometries). A wide-scale overview of the bifurcation space in all these interconnected works is needed, but even then, we observe in our work that the high-dimensional parameter space is complicated; the addition of rotation, and other relevant effects would no doubt add additional complexity.
While a large number of previous works have considered the flow of thin fluid films surrounding spherical substrates which are driven by the effect of gravity (McKinlay et al. Reference McKinlay, Wray and Wilson2023), comparatively few have considered the effect of shear as the main driving force. For instance, a similar simulation to ours was performed on a sphere (which has constant substrate curvature) by Qin et al. (Reference Qin, Xia and Gao2021), where instead of shear the driving force on the fluid was gravity. This resulted in the movement of fluid from one end of the sphere to the other, and it is possible that a symmetry-breaking instability would emerge in this formulation if the curvature of the substrate were to vary substantially. Such behaviour was not observed in the gravity-driven flow over an ellipsoid by Ledda et al. (Reference Ledda, Pezzulla, Jambon-Puillet, Brun and Gallaire2022), and would likely only emerge for different parameter regimes. The shear model used in our work, e.g. (2.6b ), is certainly an idealised model; in reality, the external flow field will depend on the shape of the viscous fluid that coats the elliptical substrate. There are several avenues for future extensions, including the use of classical potential-flow theory near a Joukowski aerofoil to obtain an estimate of the shear; and also the application of models that use boundary-integral methods to resolve the external flow field combined with the thin-film approximation to resolve the viscous fluid (Meng et al. Reference Meng, Papageorgiou and Vanden-Broeck2024). The source and sink model used in this study is also idealised. In practical, particle tracking and impacting models should be used to evaluate the liquid accretion on the substrate. In addition, a droplet break-off model can be used to replace the sink model.
If a symmetry-breaking instability were to emerge in the formulation with gravity as the driving force, rather than shear, it would also be more robust to reproduce experimentally, which is an exciting prospect.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2026.12160.
Acknowledgements
We thank the three anonymous reviewers for their constructive comments on our work.
Funding
T.P. acknowledges support from the EPSRC Centre for Doctoral Training in Statistical Applied Mathematics at Bath (EPSRC grant no. EP/S022945/1). J.S. and P.H.T. acknowledge support from EP/V012479/1.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Details on the numerical continuation procedure
In § 4 we computed branches on numerical solutions to the steady boundary-value problem formulated in § 3.1. The numerical formulation used, specifically the constraints imposed and parameters treated as unknowns, depend on the specific region of bifurcation space explored. In general, the steady thin-film equation (3.1) contains the six parameters
$(\delta , \epsilon , C, S, \lambda , \varLambda )$
. All of our presented results were for
$\epsilon =0.01$
and
$C=1$
. For the four remaining parameters, we typically treated two of these as unknowns of the problem, subject to requiring the satisfaction of the additional mass constraint
$M=\int _{0}^{2 \pi }h E^{1/2}\mathrm{d}\theta$
(below written
$\mathscr{M}$
), and amplitude constraint
$H_1=h(\pi )$
(below written
$\mathscr{A}$
).
-
(i) (Figure 3) The
$(\delta , h(\pi ))$
plane within figure 3
$(a)$
, we considered
$\{ h(\theta ),\lambda ,\varLambda \}$
as the solution to the system
$\{\text{ODE}\,(3.1), \mathscr{M}, \mathscr{A}\}$
. The boundary of the grey-shaded region is obtained by considering the unknowns
$\{ h(\theta ),\lambda ,\varLambda \}$
as the solution to the system
$\{\text{ODE}\,(3.1), \mathscr{M},\text{min}(h)=\varepsilon \}$
with
$\varepsilon = 0.01$
. -
(ii) (Figure 4) For the solution branches shown in the
$(h(\pi ),\lambda )$
plane in figure 4
$(a)$
, we considered
$\{ h(\theta ),\lambda ,\varLambda \}$
as the solution to the system
$\{\text{ODE}\,(3.1), \mathscr{M}, \mathscr{A}\}$
. The exception to this is when initially latching onto the asymmetric branch (detected via a change in sign of the Jacobian determinant), for which we considered
$\{ h(\theta ),\lambda \}$
as the solution to the system
$\{\text{ODE}\,(3.1), \mathscr{M}\}$
with
$\varLambda$
fixed to a very small non-zero value. -
(iii) (Figure 6) For the solution branches shown in the
$(S,h(\pi ))$
in figure 6
$(a)$
with
$\lambda =2$
, we considered either
$\{h(\theta ),\varLambda \}$
as the solution to the system
$\{\text{ODE}\,(3.1), \mathscr{M}\}$
, or
$\{h(\theta ),\varLambda ,S \}$
as the solution to the system
$\{\text{ODE}\,(3.1), \mathscr{M}, \mathscr{A}\}$
, depending on the local gradient of the bifurcation branch. The former was used when the gradient is close to zero, and the latter used when the gradient is close to infinity.
In all of the cases above, the parameters that we have not mentioned are assumed to take a fixed value.
The Newton iteration scheme used for each of the formulations described above is included in our supplementary material is available at https://doi.org/10.1017/jfm.2026.12160. Our numerical continuation scheme uses the local curvature of the solution branch (calculated from three previous data points) to choose the adaptive step size in parameter space.
(a)
R(θ)
x,y
0⩽θ⩽2π
θ=0
θ=π
n
θ=0
θ=π
(b)
f(θ)
h(θ,0)=1
t=1
ϵ=0.01
(ai)
C=λ=0
S>0
(bi)
C=S=0
λ>0
(ci)
C=1
δ<1
10
h(θ)
(aii
cii)
S=0
C=1
ϵ=0.01
(δ,h(π))
(a)
M=∫02πE1/2dθ
h(θ)
0.01
h(θ)
0.01
H1=h(π)
H2=h(π/2)−h(3π/2)
(b)
δ=1
H2=0
H1
(c)
δ=0.99
H1=0.5
H2={0,0.5,1,1.5,2,2.5}
[h(θ)]=0.01
(d)
S=10
C=1
ϵ=0.01
δ=0.5
λ
(h(π),λ)
(a)
b1
b5
θ=π
(b)
c1
c5
(c)
(h(π),λ)
(0.3812,0)
(2.564,0)
(c1)
(a)
Re[σ]
(a)
(b)
(c)
S
C=1
ϵ=0.01
δ=0.5
λ=2
(a)
S,h(π)
(b)
Re[σ]
(c)
(d)
7.689
S
0.5⩽δ⩽1
0⩽S⩽30
λ=2
Re[σ]
S
Re[σ]>0
δ
δ=0.8150
δ=0.5
(c)
δ
(a)
tc
δ
tc
(b)
θc
0.87⩽δ⩽1.16
∂h/∂θ
θc
q(sc)=0
δ=0.5
h(θ,0)=1
(a)
S=λ=0
C=1
O(ϵ)
(b)
S=λ=0
C=1
ϵ=0.01
(a)
t=0.04
(ai)
(aii)
(b)
t=6
(bi)
(bii)
θ
(x,y)
(aiii)
(biii)
h(θ,t)
0⩽θ⩽π
0⩽t
(a)
tc=0.04303
θc=0.8069
θ(s,τ)
λ=0
C=1
ϵ=0.01
δ=0.5
(a)
S=0.1
(b)
S=1
(c)
S=10
(d)
S=100
t=8
t=2
0⩽t⩽8
h(θ,0)=1
S>0
θ=π/2
(a)
(b)
(d)
θ=π
(c)
S=10
θ=π
t=4
(aii
dii)
h(π)
(ci)
θmax
0⩽S⩽40
δ=0.5
λ=2
t=8
λ=0
C=1
ϵ=0.01
S=0.1
S=1
S=10
8.4≲S≲26.8
π/2⩽θmax<π
Open access
