3.1. Experimental observations

In this section, we briefly present selected results from the study of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) in the presence and absence of the spanwise inlet forcing devices shown in figure 1(a). Figure 2(a) shows a film thickness distribution obtained using an absorption technique (reproduced from Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). The inlet spanwise thickness profile is amplified, and streamwise-saturated and steady rivulets are observed downstream. The saturated rivulets are periodic along the spanwise direction. There is a narrow range of attainable spacings, when the inlet is forced, around the value $\hat {L}_r=2\pi \sqrt {2}l_c^*$ (value shown in figure 2a), i.e. the most amplified wavelength in the dispersion relation of the flat film. Interestingly, even in the absence of the spanwise inlet forcing devices, the predominant spacing of the emerging rivulet structures is $\hat {L}_r$ (see figure 2b).

Figure 2.

(a) Film thickness for $\theta = 39^{\circ }$ and $h_N = 1515\ \mathrm {\mu }\text {m}$ ($u=1.5$), steady inlet forcing with the sinusoidal blade at the wavelength $\hat {L}_f = 2\pi \sqrt {2}l_c/\sqrt {\sin \theta }$, extracted from Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). The thickness is measured with the absorption method and normalized by the flat film thickness $h_N$. The in-plane distances are normalized by the reduced capillary length $l_c/\sqrt {\sin \theta }$. (b) Typical rivulet pattern in the absence of the inlet forcing devices (figure 1a), $\theta =20^\circ$.

However, far downstream in figure 2(a), oscillations appear on the rivulet profiles. These oscillations are amplified and rivulets carrying travelling lenses are observed, for these values of angle and flow rate.

3.2. Numerical observations

The aim of this section is to numerically study the emerging patterns for an initial condition that mimics the experimental conditions described in the previous section.

We consider a gravity-driven thin film of viscous Newtonian fluid flowing under a planar substrate inclined with respect to the vertical with an angle $\theta$. We introduce the following non-dimensionalization:

(3.1a–d)\begin{equation} x=\hat{x}/l_c^*; \quad y=\hat{y}/l_c^*; \quad z=\hat{z}/h_N; \quad t=\hat{t}/\tau^*, \end{equation}

where $\tau ^*={\nu l_c^2}/{h_N^3 \sin ^2(\theta ) g}$ is the characteristic time scale of the Rayleigh–Taylor instability. The numerical model for the evolution of the film thickness $h$ is the lubrication equation in which the complete expression of the curvature is retained (Ruschak Reference Ruschak1978; Wilson Reference Wilson1982; Weinstein & Ruschak Reference Weinstein and Ruschak2004):

(3.2)\begin{equation} \partial_t h +u h^2 \partial_x h +\tfrac{1}{3} \boldsymbol{\nabla} \boldsymbol{\cdot} [ h^3 \boldsymbol{\nabla} h + h^3 \boldsymbol{\nabla} \kappa]=0,\end{equation}

where $\boldsymbol {\nabla }$ operates in the $(x,y)$ directions, $u= \cot {(\theta )} \tilde {l}_c^*$ and $\tilde {l}_c^*= {l}_c^*/h_N$. The linear advection velocity $u$ corresponds to the surface film velocity at which linear interface thickness perturbations with respect to a flat condition are advected downstream (Brun et al. Reference Brun, Damiano, Rieu, Balestra and Gallaire2015). In physical quantities, an increase of the parameter $u$ implies a decrease of the flow rate (since $u$ is inversely proportional to $h_N$) or $\theta$. The curvature $\kappa$ reads

(3.3)\begin{equation} { \kappa}=\frac{\partial_{{x}{x}} {h}(1+(\partial_{{y}}{h})^2) +\partial_{{y}{y}} {h}(1+(\partial_{{x}}{h})^2) - 2\partial_{{x}{y}}{h}\partial_{{x}}{h}\partial_{{y}}{h}}{(1+(\partial_{{x}}{h})^2 +(\partial_{{y}}{h})^2)^{3/2}}. \end{equation}

The two-dimensional equation is implemented in COMSOL Multiphysics. We use the built-in finite elements method solver, exploiting cubic elements with Lagrangian shape functions and a fully implicit time solver. The largest mesh element size is half of the reduced capillary length $\tilde {l}_c^*$. The domain size is $L_x\times L_y$, where $L_x=231$ and $L_y=106$, leading to approximately $50\,000$ elements. A convergence analysis was performed, showing that convergence is achieved for this characteristic size of the elements. This characteristic element size was also validated by the experimental and numerical comparisons in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). The equations are solved for the variables ($h,\kappa$). For all the considered cases, periodic boundary conditions are used.

