2.1. The generalized thin film equation

The most frequently used version of the thin film equation (Oron, Davis & Bankoff Reference Oron, Davis and Bankoff1997; Craster & Matar Reference Craster and Matar2009; Blossey Reference Blossey2012) is that for a layer of viscous liquid on a solid substrate. It describes how viscous motion is driven by pressure gradients. In the lubrication approximation, in which the pressure $p(x,t)$ is taken to be constant over the layer, the equation for the film thickness $h(x,t)$ becomes

(2.1)\begin{equation} h_t = \tfrac{1}{3}(h^3 p_x)_x. \end{equation}

Here subscripts refer to differentiation with respect to the variable. For simplicity, we restrict ourselves to one dimension, and generalize to higher dimensions later. If the fluid is allowed to slip over the solid surface, as is the case for entangled polymer solutions (Münch, Wagner & Witelski Reference Münch, Wagner and Witelski2005; Blossey Reference Blossey2012), the mobility $h^3$ changes its power to $h^2$. The mobility $h$ is observed for the break-up of a thin neck inside a Hele-Shaw cell (Constantin et al. Reference Constantin, Dupont, Goldstein, Kadanoff, Shelley and Zhou1993).

If it were for surface tension alone, a flat film would be a state of minimum energy, and the film would relax back to it, even if perturbed (Benzaquen et al. Reference Benzaquen, Fowler, Jubin, Salez, Dalnoki-Veress and Raphael2014). However, the liquid–gas interface is often attracted by the substrate (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009), and rupture can occur. This is modelled by an additional disjoining pressure, and in the case of non-retarded van der Waals forces (Blossey Reference Blossey2012),

(2.2)\begin{equation} p ={-}\gamma h_{xx} - \frac{A}{6{\rm \pi} h^3}, \end{equation}

where $\gamma$ is the surface tension and $A>0$ the Hamaker constant. Other exponents for the van der Waals attraction are also possible, for example, $p \propto h^{-4}$ for retarded van der Waals forces (de Gennes Reference de Gennes1985; Blossey Reference Blossey2012).

Without aiming to model a particular system, but noting that exponents may differ according to experimental circumstances, we now allow the two exponents to vary freely. We also choose units of $h$, $x$ and $t$ such that prefactors become unity, to arrive at the equations

(2.3)\begin{equation} h_t - (h^n p_x)_x = 0, \quad p ={-}h_{xx} + \frac{h^{{-}m}}{m}, \end{equation}

which have been studied numerically for a long time (Williams & Davis Reference Williams and Davis1982; Sharma et al. Reference Sharma, Kishore, Salaniwal and Ruckenstein1995). The classical case of non-retarded van der Waals forces, and a Newtonian liquid on a solid surface, corresponds to $n=m=3$. Alternatively, this ‘generalized lubrication equation’ (Eggers & Fontelos Reference Eggers and Fontelos2009; Dallaston et al. Reference Dallaston, Fontelos, Tseluiko and Kalliadasis2018) may be written in the compact form

(2.4)\begin{equation} h_t + ( h^nh_{xxx} + h^{n-m-1} h_x)_x = 0. \end{equation}

In the case of a two-dimensional film, which is the main subject of this paper, the lubrication equation $h_t + {\boldsymbol {\nabla }}\boldsymbol {\cdot }{\boldsymbol {f}} = 0$ now has the flux ${\boldsymbol {f}} = -h^n{\boldsymbol {\nabla }} p$, and the pressure is $p = -\triangle h + h^{-m}/m$. Here ${\boldsymbol {\nabla }}$ is the gradient operator and $\triangle \equiv \nabla ^2$. Once more we have chosen units to normalize coefficients to unity. Now the generalized lubrication equation (Bertozzi & Pugh Reference Bertozzi and Pugh1994) becomes

(2.5)\begin{equation} h_t ={-}{\boldsymbol{\nabla}}\boldsymbol{\cdot}(h^n {\boldsymbol{\nabla}} \triangle h + h^{n-m-1}{\boldsymbol{\nabla}} h). \end{equation}

Assuming that $h(x,y,t)$ depends on a single spatial variable, one of course recovers (2.4).

To test the stability of a flat film of thickness $h_{ref}$, we consider small perturbations of size $\epsilon$ and linearize

(2.6)\begin{equation} h(x,y,t) = h_{ref} + \epsilon \exp({\omega t + \textrm{i} (k_x x + k_y y)}), \end{equation}

where $\omega$ is the growth rate of the perturbation. Inserting (2.6) into (2.5), we obtain the dispersion relation

(2.7)\begin{equation} \omega = h_{ref}^{n-m-1} k^2(1 - h_{ref}^{m+1} k^2), \quad k^2 = k_x^2+k_y^2. \end{equation}