Experimentally, in the absence of the spanwise inlet perturbation device described in figure 1(a), the rivulet spacing is the one dictated by the most amplified mode in the flat film dispersion relation, i.e. ${L}_r=2\pi \sqrt {2}$ (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). We numerically study the nonlinear time evolution when a streamwise-invariant sinusoidal initial condition is considered. We choose as initial condition a sinus of wavelength ${L}_r=2\pi \sqrt {2}$:

(3.4)\begin{equation} h(x,y,t=0)=\bar{h}_N\left(1+{A}\cos\left(\frac{2\pi y}{L_r}\right)\right), \end{equation}

where ${A}=10^{-2}$, and $\bar {h}_N=0.54$ is the initial value of the thickness that gives, for a pure streamwise saturated structure, the same local flow rate in the streamwise direction as a flat film of thickness $h=1$ (§ 5.3 in Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020).

We introduce the moving reference frame at the linear advection velocity $u$$(\xi =x-u t,y)$. Figure 3 shows the evolution of the thickness with time for (a) $u=5.45$ and (b) $u=1.5$. For visualization purposes, we focus in the region $\xi \in [-8\pi \sqrt {2},8\pi \sqrt {2}]$ and $y \in [-6\pi \sqrt {2},6\pi \sqrt {2}]$. In both cases, the streamwise invariant initial condition is amplified and reaches, at $t=800$, a saturated state in the streamwise direction. For (a) $u=5.45$, we do not observe any further evolution of the pattern for $t>800$. For (b) $u=1.5$, at $t=800$ the rivulet profiles saturate. For $t>800$, however, streamwise thickness perturbations grow, and at $t=1200$ the flow is characterized by lenses travelling on the rivulets.

Figure 3.

Nonlinear response in the case of a streamwise-invariant sinusoidal initial condition, for (a) $u=5.45$ and (b) $u=1.5$. From left to right: $t=1000$, $t=1200$. Results are reported in the moving reference frame at the linear advection velocity $(\xi =x-u t,y)$.

The streamwise-invariant sinusoidal initial condition is amplified leading to a rivulet pattern saturated in space and time, periodic along the spanwise direction. The absence (respectively presence) of observable streamwise perturbations on the rivulet profiles at high (respectively low) values of $u$ suggests that the stability of the rivulet profile to streamwise perturbations may be directly related to the linear advection velocity.

The experimental observations of predominant spanwise-periodic rivulet patterns and the occurrence of lenses on the rivulets are confirmed by the nonlinear simulations with periodic boundary conditions. In the following, we aim at rationalizing the emergence of predominant rivulet structures and their destabilization.

4.1. Base flow

In this section, we define the saturated rivulet profile $H_r(y)$, serving as base flow for the local stability analysis. The numerical base flow is the large-time solution ($t=10\ 000$) of the one-dimensional model presented in § 5 of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). The profile, of periodic wavelength ${L}_r$, is given by a one-dimensional model in which the flow rate in the streamwise direction coincides with the one of a flat film of thickness $h= 1$, leading to a mean value $\bar {h}_N=0.54$ of the thickness of the rivulet. The numerical procedure revealed that the rivulet profile is slowly saturating to a steady state $H_r(y)$. In figure 4, we report the numerical periodic profile at $t=10\ 000$ (solid line) used for the stability analysis. The rivulet is characterized by a central lobe of large thickness that saturates to a steady profile described by the pendulum equation (red circles in figure 4), while the side lobes (of very low thickness) are slowly draining with a power law $t^{-1/2}$ (Lister et al. Reference Lister, Rallison and Rees2010). It is remarkable that, with the considered non-dimensionalization, the profiles are independent of $u$, i.e. there is a unique rivulet shape (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). The numerical profile agrees well with the experimental results (dots in figure 4) and can therefore be safely used as base flow $H_r(y)$ for the stability analysis.

Figure 4.

Periodic rivulet profile (black line) used for the stability analysis, compared with the results of the pendulum equation of § 5.4 of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) (red circles), and with the experimental results for three central rivulets (grey dots), from Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020), for 10 transverse measurements at two different streamwise locations, at $\theta =39^\circ$ and different $h_N$. The red dashed line denotes the mean thickness $\bar {h}_N$ of the rivulet.

4.2. Dispersion relation

Following the classical approach of the local stability analysis, we consider as a base state the single, spanwise-periodic and steady rivulet $H_r(y)$ described in § 4.1. The quasi-steadiness of the rivulet profile allows us to neglect the slow evolution of the side lobes at long times and thus to consider a normal mode expansion in time and along the direction in which the base state is invariant, i.e. the streamwise direction $x$ (Schmid, Henningson & Jankowski Reference Schmid, Henningson and Jankowski2002). The spanwise periodicity governing the base state $H_r(y)$ is also enforced on the perturbation. The following normal mode decomposition is therefore used:

(4.1)\begin{equation} h(x,y,t)=H_r(y)+\varepsilon \tilde{\eta}(x,y,t)=H_r(y)+\varepsilon \eta(y)\, \textrm{e}^{\mathrm{i}(k_x x -\omega t)}, \quad {\varepsilon \ll 1,} \end{equation}

where $\tilde {\eta }$ is the thickness perturbation with respect to the base flow profile $H_r(y)$. By considering the two-dimensional nonlinear equation (3.2) and introducing the normal mode decomposition (4.1), one obtains, up to $ {O}(\varepsilon )$

(4.2)\begin{align} &\varepsilon \partial_t \tilde{\eta} + \varepsilon u H_r^2\partial_x \tilde{\eta} + \frac{1}{3} \partial_y \left[ H_r^3\left(\frac{\mathrm{d}H_r}{\mathrm{d} y} +\frac{\mathrm{d}\kappa_{(0)}}{\mathrm{d} y}\right)+ \varepsilon H_r^3 \partial_y\tilde{\eta} \right.\nonumber\\ &\quad \left. +\,\varepsilon H_r^3 \partial_y \tilde{\kappa}_{(1)} + 3\varepsilon H_r^2\left(\frac{\mathrm{d}\kappa_{(0)}}{\mathrm{d} y} + \frac{\mathrm{d}H_r}{\mathrm{d} y}\right) \tilde{\eta} \right] + \frac{\varepsilon}{3} \partial_x \left[ H_r^3 (\partial_x \tilde{\kappa}_{(1)} +\partial_x \tilde{\eta}) \right] =0, \end{align}

where $\kappa _{(0)}$ is the base flow curvature, i.e. (3.3) evaluated for the base flow $H_r(y)$, $\kappa _{(0)}= ({\mathrm {d}^2H_r}/{\mathrm {d}y^2})/{(1+({\mathrm {d}H_r}/{\mathrm {d}y})^2)^{3/2}}$. Furthermore, $\tilde {\kappa }_{(1)}$ is the first-order term of the curvature, i.e. the Jacobian of the curvature evaluated in the base flow and applied to $\tilde {\eta }$ ($\tilde {\kappa }_{(1)}= [\partial _{\tilde {\eta }}\kappa (H_r)]{\tilde {\eta }}$). The full expression of the operator $\partial _{\tilde {\eta }}\kappa (H_r)$ is reported in appendix A. Deriving this expression with respect to $x$ and $y$, we obtain $\partial _x \tilde {\kappa }_{(1)}=\mathrm {i}k_x {\kappa }_{(1)}(y)\exp ({\mathrm {i}(k_x x -\omega t)})$ and $\partial _y \tilde {\kappa }_{(1)}=({\mathrm {d} {\kappa }_{(1)}}/{\mathrm {d}y})(y)\exp ({\mathrm {i}(k_x x -\omega t)})$.

At $ {O}(1)$ the base flow equation is recovered, while at $ {O}(\varepsilon )$ one obtains the following evolution equation for the perturbation:

(4.3)\begin{align} &-\mathrm{i} \omega \eta+\mathrm{i}k_xu H_r^2\eta + \frac{1}{3} \frac{\mathrm{d}}{\mathrm{d} y} \left[ 3H_r^2 \left(\frac{\mathrm{d}H_r}{\mathrm{d} y}+ \frac{\mathrm{d}\kappa_{(0)}}{\mathrm{d} y}\right)\eta \right.\nonumber\\ &\quad +\left. H_r^3 \left(\frac{\mathrm{d}{\kappa}_{(1)}}{\mathrm{d} y}+ \frac{\mathrm{d}\eta}{\mathrm{d} y}\right) \right] - \frac{1}{3} k_x^2\left[ H_r^3 \left({\kappa}_{(1)} +\eta\right) \right] =0, \end{align}