One observes that wavenumbers $k < h_{ref}^{-(m+1)/2}$ are unstable, while the highest growth rate is achieved for $k = h_{ref}^{-(m+1)/2}/\sqrt {2}$. In our numerical simulations below, we will consider films sufficiently thin, so that the smallest wavenumber in the system is linearly unstable. After an initial exponential growth of perturbations, the problem becomes nonlinear and eventually leads to rupture, where the film thickness goes to zero.

2.2. One-dimensional similarity solutions and their stability

We begin with the simplest case of a single spatial variable, corresponding to rupture of a one-dimensional ridge. We consider (2.4) and look for similarity solutions of the form

(2.8)\begin{equation} h(x,t) = t'^{\alpha} H(\xi),\quad \xi = \frac{x'}{t'^{\beta}}, \end{equation}

where $t' = t_0 - t$ and $x' = x - x_0$. Inserting this into (2.4), and balancing the powers of $t'$ that arise from each of the three terms, we find a unique solution for the exponents in terms of the parameters $n$ and $m$,

(2.9a,b)\begin{equation} \alpha = \frac{1}{2 + 2m - n}, \quad \beta = \frac{1 + m}{4 + 4m - 2n}. \end{equation}

Using these values, we obtain an equation for the similarity profile $H(\xi )$,

(2.10) \begin{equation} -\alpha H + \beta\xi H_{\xi} ={-}[H^n H_{\xi\xi\xi} + H^{n-m-1}H_{\xi}]_{\xi}. \end{equation}

As observed for the case of $n=m=3$ (Zhang & Lister Reference Zhang and Lister1999a), (2.10) has an infinite sequence of solutions, which we will denote by $H_i(\xi )$, $i=0,1,2,\dots$. However, only one of them, the ‘ground state’ solution $H_0(\xi )$, is observed in simulation. This raises the question of the stability of such solutions, which can be studied by rewriting (2.4) in self-similar form. This is achieved by the transformation (Giga & Kohn Reference Giga and Kohn1985)

(2.11)\begin{equation} h(x,t) = t'^{\alpha} H\left(\xi,\tau \right), \quad \xi = x'/t'^{\beta}, \quad \tau ={-}\ln t', \end{equation}

so that the similarity equation becomes

(2.12) \begin{equation} H_{\tau} = \alpha H - \beta\xi H_{\xi} -[H^n H_{\xi\xi\xi} + H^{n-m-1}H_{\xi}]_{\xi}. \end{equation}

The advantage of this ‘dynamical system’ description (Eggers & Fontelos Reference Eggers and Fontelos2015) is that solutions of (2.10) are now fixed points of (2.12), making stability much easier to study. To investigate the neighbourhood of a similarity solution, (2.12) has to be solved subject to the far-field condition (Eggers & Fontelos Reference Eggers and Fontelos2015; Witelski & Bernoff Reference Witelski and Bernoff2000)

(2.13)\begin{equation} H_{\tau} \sim \alpha H - \beta\xi H_{\xi}, \quad |\xi| \rightarrow \infty, \end{equation}

which ensures matching to a slowly evolving far-field solution.

To investigate the linear stability of solutions $H_i(\xi )$ of (2.10), we put $H(\xi ,\tau ) = H_i(\xi ) + \delta e^{\nu \tau } P(\xi )$, and linearize (2.12) in $\delta$. This results in the eigenvalue problem

(2.14) \begin{align} \nu P &= \alpha P - \beta\xi P_{\xi} - [H_i^n P_{\xi\xi\xi} + n H_i^{n-1}H_{i,\xi\xi\xi}P \nonumber\\ &\quad +\, H_i^{n-m-1}P_{\xi} + (n-m-1)H_i^{n-m-2}H_{i,\xi}P]_{\xi} \end{align}

for eigenvalues $\nu = \nu _R + \nu _I i$, which is subject to the far-field condition

(2.15)\begin{equation} \nu P \sim \alpha P - \beta\xi P_{\xi}, \quad |\xi| \rightarrow \infty, \end{equation}

which follows from (2.13). If all eigenvalues have negative real parts, as it is typically found for the ground state $i=0$ (Bernoff, Bertozzi & Witelski Reference Bernoff, Bertozzi and Witelski1998; Eggers Reference Eggers2012), the solution is stable. The higher-order solutions $i>0$ all have at least one eigenvalue with a positive real part, and are thus unstable. This excludes two positive eigenvalues of (2.14) (Giga & Kohn Reference Giga and Kohn1985; Eggers & Fontelos Reference Eggers and Fontelos2015), which are always present owing to the invariance of the equations under space and time translations, while (2.8) assumes a fixed reference point $x_0,t_0$. A perturbation will produce a shift in the spatial and temporal position of the singularity, leading to the solution being driven away from (2.8). The corresponding positive eigenvalues are $\nu _T = 1$ and $\nu _X = \beta$ for the temporal and spatial shift, respectively. However, these eigenvalues do not describe any physical instability, and can thus be discounted.

The stability of (2.10) was investigated by Witelski & Bernoff (Reference Witelski and Bernoff2000) for the case $n=m=3$, ignoring eigenvalues with an imaginary part. We will see below that the most dangerous eigenvalues (those with the largest real part, apart from $\nu _T$ and $\nu _X$) are in fact complex. Consistent with other systems, the ground state similarity solution was found to be stable, all other solutions unstable. This remains to be true if the missing complex solutions are accounted for (Dallaston et al. Reference Dallaston, Fontelos, Tseluiko and Kalliadasis2018). We begin by studying singular solutions of the one-dimensional equation (2.4) only.

2.3. Regular and irregular motion

In Dallaston et al. (Reference Dallaston, Tseluiko, Zheng, Fontelos and Kalliadasis2017), starting from the known solution for $n=m=3$, new solutions were found by numerical continuation, holding the mobility exponent fixed at $n=3$, but lowering $m$ from its typical value $m=3$, making the interaction more long ranged. The stability of this branch of solutions was studied in Dallaston et al. (Reference Dallaston, Fontelos, Tseluiko and Kalliadasis2018). Here we extend the analysis of Dallaston et al. (Reference Dallaston, Tseluiko, Zheng, Fontelos and Kalliadasis2017, Reference Dallaston, Fontelos, Tseluiko and Kalliadasis2018) to the whole interval $1\le n\le 3$, as represented in the phase diagram shown in figure 1. We focus on the stability of the ground state solution $H_0(\xi )$, ignoring higher-order branches $H_i$, $i\ge 1$, which are always unstable. In the upper part of the diagram, labelled regular, symmetric, one finds symmetric (stable) ground state solutions, an example of which is shown in figure 2(a), for $n=3,m=2$, and marked by a star in figure 1. These solutions are similar to those found originally for $n=m=3$.

Figure 1.

Phase diagram of one-dimensional singular solutions of (2.4). Below the symmetric Hopf bifurcation (solid blue line), regular rupture solutions become unstable, and periodic, or even more complex dynamics appear. A second antisymmetric bifurcation occurs just below (red dashed line). For small values of $n$ below $\approx 1.87$, regular, asymmetric solutions are seen below a Hopf bifurcation in the reverse direction (orange dot–dashed line). Examples, seen in figure 2, are marked by stars.

Figure 2.

Three examples of similarity solutions (solutions to (2.12)), identified by stars in figure 1. (a) A regular, symmetric similarity solution for $n=3,m=2$. (b) A regular, asymmetric solution for $n=1.5,m=0.2$. (c) A periodic solution for $n=3,m=1.3$, from the antisymmetric branch. The periodic solution has a period of $T=3.7$, and four profiles are shown, at the times indicated in (d), where we plot the minimum scaled thickness $H_{min}$ as a function of the scaled time $\tau$; note that profile 1 is the mirror image of profile 3.

Of particular interest are the eigenvalues with the two largest real parts (both of which are complex), the larger one of which has a symmetric eigenfunction, the smaller (more negative) eigenvalue corresponding to an antisymmetric eigenfunction. As $m$ is lowered from $3$, the real parts increase, to produce Hopf bifurcations when $\nu _R=0$ at finite $\nu _I\ne 0$; this signals the birth of a periodic orbit (Drazin Reference Drazin1992), while the original fixed point solution $H_0(\xi )$ has become unstable. This happens first for the symmetric eigenfunction at a value $m_s$, and then for the antisymmetric eigenfunction, at $m_a$. We refer to the corresponding bifurcations as the symmetric and antisymmetric bifurcations, respectively, and the solution branches originating from those bifurcations as the symmetric and antisymmetric branches. Extending our numerical procedure to search for periodic solutions of (2.12) (see the Appendix A for details), we can continue the periodic solution branches from the Hopf bifurcation, where they are created. The symmetry of the periodic solutions along the symmetric and antisymmetric branches corresponds to the eigenfunctions which become unstable at the bifurcation. On the symmetric branch, the profiles are symmetric for all times: $H(\xi ,\tau ) = H(-\xi ,\tau )$. On the antisymmetric branch, on the other hand, the second half of the period is the mirror image of the first half: $H(\xi ,\tau ) = H(-\xi ,\tau + T/2)$.

The search for Hopf bifurcations of the ground state solution has been carried out for $1\le n\le 3$, which results in bifurcation curves $m_s(n)$ and $m_a(n)$, which are drawn as the solid blue and dashed red lines in figure 1, respectively. Below the solid blue line, the symmetric fixed point solution has become unstable, and one should expect irregular behaviour: periodic, or perhaps even more complicated. As seen in figure 1, the symmetric Hopf bifurcation line is close to linear, and well described by the equation $m = 0.8(n-1)$. As a good rule of thumb, we can therefore say that

(2.16)\begin{equation} m > 0.8(n-1):\quad \text{regular dynamics} \end{equation}

is the condition for a regular fixed point solution; in the opposite case, irregular behaviour may be seen.

We have not yet fully described the bifurcation structure below the antisymmetric Hopf line, but arguments given in Appendix A.3 indicate that only the antisymmetric periodic solution branch should be observed. This is confirmed by numerical simulations of the one-dimensional equation (2.3), whose periodic solutions always exhibit the symmetry of the antisymmetric branch.

In figure 2(c,d) we show a typical example of such a solution from the antisymmetric branch, which is found by continuing the periodic solution branch from the antisymmetric Hopf bifurcation, as explained in Appendix A. It corresponds to a solution $H(\xi ,\tau )$ in similarity space, which is periodic in $\tau$ with period $T$, after which it returns to the same profile. In addition, the profiles obey the symmetry $H(\xi ,\tau ) = H(-\xi ,\tau + T/2)$. Thus, if profiles are recorded at discrete times $\tau _n = \tau _0 + nT$, where $n$ is an integer, one observes a self-similar shrinking of the solution, hence, the name ‘discrete self-similarity.’ In figure 2(c) we show a sequence of four profiles at equal steps of $\tau$, after one returns to the original profile.

The difference between a fixed point solution and a periodic solution is illustrated in figure 3, where direct numerical simulations of the one-dimensional thin film equation (2.4) are shown. On the left, we show the regular case: the same symmetric profile appears on smaller and smaller scales. On the right, we show the irregular case of a periodic orbit (discrete self-similarity) for the same parameter values as those shown in figure 2(c,d). The profile undergoes a sequence of changes, which repeat themselves after a period $T$ in logarithmic time.

Figure 3.

Simulations of (2.3) with $n=3,m = 2$ (a, regular case), and $n=3,m = 1.3$ (b, irregular case), plotting $\log _{10}h(x,t)$ as the singularity is approached. Only the neighbourhood of the singularity is shown, and a new profile is recorded each time the minimum thickness $h_{min}$ has decreased by a factor of $0.8$.

In Dallaston et al. (Reference Dallaston, Fontelos, Tseluiko and Kalliadasis2018) it is shown with the example of $n=3,m=1$, that solutions can also become non-periodic. While in the discretely self-similar case the profile continually undergoes the same sequence of instabilities, in the non-periodic case new structures are always created. While the periodic case can be seen as creating a monofractal (similar to a Cantor set Halsey et al. Reference Halsey, Jensen, Kadanoff, Procaccia and Shraiman1986), where the same pattern appears on smaller and smaller scales, the non-periodic case corresponds to a multifractal, which consists of a superposition of sets with different scaling exponents. Apart from some specific examples, we have not yet delineated regions of more complicated behaviour in the phase diagram.

However, as seen in figure 1, not all solutions below the blue bifurcation line are irregular. For $n\lesssim 1.87$, below the orange dot–dashed bifurcation line, fixed point solutions become once more stable, but are now asymmetric, an example of which is shown in figure 2(b). As detailed in Appendix A, these solutions were found from a continuation procedure from the original symmetric ground state solution. Beyond the symmetric and antisymmetric Hopf bifurcations, a pitchfork bifurcation leads to a pair of unstable, asymmetric solutions with complex eigenvalue; the asymmetric pair is the mirror image of one another. Continuing this branch, the real part of the eigenvalue eventually goes through zero and becomes negative in another Hopf bifurcation, but going in the opposite direction, and producing a stable, asymmetric solution.

We have confirmed the structure of the one-dimensional phase diagram in figure 1 by performing numerical simulations of (2.3) with a spatial resolution of up to 14 orders of magnitude (two examples of which are shown in figure 3). We have also performed simulations of the rescaled equation (2.12), as explained in more detail in Appendix A. There remains some uncertainty as to the exact structure of the diagram in its lower half, below $m \approx 0.4$. There, the dynamics is characterized by very long transients, which make the distinction between regular and irregular behaviour difficult. We could identify some return to fixed point behaviour for sufficiently small $m$, even to the right of the orange bifurcation line. However, we must leave the exact delineation of these boundaries to a future publication.

2.4. Higher-dimensional solutions

We now come to the main aim of this paper, which is to understand the structure of solutions in two and higher dimensions, using the general ansatz (1.1). We first analyse the case of $\bar {\beta } = \beta$, leading to pointlike solutions. Then we analyse the case $\bar {\beta } < \beta$, which results in quasi-one-dimensional solutions.

2.4.1. Pointlike solutions

Inserting (1.1) with $\bar {\beta } = \beta$ into the generalized thin film equation (2.5), we obtain the similarity solution

(2.17)\begin{equation} -\alpha H + \beta\xi H_{\xi} + \beta\eta H_{\eta} ={-}{\boldsymbol{\nabla}}\boldsymbol{\cdot}(H^n {\boldsymbol{\nabla}} \triangle H + H^{n-m-1}{\boldsymbol{\nabla}} H), \end{equation}

where ${\boldsymbol {\nabla }} = (\partial _{\xi },\partial _{\eta })$ and $\triangle = \partial _{\xi }^2 + \partial _{\eta }^2$. In view of spatial isotropy, it is natural to look for solutions of (2.17) which are axisymmetric (Zhang & Lister Reference Zhang and Lister1999b; Witelski & Bernoff Reference Witelski and Bernoff2000). The scaling of the exponents (2.9a,b) is the same as in the one-dimensional case. This leads to the radially symmetric solution

(2.18)\begin{equation} h(x,y,t) = t'^{\alpha} H(\rho), \quad \rho = r / t'^{\beta}, \end{equation}

for which the similarity equation becomes in two dimensions

(2.19)\begin{equation} -\alpha H + \beta\rho H_{\rho} + \frac{1}{\rho} \left[ \rho H^n \left(\frac{\left(\rho H_{\rho}\right)_{\rho}}{\rho}\right)_ {\rho} + H^{n-m-1} H_{\rho}\right]_{\rho} = 0. \end{equation}

Solutions to (2.19) have been computed for $n=m=3$ (Zhang & Lister Reference Zhang and Lister1999b), and have been found to be stable to radial as well as non-axisymmetric perturbations (Witelski & Bernoff Reference Witelski and Bernoff2000). However, the stability analysis again excluded complex eigenvalues. We will see below that axisymmetric solutions are in fact stable only for a range of $n,m$ values.

It is worth noting that pointwise similarity solutions are not necessarily axisymmetric. An example are optical caustics, which are described by the eikonal equation (Nye Reference Nye1999). In that case, radially symmetric solutions would correspond to perfect focusing, and would be the least generic situation. Instead, the most stable solutions, which take the fewest number of adjustable parameters to find, are the elliptic and hyperbolic umbilic. They can be found by solving the analogue of (2.19) for the eikonal equation (Eggers Reference Eggers2020).

2.4.2. Quasi-one-dimensional solutions

The other possible solution of the form (1.1) is the case $\bar {\beta } < \beta$, which means that the solution is varying slowly in the transversal direction; the direction in which the solution is most singular (largest derivatives), we choose as the $x$-direction. We have investigated solutions of this type previously in the case of the eikonal equation (Eggers et al. Reference Eggers, Hoppe, Hynek and Suramlishvili2015), the dispersionless Kadomtsev–Petviashvili equation (Grava et al. Reference Grava, Eggers and Klein2016) and shock formation in the compressible Euler equation (Eggers et al. Reference Eggers, Grava, Herrada and Pitton2017). Here we point out the generality of the approach, which is independent of the particular structure of the equation, and apply it to the generalized thin film equation.

Owing to the slow variation in the $y$-direction, we can look at the solution as a superposition of one-dimensional solutions (2.8). Since the eigenvalue $\nu _T$ is usually the most singular one, we consider a corresponding shift in the singularity time of the one-dimensional solution. As $y$ is varied, the singular time $t_c(y)$ effectively varies (Pomeau et al. Reference Pomeau, Le Berre, Guyenne and Grilli2008; Grava et al. Reference Grava, Eggers and Klein2016). Setting the origin of the transversal variable such that the singularity occurs for $y=0$ first, and expanding, we have

(2.20)\begin{equation} t_c = t_0 + ay^2 + O(y^3); \end{equation}

we must have $a>0$ so that $t_c \ge t_0$ always. For the same reason, there can be no linear term. Replacing $t'$ in (2.8) by $t_c - t = t' + a y^2 + O(y^3)$ implies the scaling $y^2 \propto t'$ and, thus, $\bar {\beta } = 1/2$. Of course, $\bar {\beta } = 1/4$, corresponding to $t_c - t = t' + a' y^4 + O(y^5)$ is another possibility, but that would be a non-generic situation. The singularity proceeds in the $x$-direction, but is ‘unfolded’ in the transversal ($y$) direction, a situation well understood in the case of the eikonal equation (Nye Reference Nye1999; Eggers et al. Reference Eggers, Hoppe, Hynek and Suramlishvili2015).

Thus, we write for the similarity solution describing the unfolding,

(2.21)\begin{equation} h(x,y,t) = t'^{\alpha} H\left(\frac{x'}{t'^{\beta}},\frac{y'}{t'^{1/2}}\right) \equiv t'^{\alpha} H\left(\xi,\eta\right), \end{equation}

which is consistent if $\beta > 1/2$. In that case the $h_x = t'^{\alpha -\beta }H_{\xi } \gg h_y = t'^{\alpha -1/2}H_{\eta }$, so that the right-hand side of (2.5) becomes to leading order ${ -t'^{\alpha -1}[H^n H_{\xi \xi \xi } + H^{n-m-1}H_{\xi }]_{\xi }}$, and the similarity equation of the thin film problem is

(2.22)\begin{equation} -\alpha H + \beta\xi H_{\xi} + \frac{\eta}{2} H_{\eta} ={-}[H^n H_{\xi\xi\xi} + H^{n-m-1}H_{\xi}]_{\xi}. \end{equation}

Combining the condition $\beta > 1/2$ for a quasi-one-dimensional singularity with the formula (2.9a,b) for $\beta$, we obtain

(2.23)\begin{equation} n > 1 + m, \quad \text{quasi-one-dimensional} \end{equation}

as a necessary condition for a quasi-one-dimensional singularity. The same stability boundary is obtained by considering a small $y$-dependent perturbation to a one-dimensional solution (Witelski & Bernoff Reference Witelski and Bernoff2000), realising that after a shift in the $y$-direction, it is enough to consider quadratic perturbations to leading order. Indeed, for the case $n=m=3$, for which a pointlike solution should be observed according to (2.23), Witelski & Bernoff (Reference Witelski and Bernoff2000) describe numerically a transition of line-like and ring-like solutions toward point rupture.

With the above insight, we use (2.22) directly to find the most general unfolding of a one-dimensional solution $H^{(1)}(\xi )$ of (2.10). We claim that

(2.24)\begin{equation} H(\xi,\eta) = (1 + A\eta^2)^{\alpha} H^{(1)}\left(\frac{\xi + B\eta^{2\beta}}{\left(1 + A\eta^2\right)^{\beta}}\right) \end{equation}

is a solution of (2.22) for arbitrary constants $A,B$; these constants are adjustable parameters, which depend on initial conditions. Indeed, putting $\zeta = (\xi + B\eta ^{2\beta })/ (1 + A\eta ^2)^{\beta }$, the right-hand side (RHS) of (2.22) becomes

(2.25) \begin{equation} \text{RHS} ={-}(1 + A\eta^2)^{\alpha-1} [H^n H_{\zeta\zeta\zeta} + H^{n-m-1}H_{\zeta}]_{\zeta}, \end{equation}

while we compute

(2.26) \begin{align} -\alpha H + \beta\xi H_{\xi} &= (1 + A\eta^2)^{\alpha-1} \left(-\alpha H^{(1)} + \beta\xi H^{(1)}_{\xi} + \frac{\eta}{2} H^{(1)}_{\eta}\right)\nonumber\\ &\quad + (1 + A\eta^2)^{\alpha-1}[-\alpha A \eta^2 H^{(1)} + \beta A \eta^2 \zeta H^{(1)}_{\zeta}]\nonumber\\ &\quad - (1 + A\eta^2)^{\alpha-\beta} \beta B \eta^{2\beta} H^{(1)}_{\zeta}\nonumber\\ &= \text{RHS} - \frac{\eta}{2} H_{\eta}, \end{align}

using the fact that $H^{(1)}$ is a solution of the one-dimensional problem (2.10).

We stress that the unfolding (2.24) will provide solutions for any isotropic problem (not just the thin film equation), as long as the condition $\beta > 1/2$ is satisfied. The generality of (2.24) can also be understood from noting that since $H^{(1)}(\xi )$ does not contain any free parameters, the only way a dependence on $y$ can come in is through $t_0(y)$ and $x_0(y)$. Then to leading order $t_c - t = t' + A y^2 = t'(1 + A\eta ^2)$, and if we replace $t'$ in

(2.27)\begin{equation} h = t'^{\alpha}H^{(1)}\left(\frac{x'}{t'^{\beta}}\right) \end{equation}

by $t'(1 + A\eta ^2)$, we obtain the terms in (2.24) involving $A$. Next, we allow $x_c(y)$ to depend on $y$, and expand in $y$,

(2.28)\begin{equation} x - x_c = x' - \sum_{i=1}^{\infty} b_i y^i, \end{equation}

so that

(2.29)\begin{equation} \frac{x - x_c}{t'^{\beta}} = \xi-\sum_{i=1}^{\infty} b_i \eta^i t'^{i/2-\beta}. \end{equation}

This only leads to a finite result if $2\beta$ is an integer, so that in fact the contribution $B\eta ^{2\beta }$ in (2.24) is regular. All coefficients $b_i$ with $i/2-\beta < 0$ must vanish (otherwise they would blow up), and so the leading contribution involving $B$ is that of (2.24).

The form of the solution (2.24) is a general feature, independent of the structure of the equation. In the case of additional symmetries, it might even have a more general form. For example, considering the cusp singularity (Nye Reference Nye1999) of the eikonal equation, and using the equivalent of (2.24), one finds a solution with three independent parameters. Only after invoking an additional invariance of the eikonal equation, does one recover the required four parameters.

In the case of the dynamical system (2.12), important for describing the irregular behaviour discussed in § 2.3, the equation extending (2.22) is

(2.30)\begin{equation} H_{\tau} = \alpha H - \beta\xi H_{\xi} - \frac{\eta}{2} H_{\eta} -\left[H^n\left(H_{\xi\xi} - \frac{1}{m H^m}\right)_{\xi}\right]_{\xi}; \end{equation}

now if $H^{(1)}(\xi ,\tau )$ is a solution of (2.12),

(2.31)\begin{equation} H(\xi,\eta) = (1 + A\eta^2)^{\alpha} H^{(1)}\left(\frac{\xi + B\eta^{2\beta}}{(1 + A\eta^2)^{\beta}}, \frac{\tau}{1 + A\eta^2} \right) \end{equation}

solves (2.30). Note that depending on $y$, the solution is found in different phases of the evolution. In the case of a non-periodic solution, this will be even more complicated.

3.1. Matlab

In this scheme, constructed with Matlab, we have taken advantage of the symmetry of the problem to simulate just a quarter of the domain: $[0,0.5]\times [0,0.5]$. The two-dimensional generalized thin film equation (2.5) leads to a nonlinear equation for $h$ in terms of $h_t$ and partial derivatives in the spatial variables,

(3.2)\begin{align} h_t &= h^n[h_{yyyy} + 2h_{xxyy} + h_{xxxx} + (h_{yy}+ h_{xx})/h^{m + 1} - (h_y^2+h_x^2)(m + 1)/h^{m + 2} \nonumber\\ &\quad + h_yh^{n - 1}n(h_{yyy} + h_{xxy} + h_y/h^{m + 1}) + h_xh^{n - 1}n(h_{yyx} + h_{xxx} + h_x/h^{m + 1})]. \end{align}

The time derivative is discretised using second-order backward differences while second-order central differences were employed to discretise the spatial derivatives. The method is fully implicit; to solve the nonlinear system resulting from the discretisation, a Newton–Raphson technique was used, where the required Jacobian matrix is obtained by combining analytical functions and collocation matrices (see Herrada & Montanero (Reference Herrada and Montanero2016) for further details). This allows us to take advantage of the sparsity of the resulting matrix, to reduce the computational time to invert it.

A variable time step $\delta t$, based on the change in the minimum of $h$, was used for the time integration. The minima $h_{min}$ and $h_{min,1}$ at the current and previous time steps, respectively, as well as the last time step, $\delta t_1$, are used to set

(3.3)\begin{equation} \delta t_{try}=(\alpha h_{min}-h_{min,1})/(h_{min}-h_{min,1})\delta t_1-\delta t_1. \end{equation}

Here $\alpha$ is a parameter to control the variation of the time step. The step $\delta t_{try}$ is used when the solution approaches the singularity, while a fixed time step $\delta t_{fixed}$ is used at the beginning of the simulation. To that end, the time step was set as $\delta t=\textrm {min}(\delta t_{try},\delta t_{fixed})$. Most of the simulations were carried out with $\alpha =0.9$ and $\delta t_{fixed}=0.001$.

All simulations were started with a uniform grid with a grid size $\Delta x = \Delta y = 0.002$. At a time determined individually for each case, the solution was interpolated to a new grid, concentrated around ($x_{min},y_{min}$), where $h$ reaches its minimum value. Near this point a uniform grid with spacing $\Delta x = \Delta y = \varDelta _{min}$ is used. Starting at this point and after each $N$ grid points, the grid size was doubled in each direction. The simulations were carried out with $\varDelta _{min}=5\times 10^{-5}$ and $N=40$. The black lines in figure 4(a) show an example of this finer mesh, which was used for $h_{min}\le 0.09$.

Figure 4.

A comparison of the two numerical codes for the example of figure 5. (a) Thin black lines show the ‘fine’ grid of the Matlab code, used for $h_{min}\le 0.09$; the thick red lines show one quadrant of the most refined grid of the Basilisk code, for $h_{min} = 2.6 \times 10^{-3}$. The Matlab code reached $h_{min} = 6 \times 10^{-4}$. (b) Plot of $h_{min}(t)$ for both codes.

3.2. Basilisk library

Here the nonlinear system of equations is solved using a Newton–Raphson scheme, implemented using the multigrid solver of the free C-language library Basilisk (Popinet Reference Popinet2015). This allows us to make use of the sophisticated adaptive grid refinement schemes available for Basilisk, and to control the discretisation error directly. Basilisk provides tools for the numerical solution of partial differential equations (PDEs), discretised spatially with adaptive Cartesian meshes.

Setting $p = -h_{xx} - h_{yy}$, (2.5) reads in components as

(3.4) \begin{equation} h_t = (h^n p_x)_x + (h^n p_y)_y - (h^{n-m-1} h_x)_x - (h^{n-m-1} h_y)_y. \end{equation}

Starting from a trial solution $(p^*,h^*)$, we set $p = p^* + \delta p$ and $h = h^* + \delta h$, with $(p,h)$ the converged solution, and $(\delta p,\delta h)$ a small correction. Expanding $h^n = (h^* + \delta h)^n\approx (h^*)^n + n (h^*)^{n-1}(h-h^*)$, we obtain a linear system equivalent to (3.4), up to terms quadratic in $\delta h, \delta p$,

(3.5)\begin{align} \left. \begin{gathered} p + h_{xx} + h_{yy} = 0 , \\ (\alpha p_x)_x + (\alpha p_y)_y + (\beta h_x)_x + (\beta h_y)_y + (\gamma_x h)_x + (\gamma_y h)_y = h_t + (\gamma_x h^*)_x + (\gamma_y h^*)_y , \end{gathered} \right\} \end{align}

with

(3.6)\begin{align} \alpha = (h^*)^n, \: \beta ={-}(h^*)^{n-m-1} \quad \text{and} \quad \gamma_{x,y} = n(h^*)^{n-1} p^*_{x,y} -(n-m-1)(h^*)^{n-m-2} h^*_{x,y}. \end{align}

The system (3.5) is solved for $h,p$ up to convergence, using a multigrid method. Convergence is checked by monitoring $\delta h$ a posteriori.

Spatial derivatives are discretised using second-order differences, while the time derivative has been discretised using first- or second-order backward differences; no significant difference was found between both time discretisation methods. The time step has been set according to two different strategies; it is either set proportional to the minimum height, or by estimating the time truncation error, calculated by comparing the solution obtained with second- and first-order time derivatives. The equations are very stiff, and prone to numerical divergence if the time step is not small enough. Therefore, both strategies require some fine tuning of parameters, in order to approach the singularities as closely as possible. The spatial discretisation is adapted automatically by monitoring the inverse of the height field, $h_i = 1/h$. In the vicinity of the singularities grids as small as $\Delta x = \Delta y = 3 \times 10^{-5}$ are used while the grid size grows to $\Delta x = \Delta y = 7.8 \times 10^{-3}$ far away from the singularities.

The main advantage of the Basilisk code is its automatic adaptivity of the mesh which can handle several singularities simultaneously, which will be illustrated in figure 12 of the discussion. As a drawback, the discrete set of linearized equations (3.5) is difficult to solve with Basilisk's iterative multigrid solver. The multigrid solver is extremely delicate to handle in this simulation and is prone to diverge if the trial solution is not close enough to the true solution (either because the time step is too large and/or the interpolation in the remeshing is not accurate enough). This problem does not arise in the Matlab code, given that the linear system is solved exactly by the direct inversion of a sparse matrix. As a constraint, the Matlab code can only handle a single singularity, since the stretching of the grid cannot deal with multiple regions of refinement.

We compare the two codes in figure 4, which use the same parameters as in figure 5, for which pointlike behaviour is observed. On the right it is shown that the minimum thickness $h_{min}$ as a function of time agrees closely between the two simulations. However, with the Basilisk code we are only able to integrate to $h_{min} = 2.6\times 10^{-3}$, owing to the multigrid solver failing to converge. The adaptive grid that was generated at this point is shown as the thick red lines in figure 4(a). Up to this point the two codes agree, as shown by plotting $h_{min}(t)$ on the right. With the Matlab code, we were able to integrate down to $h_{min} = 6\times 10^{-4}$, where for $h_{min} \le 0.09$ the finer, graded grid was used, shown as black lines on the left of figure 4.

Figure 5.

Simulation of (2.5) with $n=2,m=1.5$, ($\alpha \approx 0.33,\beta \approx 0.42$), and initial condition (3.1), using $\epsilon _1 = 0.05, \epsilon _2 = 0.03$ and $h_{ref} = 0.2$. (a) A perspective plot of $1/h$ at $\tau = 10.1$ demonstrates the pointlike character. Plots (b,c) show cuts in the $x$ and $y$ directions, respectively, for the values of $\tau = -\ln t'$ shown. Profiles are collapsed according to (2.18), and agree with a solution of (2.19) (dot–dashed line), demonstrating axisymmetry.