which is the dispersion relation $D_r(\omega , k_x)=0$. The base flow $H_r(y)$ can be perturbed by (i) imposing the streamwise wavenumber $k_x\in \mathbb {R}$ and looking at the temporal evolution through the complex frequency $\omega \in \mathbb {C}$ (temporal stability analysis) or (ii) imposing a temporal forcing of real frequency $\omega$ and looking at the spatial amplification of the perturbation, embodied by the complex spatial wavenumber $k_x\in \mathbb {C}$ (spatial stability analysis).

The numerical implementation of (4.3) is performed in MATLAB by a spectral collocation Fourier method. Once discretized, the eigenfunction problem (4.3) becomes an eigenvalue problem. The temporal and spatial stability analyses are respectively solved using the built-in MATLAB functions eig and polyeig. Numerical convergence is achieved for 100 collocations points. A preparatory analysis on the numerical rivulet profile $H_r(y)$ used as base flow for the stability analysis revealed a variation of the eigenvalues of the order of the numerical discretization, as long as $t>5000$.

4.3. Temporal stability analysis

In this section, we report the results for the temporal stability analysis. Positive (respectively negatives) values of the temporal growth rate $\mathrm {Im}(\omega )=\omega _i$ denote unstable (respectively stable) wavenumbers. A preliminary analysis on the spectrum revealed that all the eigenvalues have negative $\omega _i$ for all $k_x$, except one that is analysed in the following.

In figure 5(a) we report the variation of $\omega _i$ with $k_x$, for different values of $u$. The dispersion relations are characterized by a local maximum associated with the dominant wavenumber, and by a value of the wavenumber beyond which the temporal growth rate is negative (the cutoff wavenumber), i.e. perturbations with wavenumber larger than the cutoff are damped. Rivulets are strongly stabilized as the value of $u$ increases. For $u=1$ the growth rate $\omega _i$ presents its maximum value at a dominant wavenumber close to $k_x=0.56$, while the cutoff wavenumber $k_x^{cut}=0.8$. An increase of $u$ quickly quenches large wavenumbers. Both the dominant growth rate and the cutoff wavenumber decrease. At $u=5$, $k_x^{cut} \sim 10^{-2}$, with ${\rm max}(\omega _i) \sim 10^{-3}$. For these values of $u$, the unstable wavelengths are of the order of one hundred reduced capillary lengths. The real frequency $\mathrm {Re}(\omega )=\omega _r$ increases slightly less than linearly with $k_x$ (figure 5b). The resulting phase velocities $\omega _r/k_x$ increase as $u$ increases.

Figure 5.

(a) Temporal growth rate $\omega _i$ and (b) real frequency $\omega _r$ as functions of the streamwise wavenumber $k_x$, from the temporal stability analysis, for $u=1$ (blue line), $u=1.5$ (red line), $u=2$ (yellow line), $u=2.5$ (purple line), $u=3$ (green line), $u=5$ (light blue line).

In figure 6(a) we show the real (dashed-dotted line) and imaginary (dashed line) parts of the mode $\eta (y)$ for the dominant wavenumber $k_x=0.5$, normalized by the maximum modulus ${\rm max}(|\eta |)$, for $u=1.5$. The mode is non-zero only in the steady central lobe region. For the same value of $u$, in figure 6(b) we report a three-dimensional plot of the linear combination of the base flow $H_r(y)$ (extended in the $x$ direction along which it is invariant) with the mode at the dominant wavenumber (normalized by the maximum modulus), i.e. $h(x,y)=H_r(y)+{A}\,\mathrm {Re}(\eta (y) \exp (\mathrm {i} k_x x))$, with ${A}=0.25$ an arbitrary amplitude for visualization purposes. The resulting pattern is characterized by rivulets that carry lenses. The temporal dependence of the mode, which is not represented in figure 6(b), is characterized by a growing amplitude $\exp (\omega _i t)$ and by an oscillating behaviour $\exp (\mathrm {i} \omega _r t)$. The presence of a non-zero real part of $\omega$ (figure 5b) implies that the perturbations are oscillating in time at fixed locations. This effect is related to the advection as lenses are travelling along the streamwise direction.

Figure 6.

Temporal stability analysis, $u=1.5$. (a) Real (solid line) and imaginary (dashed line) parts of the eigenvector $\eta (y)$, for the dominant wavenumber $k_x=0.5$, normalized by the maximum modulus. (b) Linear combination of the base flow $H_r(y)$ (extended in the $x$ direction along which it is invariant) with the mode at the dominant wavenumber (normalized with the maximum modulus), i.e. $h(x,y)=H_r(y)+{A}\,\mathrm {Re}(\eta (y) \exp (\mathrm {i} k_x x))$. $A=0.25$ is an arbitrary amplitude for visualization purposes.

The stability analysis reveals the occurrence of a secondary instability of the saturated and one-dimensional rivulets, which is located in the steady central lobe and leads to a pattern characterized by lenses that travel on the rivulets. Nevertheless, an increase in the advection $u$ induces a very strong stabilization and only very large wavelengths remain slightly unstable. The stabilization is related to the advection term. In particular, perturbations in regions of different thickness experience different advection velocities, proportional to $u H_r^2$ (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). Regions of higher thickness travel faster than regions of lower thickness, leading to a steepening of the interface profile and eventually to a capillary levelling of perturbations. This steepening–levelling mechanism is all the more pronounced as $u$ is large. Small wavelengths, which present high interface gradients, are progressively stabilized with $u$, leading to a cutoff wavelength of the order of $10^2 {l}_c^*$ at $u=5$. In the numerical simulation of figure 3(a) the resulting pattern does not show any appreciable streamwise perturbations since the cutoff wavelength ($L_c =2\pi /k_x^{cut} \approx 2 \times 10^2$) is of the order of the maximum acceptable wavelength fitting in the domain. These results are consistent with the experimental observations of Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) when large values of $u$ are considered. For $u>3$, only very large wavelengths are unstable and they are eventually suppressed because of the size of the experiment ($2 \times 10^2 l_c^* < L_x < 3 \times 10^2 l_c^*$).

4.4. Spatial stability analysis

In this section, we study the spatial stability properties of the rivulet base flow $H_r(y)$ introduced in § 4.1. The saturated rivulet profile is perturbed with a temporal harmonic perturbation of real frequency $\omega =\omega _r$ and we look for the spatial evolution of the perturbation, in terms of spatial growth rate $-\mathrm {Im}(k_x)$ and streamwise wavenumber $\mathrm {Re}(k_x)$ through the dispersion relation $D_r(k_x,\omega )$ ((4.3)). Positive values of the spatial growth rate denote unstable configurations associated with downstream propagating waves (Huerre & Rossi Reference Huerre and Rossi1998; Schmid et al. Reference Schmid, Henningson and Jankowski2002; Gallaire & Brun Reference Gallaire and Brun2017). The spectrum is characterized by only one unstable mode associated with downstream propagating waves, which is described in the following.

In figure 7(a) we report the spatial growth rate $-\mathrm {Im}(k_x)$ as a function of $\omega$. The spatial growth rate presents a behaviour similar to the temporal growth rate of § 4.3, i.e. characterized by a maximum (dominant) value and a cutoff frequency beyond which perturbations are damped. The dominant value of $-\mathrm {Im}(k_x)$ strongly decreases with $u$, while its associated dominant frequency presents a non-monotonic behaviour. The same non-monotonic behaviour is observed in the cutoff frequency. The streamwise wavenumber $\mathrm {Re}(k_x)$ (figure 7b) shows, to a good approximation, a linear dependence with $\omega$. For fixed $\omega$, the value of $\mathrm {Re}(k_x)$ decreases with $u$.

Figure 7.

(a) Spatial growth rate and (b) streamwise wavenumber as functions of $\omega$, from the spatial stability analysis, for $u=1$ (blue line), $u=1.5$ (orange line), $u=2$ (yellow line), $u=2.5$ (purple line), $u=3$ (green line). The circles identify the values of the spatial growth rate obtained by the Gaster transformation.

The results of the spatial stability analysis are compared with those of the temporal stability analysis, suitably rescaled by the Gaster transformation (Gaster Reference Gaster1962), valid for strongly convectively unstable systems (see appendix B for details). Within this approximation, from the temporal stability analysis of § 4.3 (labelled with $(T)$) we retrieve the spatial stability analysis properties (labelled with $(S)$) through the relations

(4.4a–c)\begin{equation} \omega_r(S)=\omega_r(T), \quad \mathrm{Re}(k_x(S))=\mathrm{Re}(k_x(T)), \quad \mathrm{Im}(k_x(S))=-\frac{\omega_i(T)}{\dfrac{\partial \omega_r}{\partial k_x}(T)}. \end{equation}

The results of the Gaster transformation (4.4a–c) (circles) are in good agreement with the spatial stability analysis results in figure 7(a), for $u>1$. In appendix B we report the results for $u<1$, where the Gaster transformation prediction deviates from the spatial stability analysis results.

In the following, we experimentally investigate the link between the spatial stability analysis and the observable dynamics.

6.1. Experimental observation

In this section, we introduce a qualitative visualization of the evolution of a localized perturbation in the film thickness. The experimental apparatus is set without any inlet perturbation devices shown in figure 1(a). When high inclination angles and low flow rates are considered (i.e. high values of $u$), we experimentally observe a large region characterized by a uniform flat film where thickness perturbations from the lateral boundaries of the experiment do not penetrate (Lerisson et al. Reference Lerisson, Ledda, Balestra and Gallaire2020). In this region, we trigger the destabilization with a thickness perturbation by blowing a puff of air with a syringe. The whole field is then projected on a screen via the shadowgraph technique and captured with a camera.

In figure 12 we show the evolution of the perturbation with time. The initially localized perturbation is advected away in the streamwise direction with a constant velocity and spreads in the domain. The perturbation phase lines are concentric circles in the upstream part of the response. Nevertheless, the isotropy disappears in the downstream part. The shadowgraph reveals that the phase lines tend to be parallel to the streamwise direction, the effect becoming more and more evident as the time increases.

Figure 12.

Shadowgraph visualization of an experimental impulse response, for $\theta =20^{\circ }$ and $h_N=1292\ \mathrm {\mu }$m, i.e. $u=5.45$. Time increases going to the right and each snapshot is separated by $15$ s. (a) $t=0$ s, (b) $t=15$ s, (c) $t=30$ s and (d) $t=45$ s.

The presence of phase lines aligned with the streamwise directions suggests the existence of a wavefront characterized by streamwise structures, i.e. rivulets, when the flat film is perturbed using an impulse thickness perturbation. The selection of a streamwise wavefront is not related to the boundaries of the thin film in the experiment, i.e. the rivulets selection is intrinsic.

6.2. Numerical observation

Inspired by this experimental observation, in this section we numerically simulate the impulse response, via (3.2), for the same values of angle and flow rate used in the shadowgraph of figure 12, i.e. $u=5.45$, in a double-periodic domain. The initial condition is taken in the form

(6.1)\begin{equation} h(x,y,0)=1+{A}\exp\left(-\frac{x^2+y^2}{2}\right), \end{equation}

where $A=10^{-2}$. In figure 13 we plot the time evolution of the response in the moving reference frame $(\xi =x-u t,y)$, from $t=0$ to $t=90$. In the moving reference frame, the response progressively invades the domain from the initial impulse location. At $t=30$, we observe circular phase lines. At $t=60$ the response loses its isotropy in the downstream part. At $t=90$ streamwise structures are dominant in the downstream front of the response and they are also observable upstream.

Figure 13.

Impulse response, $\theta =20^{\circ }$ and $h_N=1292\ \mathrm {\mu }$m ($u=5.45$). The time increases from left to right and the time step is $30$. Results are reported in the moving reference frame at the linear advection velocity $(\xi =x-u t,y)$. (a) $t=0$, (b) $t=30$, (c) $t=60$, (d) $t=90$.

In the moving reference frame, the response spreads from the initial impulse location, meaning that in the fixed reference frame, the response is advected downstream at the linear advection velocity $u$. The numerical evolution qualitatively agrees with the experimental observation of § 6.1. We first observe the evolution of the impulse response into an isotropic pattern. At large times, the response mostly evolve towards streamwise structures. However, the complicated form of (3.2), including nonlinear advection, hydrostatic pressure distribution and capillary effects, does not allow one to identify the physical mechanisms that lead to the emergence of streamwise structures observed in figures 12 and 13. Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020) furthermore observed that the rivulet propagation and growth are well described by the linear stability analysis of the flat film even at large amplitudes of the thickness perturbation, beyond the expected validity of the linear theory. Hereafter, we study the origin of the selection of rivulet structures by the linear and weakly nonlinear dynamics.

6.3. Linear response

Upon introduction of the decomposition $h=1+\varepsilon \eta$ ($\varepsilon \ll 1$) in (3.2), the linearized equation at $ {O}(\varepsilon )$ for the evolution of the thickness perturbation $\eta$ with respect to the flat film reads

(6.2)\begin{equation} \partial_{t} \eta+u \partial_{x} \eta+\tfrac{1}{3}[\boldsymbol{\nabla}^{2} \eta+\boldsymbol{\nabla}^{4} \eta]=0.\end{equation}

The dispersion relation is recovered by introducing the normal mode decomposition $\eta \propto \exp [\mathrm {i}(\boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {x}-\omega t)]$ , with $\boldsymbol {k}=(k_{x}, k_{y})$, where $k_x$ and $k_y$ denote respectively the streamwise and spanwise wavenumbers

(6.3)\begin{equation} \omega=u k_{x}+\frac{\mathrm{i}}{3}({k}^{2}-{k}^{4}),\end{equation}

where $k=\sqrt {k_x^2+k_y^2}$. The dispersion relation $D(\omega ,k_x,k_y)=0$ is characterized by an isotropic temporal growth rate $\omega _i$, as shown in figure 14(a). The temporal frequency $\omega _r$ is linear in $k_x$ and does not depend on $k_y$.

Figure 14.

(a) Temporal growth rate $\omega _i$ as a function of $k_x$ and $k_y$. (b) Linear impulse response, $u=5.45$. The vertical and horizontal axes are respectively the streamwise and spanwise directions. From left to right: $t=30$, $t=60$, $t=90$. Results are reported in the moving reference frame at the linear advection velocity $(\xi =x-u t,y)$.

The initial condition for the numerical simulation is the thickness perturbation $\eta (x,y,0)={{A}} \exp (-x^2/2 -y^2/2)$, where ${A}=10^{-2}$. The linear numerical simulation results for $u=5.45$, in the moving reference frame $(\xi =x-u t,y)$, are presented in figure 14(b). As time increases, the perturbation spreads in concentric circles from the initial impulse location. Similarly to the nonlinear simulation of figure 13, the response is advected away at the linear advection velocity $u$, in the fixed reference frame. The results can be rationalized considering the dispersion relation of (6.2). The wavepacket is non-dispersive since $\omega _r$ is linear in $k_x$. This means that there is no distortion of the wavepacket. Since the growth is isotropic, concentric circles invade the domain and at the same time are advected downstream with constant velocity $\omega _r/k_x=u$. Higher values of $u$ imply faster advection velocities. In the moving reference frame $(\xi ,y)$, the (6.2) reads

(6.4)\begin{equation} \partial_t \eta +\tfrac{1}{3}[ \boldsymbol{\nabla}_{\xi y}^2 \eta + \boldsymbol{\nabla}_{\xi y}^4 \eta ]=0,\end{equation}

where $\boldsymbol {\nabla }_{\xi y}$ operates in the reference frame $(\xi ,y)$. In this reference frame, the response spreads in perfectly isotropic concentric circles without being advected away. The linear dynamics agrees well with the early-time evolution of the nonlinear simulation shown in figure 13, when the amplitude of the perturbations is still very small. However, since the linearized dynamics is not able to capture the anisotropy of the pattern observed in the nonlinear simulation, we propose next a weakly nonlinear study.

6.4. Weakly nonlinear response: the Nepomnyashchy equation

We consider a weakly nonlinear model for the flow of a thin film on the underside of an inclined planar substrate. Following Kalliadasis et al. (Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012), the derivation is based on a multiple scale approach combined with an asymptotic expansion. Under the assumption of small interfacial disturbances and $u= {O}(1)$, the weakly nonlinear dynamics for a thickness perturbation $\eta$ with respect to the flat film reads

(6.5)\begin{equation} \partial_t \eta +2u \eta \partial_\xi \eta+\tfrac{1}{3} [ \boldsymbol{\nabla}_{\xi y}^2 \eta + \boldsymbol{\nabla}_{\xi y}^4 \eta ]=0, \end{equation}

where $\boldsymbol {\nabla }_{\xi y}$ operates in moving the reference frame $(\xi ,y)$. The equation is formally analogous to the Nepomnyashchy equation (Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2012). We consider the evolution of the thickness perturbation $\eta$ starting from a Gaussian impulse $\eta (\xi ,y,0)={A} \exp (-\xi ^2/2 -y^2/2)$ (${A}=10^{-2}$), in analogy with the linear simulation.

In figure 15(a) we report the thickness perturbation evolution. The initial localized perturbation spreads in the domain and is always centred in the vicinity of the initial impulse position, because of the moving reference frame. At $t=30$ the perturbation has spread isotropically in the domain. Nevertheless, at $t=60$, streamwise structures arise. At $t=90$, the streamwise structures have invaded most of the perturbation region.

Figure 15.

Case $u=5.45$. (a) Impulse response in the moving reference frame $(\xi =x-u t,y)$ from the weakly nonlinear model and (b) its two-dimensional Fourier energy spectrum. From left to right: $t=0$, $t=30$, $t=60$, $t=90$.

In figure 15(b) we show the two-dimensional Fourier energy spectrum of $\eta$, normalized by its maximum value. Since we are considering a real signal, the Fourier spectrum is symmetric with respect to the $k_x$ and $k_y$ axes. We thus report only the values in the first quadrant ($k_x>0,k_y>0$). At $t=0$, we observe the Fourier spectrum of a Gaussian impulse, which is a Gaussian centred around ($k_x=0,k_y=0$), i.e. the initial spectrum is isotropic. At $t=30$ the energy is located in a region around $\sqrt {k_x^2+k_y^2}=1/\sqrt {2}$. As time progressively increases, the energy concentrates towards $(k_x=0,k_y=1/\sqrt {2})$.

Initially, the response is characterized by an isotropic pattern, reminiscent of the linear growth that is experienced in the first stages of the perturbation growth. As the amplitude becomes sufficiently large, the spectrum shows that the energy is focusing on the axis $k_x=0$, i.e. streamwise structures are selected. The emergence of streamwise structures agrees well with the results of the fully nonlinear simulation and with the experimental observation. Moreover, the spectrum is localized around $k_y=1/\sqrt {2}$, the most amplified wavelength predicted by the flat film dispersion relation ((6.3)), and rivulet structures are growing exponentially. Thus, the dynamics of pure streamwise structures stays linear, even in the weakly nonlinear regime.

The origin of the selection of rivulet structures is identified in the weakly nonlinear advection term $2u \eta \partial _\xi \eta$, which acts in indirect manner to favour rivulet structures while damping all other orientations. The weakly nonlinear model of (6.5) is formally analogous to the linear model of (6.4), except for the weakly nonlinear advection term. It should be noticed that this term influences the dynamics of streamwise-inhomogeneous structures only, on which it is seen to have a damping effect. The nonlinear advection term embodies the difference in the perturbation advection velocity in regions of different thickness and is known to create wave steepening (Babchin et al. Reference Babchin, Frenkel, Levich and Sivashinsky1983). The emerging steep gradients are damped by surface tension effects, levelling therefore the non-streamwise structures. In conclusion, the most unstable solution in the weakly nonlinear regime is the one in which the capillary damping is reduced the most, as the term $2u \eta \partial _\xi \eta$, responsible of wave steepening, vanishes.

When only streamwise structures are present, the advection term disappears and the weakly nonlinear model is formally analogous to the linear equation in the moving reference frame (6.4). Consequently, the response of streamwise structures is linear up to second order in the perturbation.

In conclusion, the weakly nonlinear dynamics gives an insight into the origin of the emergence of rivulet structures: the latter are the only ones screened from the action of the difference in the advection. The dynamics of pure streamwise structures remains linear even in the weakly nonlinear regime, thus explaining the agreement between the linear prediction and the experimental measurements at large amplitudes observed in Lerisson et al. (Reference Lerisson, Ledda, Balestra and Gallaire2020). At late times, rivulets eventually invade the perturbation region. In the case of steady inlet forcing (figure 2) rivulets invade the whole domain and steady and streamwise saturated rivulet structures emerge downstream, as a result of the weakly nonlinear dynamics. As seen in the previous sections, rivulets may eventually destabilize through a secondary instability, resulting in travelling lenses. In both the emergence and the stability of rivulets, the differences in advection in regions of different thickness is crucial.