# Instability and rupture of sheared viscous liquid nanofilms

Vira Dhaliwal  Mechanics Division, Department of Mathematics, University of
Oslo, 0316 Oslo, Norway  Christian Pedersen  Mechanics Division, Department of
Mathematics, University of Oslo, 0316 Oslo, Norway  Kheireddine Kadri
Laboratoire PIMM, Arts et Métiers, CNRS, Hesam Université, 151 boulevard de
l’Hopital, Paris, France  Guillaume Miquelard-Garnier  Laboratoire PIMM, Arts
et Métiers, CNRS, Hesam Université, 151 boulevard de l’Hopital, Paris, France
Cyrille Sollogoub  Laboratoire PIMM, Arts et Métiers, CNRS, Hesam Université,
151 boulevard de l’Hopital, Paris, France  Jorge Peixinho  Laboratoire PIMM,
Arts et Métiers, CNRS, Hesam Université, 151 boulevard de l’Hopital, Paris,
France  Thomas Salez<EMAIL_ADDRESS>Univ. Bordeaux, CNRS, LOMA, UMR
5798, F-33405, Talence, France  Andreas Carlson<EMAIL_ADDRESS>Mechanics Division, Department of Mathematics, University of Oslo, 0316 Oslo,
Norway

###### Abstract

Liquid nanofilms are ubiquitous in nature and technology, and their
equilibrium and out-of-equilibrium dynamics are key to a multitude of
phenomena and processes. We numerically study the evolution and rupture of
viscous nanometric films, incorporating the effects of surface tension, van
der waals forces, thermal fluctuations and viscous shear. We show that thermal
fluctuations create perturbations that can trigger film rupture, but they do
not significantly affect the growth rate of the perturbations. The film
rupture time can be predicted from a linear stability analysis of the
governing thin film equation, by considering the most unstable wavelength and
the thermal roughness. Furthermore, applying a sufficiently large
unidirectional shear can stabilise large perturbations, creating a finite-
amplitude travelling wave instead of film rupture. In contrast, in three
dimensions, unidirectional shear does not inhibit rupture, as perturbations
are not suppressed in the direction perpendicular to the applied shear.
However, if the direction of shear varies in time, the growth of large
perturbations is prevented in all directions, and rupture can hence be
impeded.

††preprint: APS/123-QED

##  I Introduction

Thin liquid films are found in many biological systems such as the human tear
film [1, 2] as well as in modern micro- and nanofabrication processes such as
multilayer coextrusion of polymers [3, 4]. The stability of such films is
often an important consideration, with hole formation due to film rupture
often being undesirable [4, 5]. In other applications, the rupture and
dewetting of thin liquid films is intended, and can be manipulated in order to
fabricate patterned materials [6, 7].

It has been known for over half a century that the rupture of a thin liquid
film can be caused by the amplification of small interfacial perturbations by
long-range intermolecular van der Waals forces. In early works by Vrij [8] and
Sheludko [9], the attractive surface interaction due to the van der Waals
forces was represented by a thickness-dependent potential, and thermodynamic
approaches were used to derive a critical wavelength above which perturbations
to the flat film profile are unstable. This potential was later incorporated
in a hydrodynamic model, which enabled derivation of the growth rate of a
surface perturbation as a function of its wavelength, thus allowing the
estimation of the rupture time [10].

The hydrodynamic model has been simplified by using the lubrication
approximation to derive the so-called thin film equation, which describes the
spatiotemporal evolution of the liquid film height with a single, highly
nonlinear partial differential equation [11, 12]. This has led to numerous
analytical and numerical works which have shown that the growth of surface
perturbations can be separated into two regimes: an initial linear regime
during which the film height has not yet deviated significantly from its
initial value, and a subsequent nonlinear regime during which the growth of
perturbations is greatly accelerated [11, 13]. Zhang and Lister [14] showed
that the nonlinear late-stage dynamics of rupture due to a disjoining pressure
are governed by a similarity solution, with the minimum film thickness rapidly
decreasing according to a power law. Deviation from an idealized disjoining
pressure derived from the Lennard-Jones potential has been shown to produce
discretely self-similar solutions [15]. Such a process was also studied for
analogous elastic interfaces [16]. Beyond the initiation of film rupture, thin
film models have also been used to describe the droplet patterns formed during
dewetting [17].

Van der Waals forces leading to liquid film rupture on solid surfaces only
become significant when the film thickness is on the order of tens of
nanometers [18]. At these length scales and for common fluids at ambient
temperatures, thermally driven molecular motion may cause significant
fluctuations of the film height. Thus, it is natural to consider what role
these microscopic thermal fluctuations play during the rupture process. In
order to take these into account, a stochastic version of the thin-film
equation was derived from the Navier-Stokes equation with an additional random
stress tensor [19, 20]. This formulation has been used to describe other thin
film processes such as the spreading of a viscous bump under an elastic plate
[19, 21, 22] and the transport of solutes through nanopores [23]. Numerical
and theoretical studies using the stochastic thin film equation in two
dimensions have shown that random thermal fluctuations generate perturbations
to an initially flat film profile, which eventually coarsen and approach the
wavelength of maximum growth obtained from the linear stability analysis of
the fluctuation-free system [24, 25, 26]. Simulations have also shown that
increasing the fluctuation intensity decreases the rupture time and creates a
more nonuniform pattern of droplet sizes after dewetting [27, 25, 28].
Nevertheless, these thermal effects seem important only at the early stage of
the rupture process, with the late-stage behaviour being unaffected due to the
dominance of the van der Waals force [29].

One way to modify the dynamics of a thin liquid film is by applying a shear
flow to it – a situation widely encountered in industrial applications.
Interestingly, shear is expected to dampen the amplitude of thermal interface
fluctuations [30, 31, 32], which might in turn influence the rupture time,
i.e., the time it takes for the interface to touch down on the substrate.
Numerical studies of deterministic thin film rupture in two dimensions have
indeed shown that the presence of unidirectional shear can delay rupture when
the interface has unstable perturbations [33, 34, 35]. To the authors’
knowledge, experiments have not yet demonstrated the theorized rupture-
suppressing effect of shear, but have shown that shear changes the morphology
of the holes created when a polymer film dewets [36].

In this article, we solve numerically the stochastic thin film equation to
improve our understanding of how shear affects the rupture dynamics of
nanometric liquid films. By studying the film rupture for various combinations
of fluctuation intensity, shear rate, and film thickness, we delineate the
mechanisms by which rupture is affected in both two dimensions (2D) and three
dimensions (3D).

##  II Mathematical modeling and numerical methods

Figure 1: Schematic of the studied physical system. A thin layer of viscous
fluid with thickness $h(x,y,t)$ rests upon a flat solid substrate. A positive
disjoining pressure, $\Pi(x,y,t)$, attracts the free surface to the solid
substrate. The interface dynamics are also affected by shear and thermal
fluctuations.

###  II.1 Stochastic thin film equation

We consider the dynamics of a thin viscous liquid film as shown in Fig. (1).
The dynamics of the film is represented by the spatiotemporal evolution of its
height $h(x,y,t)$. For an initially flat film with height $h_{0}$, we define
the perturbation to the film height as $\delta h=h(x,y,t)-h_{0}$. At any point
in time, the minimum height of the film is denoted $h_{min}.$ Since we are
interested in films where the height $h$ is much smaller than the horizontal
length scale, the viscous flow profile in the fluid layer can be described by
the lubrication approximation [37]. A flux balance for a slice of fluid in the
film then gives us the thin film equation describing how $h(x,y,t)$ varies
with the pressure profile $p(x,y,t)$ [12]. In our system, the pressure
$p(x,y,t)$ has contributions arising from the Laplace pressure due to the
curvature of the free surface, the disjoining pressure $A^{*}/[6\pi
h^{3}(x,y,t)]$ resulting from the van der Waals interactions [38], where
$A^{*}$ is the Hamaker constant, and an additional stochastic stress arising
from the thermal fluctuations in the fluid. The shear force is implemented in
the tangential stress boundary condition at the film surface, taking the form
$\mu\partial\textbf{u}/\partial z\rvert_{z=h(x,y,t)}=\bm{\tau}$, where
$\textbf{u}(x,y,z,t)$ is the instantaneous velocity, $\mu$ is the dynamic
viscosity, and $\bm{\tau}=\tau\bm{e}_{\tau}$ is the specified shear stress at
the interface, with magnitude $\tau$ and oriented along the unit vector
$\bm{e}_{\tau}$. Assuming incompressible flow of a Newtonian fluid, small
Reynolds number and small slope of the film profile, $\nabla h(x,y)\ll 1$, the
thin-film equation reads:

$\displaystyle 0=\frac{\partial h(x,y,t)}{\partial
t}+\nabla\cdot\left[\underbrace{\frac{\gamma}{3\mu}h^{3}(x,y,t)\nabla\nabla^{2}h(x,y,t)}_{\text{Laplace
pressure}}+\underbrace{\frac{A^{*}}{6\pi\mu h(x,y,t)}\nabla
h(x,y,t)}_{\text{Disjoining pressure}}\right]$
$\displaystyle+\underbrace{\frac{\tau}{\mu}\bm{e}_{\tau}\cdot h(x,y,t)\nabla
h(x,y,t)}_{\text{Shear
force}}+\underbrace{\sqrt{\frac{k_{B}T}{6\mu}}\nabla\cdot\left[h^{3/2}(x,y,t)\bm{\eta}(x,y,t)\right]}_{\text{Stochastic
force}},$ (1)

where $\gamma$ is the surface tension coefficient, $k_{B}$ is Boltzmann’s
constant and $T$ is the temperature. The last term in the thin-film equation
accounts for the thermal fluctuations in the fluid, where $\bm{\eta}(x,y,t)$
is a random vector in the $(x,y)$-plane, the two components, $\eta_{i}$ with
$i=x,y$, of which being independent delta-correlated spatiotemporal Gaussian
noises with null averages [26, 19], i.e. $\langle\eta_{i}(x,y,t)\rangle=0$ and
$\langle\eta_{i}(x,y,t)\eta_{j}(x^{\prime},y^{\prime},t^{\prime})\rangle=\delta_{ij}\delta(x-x^{\prime})\delta(y-y^{\prime})\delta(t-t^{\prime})$,
where $\delta_{ij}$ is the Kronecker symbol and $\delta$ is the Dirac
distribution.

We nondimensionalize Eq. (1) by introducing the scaling relations:

$h=\bar{h}\sqrt{\frac{A^{*}}{2\pi\gamma}},\hskip
10.0ptx=\bar{x}\sqrt{\frac{A^{*}}{2\pi\gamma}},\hskip
10.0pty=\bar{y}\sqrt{\frac{A^{*}}{2\pi\gamma}},\hskip
10.0ptt=\bar{t}\sqrt{\frac{9A^{*}\mu^{2}}{2\pi\gamma^{3}}},$

$\bm{\eta}=\bar{\bm{\eta}}\left(\frac{8\pi^{3}}{9}\frac{\gamma^{5}}{A^{*3}\mu^{2}}\right)^{1/4},\hskip
10.0pt\bm{\tau}=\bar{\tau}\bm{{e_{\tau}}}\sqrt{\frac{2\pi\gamma^{3}}{9A^{*}}}$

where the bars indicate dimensionless quantities. Both vertical and horizontal
lengths are nondimensionalized by a physical length scale,
$\sqrt{A^{*}/(2\pi\gamma)}$, representing the characteristic film thickness at
which the disjoining and capillary effects balance. We are mainly interested
in the behavior of films with dimensionless initial thickness,
$\bar{h}_{0}\approx 1$. Films with $\bar{h}_{0}\gg 1$ will not be
significantly affected by the long-range van der Waals forces, while films
with $\bar{h}_{0}\ll 1$ will rupture almost instantaneously. The time is
normalized by the aforementioned length scale divided by the capillary
velocity $\gamma/(3\mu)$. The dimensionless stochastic thin film equation thus
reads:

$\displaystyle 0=\frac{\partial h(x,y,t)}{\partial
t}+\nabla\cdot\left[h^{3}(x,y,t)\nabla\nabla^{2}h(x,t)+\frac{1}{h(x,t)}\nabla
h(x,y,t)\right]$ $\displaystyle+B\bm{e}_{\tau}\cdot h(x,y,t)\nabla
h(x,y,t)+Q\nabla\cdot\left[h^{3/2}(x,y,t)\bm{\eta}(x,y,t)\right],$ (2)

where we have dropped the bars for simplicity, and introduced two
dimensionless numbers:

$B=\bar{\tau}=\tau\sqrt{\frac{9A^{*}}{2\pi\gamma^{3}}},$ (3)
$Q=\sqrt{\frac{\pi k_{B}T}{A^{*}}},$ (4)

representing the dimensionless shear force, and the dimensionless thermal
roughness (or the ratio between thermal and disjoining energies),
respectively. Along with the initial film height $h_{0}$, these parameters
define the thin-film dynamics.

We note that, when we model thin films in 2D (i.e. when there is invariance in
the $y$-direction), the fluctuation vector $\bm{\eta}$ needs to be adapted due
to the inherently three-dimensional nature of thermal fluctuations. Eq. (II.1)
remains valid, but $\bm{\eta}(x,t)$ is now a random vector in the
$x$-direction only, the single component $\eta$ of which being spatiotemporal
Gaussian noise with $\langle\eta(x,t)\rangle=0$ and
$\langle\eta(x,t)\eta(x^{\prime},t^{\prime})\rangle=\delta(x-x^{\prime})\delta(t-t^{\prime})$.
As a consequence, $Q$ is modified in 2D, as:

$Q_{2D}=\sqrt{\frac{\pi k_{B}T}{\bar{w}A^{*}}},$ (5)

where $\bar{w}=w\sqrt{\frac{2\pi\gamma}{A^{*}}}$, with $w$ a new length scale
representing the width in the y-direction.

###  II.2 Finite element solver

We solve the dimensionless stochastic thin film equation (Eq. (II.1)) using
the finite element method. In all simulations, the initial condition is a flat
film with dimensionless thickness $h_{0}$. Both in 2D and 3D, periodic
boundary conditions are imposed in the horizontal directions. The domain size
in the horizontal directions is always much greater than the thickness of the
film to ensure the validity of the lubrication approximation.

To solve Eq. (2), the order of the partial differential equation is first
reduced by introducing the film curvature $\nabla^{2}h(x,y,t)$ as a separate
variable. Eq. (II.1) then simplifies to a system of two coupled second order
partial differential equations. These are then expressed in a weak form where
boundary terms disappear due to the periodic boundary conditions. The scalar
fields $\nabla^{2}h(x,y,t)$ and $h(x,y,t)$ are discretized with linear
elements and solved using a Newton solver from the FEniCS library [39]. For
the 2D case, the simulations are realized on the domain $x=[0,65]$ with an
equidistant grid spacing $\Delta x=0.01$. In the 3D case, the domain is
$x\times y=[0,64]\times[0,64]$, and the grid spacings are $\Delta x=\Delta
y=0.16$ for simulations used to generate the contour plots. For the data shown
in other figures where many iterations are averaged, the grid spacing is
increased to $\Delta x,\Delta y=0.64$ in order to shorten the simulation time,
and we have checked that the results are insensitive to the change in spatial
resolution. Due to the fact that Eq.(II.1) involves a travelling wave, a
second-order Crank-Nicholson scheme is required for numerical integration in
order to prevent numerical dispersion for large values of $B$. The size of the
time step for all simulations shown is $\Delta t=0.003$, with the exception of
certain simulations involving $B>30$ as well as the adaptive time step
simulations described later in this section and in appendix A.

The stochastic term $\bm{\eta}(x,y,t)$ is implemented in python by assigning
random numbers using the “normal” function in the “random” class of NUMPY
[40]. The values are drawn from a Gaussian distribution with zero mean and a
variance of $1/(\Delta x\Delta t)$ in 2D and $1/(\Delta x^{2}\Delta t)$ in 3D.
At every time step, each component of $\bm{\eta}(x,y,t)$ is assigned a new
value at every point in the mesh. Due to the stochastic nature of
$\bm{\eta}(x,y,t)$, the film dynamics and rupture time vary somewhat between
individual numerical simulations, despite all the input parameters being
identical. In order to obtain statistically robust results for the wavelength
and rupture time, we repeat the simulations and average the results $N=10$
times. Nevertheless, the variance of the rupture time and wavelength is
generally so small that the trends we describe can be observed even for a
single run.

For the non-stochastic simulations of the deterministic version of Eq. II.1,
an initial sinusoidal perturbation of amplitude $\delta h_{0}=0.001$ is
imposed. In order to maintain the perdiodic boundary condition, the extent of
the domain is set to 8 times the wavelength of the perturbation.

Previous works on shear-free films as well as the initial simulations of our
system indicate that the highly nonlinear van der Waals force is completely
dominant during the final stage of rupture when the minimum film height
$h_{min}(t)=min(h(x,y,t))$ approaches zero [14, 29]. Due to the accelerated
dynamics in this regime, a much smaller time step is required in order to
capture the final moments before film rupture. To address this, we perform
simulations with an adaptive time step, the results of which are shown in
appendix A. In those simulations, the time step starts at $\Delta t=0.1$ but
is gradually reduced as rupture accelerates, reaching a minimum value of
$\Delta t\approx 10^{-8}$. The grid spacing is also decreased to $\Delta
x=0.001$ in order to resolve more of the details near the rupture point. As
shown in appendix A, we recover the expected $h_{min}\sim(t_{R}-t)^{1/5}$
power law of Zhang and Lister [14], and also observe that although this stage
of rupture accounts for only a tiny fraction of the total duration of the film
dynamics until rupture, it makes up most of the computational time due to the
reduction in the time step. Our simulations with the adaptive timestep also
show that the late stage before rupture seems to be independent of the thermal
fluctuations and shear. We also note that close to rupture, lubrication theory
cannot be used to describe the local flow [41]. Since we are mainly interested
in understanding the instability growth rather than the details in the already
characterized final instants before rupture, we decide to fix the time step
and stop the simulation when $h_{min}$ reaches an arbitrary threshold
$h^{*}=0.33$, which ensures that we capture the entirety of the early-stage
dynamics without using unnecessary computational resources. We denote the time
at which this occurs as the rupture time $t_{R}$.

###  II.3 Linear stability analysis

From the numerical solution of Eq. (II.1), as will be described below, it is
clear that the thermal fluctuations only play a major role in the evolution of
$h_{min}$ during the first few time steps. As shown in appendix A, it is also
clear that the late stage of pre-rupture film dynamics, when $h_{0}-h_{min}\gg
0.1$, is very short and as such does not make a significant contribution to
the measured rupture time. Thus we presume that the rupture time of the film
is primarily determined by how quickly perturbations grow while they are still
too small for nonlinearities to dominate. Here, we perform a linear stability
analysis on the deterministic version of Eq. (II.1) (i.e. with $Q=0$) [34,
35]. We introduce a small sinusoidal perturbation to an initially flat film,
of the form

$h(x,t)=h_{0}+h^{\prime}e^{ikx+\omega t},$ (6)

where $k$ is the dimensionless angular wave number, $\omega$ is the
dimensionless complex growth rate, and $h^{\prime}$ is the amplitude of the
perturbation with $h_{0}\gg h^{\prime}$. Inserting Eq. (6) into the
deterministic version of Eq. (II.1) and linearising the equation gives us the
following dispersion relation

$\omega=\left({\frac{k^{2}}{h_{0}}-k^{4}h_{0}^{3}}\right)-ik{Bh_{0}},$ (7)

with growth rate $\frac{k^{2}}{h_{0}}-k^{4}h_{0}^{3}$ and wave speed $Bh_{0}$.
From Eq. (7), we expect that the effect of the shear within the linear regime
is simply to cause a horizontal translation in the direction of shear with
speed $Bh_{0}$. Moreover, only perturbations with a wavelength above a
critical wavelength $\lambda_{c}=2\pi h_{0}^{2}$ will grow in time, and the
fastest growing mode has dimensionless wavelength

$\lambda_{d}=2\pi\sqrt{2}h_{0}^{2}.$ (8)

The growth rate associated with this dominant wavelength is

$\omega_{d}=\frac{1}{4h_{0}^{5}}.$ (9)

Presuming that the fluctuations trigger perturbations across the wavelength
spectrum, we expect that $\lambda_{d}$ will be the dominant mode observed in
our numerical results, and thus govern the change in minimum height of the
film. If we only consider the minimum height in Eq. (6), presume that the
growth rate of this corresponds to Eq. (9), and divide by $h_{0}$, we then get
the following expression for how the minumum height of the film varies with
time:

$\frac{h_{min}(t)}{h_{0}}=1-\frac{h^{\prime}}{h_{0}}e^{\omega_{d}t}.$ (10)

We consider rupture to occur at the time $t_{R}$ when $h_{min}$ reaches an
arbitrary threshold value, $h^{*}=0.33$, at which the van der Waals forces
become dominant:

$\frac{h^{*}}{h_{0}}=1-\frac{h^{\prime}}{h_{0}}e^{\omega_{d}t_{R}}.$ (11)

Solving this equation for $t_{R}$ and inserting our expression for the
dominant growth rate in Eq. (9), we then obtain the following estimate for the
film rupture time

$t_{R}=4h_{0}^{5}\ln\left(\frac{h_{0}-h^{*}}{h^{\prime}}\right).$ (12)

##  III Results and Discussion

###  III.1 Two-dimensional case

$(a)$

$(b)$

$(c)$

$(d)$

Figure 2: $(a)$ Profiles of an unstable 2D thin film with $Q=10^{-3}$, $B=0$,
and $h_{0}=1$ at various times. We stop the simulations when the minimum film
height $h_{min}$ reaches an arbitrary threshold value $h^{*}=0.33$
(represented by the dotted black line), as described in section III.1. $(b)$
Time evolution of the minimum film height, $h_{min}$, plotted for simulations
with initial height $h_{0}=1.0$ and different fluctuation intensities, $Q$.
$(c)$ Amplitude of the height perturbation as a function of time for the same
simulations. The dotted black line has a slope corresponding to the prediction
of Eq. (10) with arbitrary amplitude. The coloured dotted lines on the left
indicate the value of $Q$ corresponding to each simulation. The three regimes
mentioned in the text have been identified by eye as follows: the fluctuation-
dominated regime at the begin of the simulations has been shaded red, while
the area shaded in blue represents the nonlinear regime at the end of the
simulation. $(d)$ Mean rupture time $t_{R}$ of a film with $B=0$ as a function
of $Q$ for different values of $h_{0}$. The dashed lines represent Eq. (13)
with no free parameter for each value of $h_{0}$. Error bars represent the
standard deviation for a set of $N=20$ simulations for each data point.

An example of the numerical solution of Eq. (II.1) without shear ($B=0$) in 2D
is shown in Fig. 2$(a)$ in the form of height profiles along $x$ at different
times. Thermal fluctuations eventually give rise to sinusoidal-like
perturbations with a certain wavelength, which then grow until the film
ruptures. For each film height $h_{0}$, the wavelength of the unstable film
profile that develops closely matches the dominant wavelength $\lambda_{d}$
predicted in Eq. (8) and the rate of change of $h_{min}(t)$ closely matches
the value predicted by Eq.(9).

In Fig. 2(b-c), we can see how $h_{min}$ changes in time during individual
simulation runs for different values of $Q$. As discussed above, $Q$ is
proportional to the dimensionless thermal roughness, and we here present
results as $Q$ varies from relatively small values until the limit where
thermal fluctuations are so large that they rapidly cause rupture even without
a disjoining pressure. Our results show that the primary effect of thermal
fluctuations is to instigate a perturbation of a characteristic size $Q$
during the first time steps. Fluctuations with approximately this size
naturally continue to occur at all times, but after $t\approx 10$, the minimum
height begins to decrease exponentially as one would expect in the
deterministic case, and the fluctuations become insignificant in comparison to
the growing sinusoidal perturbation. As shown in Fig. 2$(b)$, the time
evolution of $h_{min}$ as the film approaches rupture is almost identical for
different values of $Q$; the strength of the fluctuations seems to only
determine how far back in time from the moment of rupture the curve starts at.
In Fig. 2$(c)$ one can also see quite clearly that there are three distinct
regimes during rupture: first a fluctuation-dominated phase that initiates a
perturbation (shaded red in Fig. 2$(c)$), then a period of exponential growth
of the roughly sinusoidal perturbation (unshaded in Fig. 2$(c)$), followed by
a phase where the growth toward rupture is greatly accelerated due to the
nonlinear effects of the disjoining pressure term in Eq. (II.1) (shaded blue
in Fig. 2$(c)$). The predominance of the exponential growth period during the
rupture process seen in Fig. 2$(c)$ suggests that the rupture time of a film
can be predicted by Eq. (12). We also see from Fig. 2$(c)$ that the initial
perturbation size, $h^{\prime}$, is approximately equal to $Q$. Setting
$h^{*}=0.33$, which appears to be a reasonable bound of the linear domain (i.e
with exponential growth of the instability), we can rewrite Eq. (12) as

$t_{R}=4h_{0}^{5}\ln\left(\frac{h_{0}-0.33}{Q}\right).$ (13)

In 2$(d)$ we present the numerically measured rupture time from repeated
simulations of our system as a function of $Q$ for various values of $h_{0}$.
The predicted rupture time from Eq. (13), represented by the straight dashed
lines, provides an excellent fit to the data with no adjustable parameters.

$(a)$

$(b)$

$(c)$

$(d)$

Figure 3: $(a-c)$ Profiles of sheared 2D thin films with $Q=10^{-3}$ and
$h_{0}=1$ at different times. The shear is from left to right in each case,
with a non-dimensional size $B$ specified at the top of each figure. The
dotted black line at the bottom indicates the threshold value $h^{*}=0.33$ at
which the simulations are interrupted. $(a)$ $B=9.0$. $(b)$ $B=12.6$ $(c)$
$B=24.0$ $(d)$ Time evolution of the height perturbation $h_{0}-h_{min}(t)$
for different values of dimensionless shear stress $B$. The dotted black line
indicates an exponential behaviour as predicted by Eq. (10).

$(a)$

$(b)$

$(c)$

$(d)$

Figure 4: $(a)$ Rupture time $t_{R}$ as a function of $B$ for two different
values of $Q$. Rupture times are normalized by the rupture time for $B=0$.
Error bars represent the standard deviation for a set of 20 simulations for
each data point. $(b)$ Final value of the minimum film height as a function of
$B$ for values above $B_{crit}$. Error bars represent the standard deviation
for a set of 20 simulations for each data point.$(c)$ Rupture time as a
function of $B$ for deterministic thin films with initial perturbations of
varying wavelength $\lambda$. $(d)$ Plot of the critical shear required to
suppress rupture within the simulated time period as a function of $\lambda$.
The top and bottom of the error bars show the lowest shear that did not
rupture and the largest shear that did rupture, respectively.

We now turn our attention to the effects of shear on the previous
phenomenology. Eq. (II.1) is simulated with $B\neq 0$ and simulations are
stopped either at rupture when $h\leq 0.33$, or when the simulation has
reached a maximum time of $t=100$. In Fig. 3$(a-c)$, we see snapshots of the
film profile under shear with 3 different values of $B$. As predicted by the
linear stability analysis (section II.3), the constant shear term leads to a
simple horizontal translation of the perturbation in the shear direction, with
speed $Bh_{0}$. We observe that the perturbations remain roughly sinusoidal-
like during the exponential growth stage, and maintain a wavelength close to
$\lambda_{d}$ as well as a growth rate of approximately $\omega_{d}$. In the
nonlinear regime, however, they start to deform as the van der Waals forces
become stronger, which gives rise to characteristic asymmetric shapes in the
nonlinear stage. For $B\gtrsim 5$, the change in shape is noticeable, and is
associated with a short period of decreased perturbation growth, which leads
to a slight increase in rupture time, as shown in Figs. 3(a) and 3(d). As
shear is increased to around $B\gtrsim 11$, however, the perturbations begin
to stabilize when they reach a certain size, and a seemingly stable wave seems
to form for a period of time, before it eventually ruptures, as can be seen in
Figs. 3(b) and 3(d). For shear stronger than a critical value of $B\approx
12.5$, this stable translating perturbation does not rupture within the
simulation time. We denote this critical value as $B_{crit}$. For
$B>B_{crit}$, the perturbation waveform seemingly propagates indefinitely as
shown in Figs. 3(c) and 3(d).

In Fig. 3(d) it is apparent that imposing a shear below $B_{crit}$ can delay
rupture to some extent. This effect is quantified in Fig. 4(a), where it can
be seen that imposing $B>5$ leads to a modest increase in rupture time, and
that this effect is dramatically increased as $B$ approaches $B_{crit}$. Fig.
3(d) also suggests that beyond $B_{crit}$, the amplitude of the stable
travelling wave generated under shear is dependent on the size of $B$. This is
indeed the case, as is shown in Fig. 4(b). Just above $B_{crit}$, the minimum
height of the final film profile is around $h_{min}=0.65$, but this increases
as $B$ is increased, eventually converging towards a complete suppression of
the perturbation as B approaches infinity.

Interestingly, we note that the observed value of $B_{crit}\approx 12.5$ is
larger than the value of 9.7 observed for the deterministic case with an
initial perturbation of wavelength $\lambda_{d}$ by Davis et al. [34]. In Fig.
4(c) we plot the increase in rupture time with $B$ for deterministic
simulations (i.e. Eq. (II.1) with $Q=0$ and a sinusoidal initial perturbation)
for a few different perturbation wavelengths. In this plot we note that the
critical shear for the dominant mode $\lambda_{d}$ is indeed around 9.7, but
that larger and smaller wavelengths require more shear in order to suppress
rupture. We propose that fluctuations enable rupture even beyond the critical
shear of the dominant mode by triggering secondary modes once the dominant
mode begins to stabilize. This is supported by the fact that the critical
shear observed in the stochastic case approximately matches the largest values
of $B_{crit}$ for large and small wavelengths, as is shown in Fig. 4(d)

In Fig. 3, we observe that the size of the initial decrease $h_{0}-h_{min}(0)$
in film height due to the stochastic term seems to be of size $Q$ regardless
of the strength of the shear. This is rather unexpected [31] if the size of
the initial perturbation is interpreted as the equilibrium thermal roughness.
We note however, that the sinusoidal perturbations that lead to rupture form
within the first few time steps of the simulation, long before the interface
settles to an equilibrium roughness. This is discussed in more detail in
appendix B, where we also show that the equilibrium thermal roughness is
decreased by $B$ when there is no disjoining pressure.

###  III.2 Three-dimensional case

Figure 5: Contour plots illustrating the height $h(x,y,t)$ of 3D thin films
for different times, with $Q=10^{-3}$. The color represents the deviation
$\delta h(x,y)=h(x,y)-h_{0}$ from the initial film thickness $h_{0}$. $(a-c)$
No-shear ($B=0$) case, where $(c)$ shows the final time step at which the film
ruptures. $(d-f)$ Shear ($B=30$, along $y$-axis) case, where (f) shows the
final time step at which the film ruptures. Note that the direction of shear
(as indicated by the black arrows in panel $(d)$) is not the same as in Fig.
1.

Although our results in two dimensions shed light on 2D film rupture, it is
not clear how these effects translate to physically realistic 3D films. To
quantify the effect of spatial dimensions, we solve Eq. (II.1) to obtain 3D
film profiles, obtaining the results shown in Fig. 5. In Fig. 5$(a-c)$, we
note that the shear-free rupture process is both qualitatively and
quantitatively similar to the observations in two dimensions (see Fig.
2$(a)$): the stochastic fluctuations give rise to small perturbations that
coarsen until a characteristic size given by $\lambda_{d}$ and subsequently
grow over time. The surface depressions that lead to rupture have a circular-
like cross-section, but grow with the same growth rate as in the 2D case (Eq.
(9)). The rupture time is still fairly well described by Eq. (12).

$(a)$

$(b)$

Figure 6: $(a)$ Time evolution of the height perturbation $h_{0}-h_{min}(t)$
for different values of $B$ in 3D simulations, with unidirectional shear and
for $Q=10^{-3}$,$h_{0}=1$. The dotted black line line indicates exponential
behaviour with a growth rate given by Eq. (9). $(b)$ The rupture time as a
function of shear magnitude, where the error bars represent the standard
deviation for a set of 5 simulations for each value of $B$.

In panels $(d-f)$ of Fig. 5, we can see how a strong unidirectional shear in
the $x$-direction affects the three-dimensional film rupture process. In this
case, the shear force has a suppressing effect on perturbations along the
shear direction, but perturbations grow freely in the transverse direction. As
the perturbations grow to a significant size, they begin to align with the
direction of shear, and eventually form a clear pattern of ridges. These
ridges then deepen as the film rupture process proceeds. In Fig. 6$(a)$, we
show how the minimum film height evolves in time for 3D simulations of rupture
with unidirectional shear of varying strength. As in the 2D case, the growth
rate in the exponential regime is close to $\omega_{d}$ predicted by the
linear stability analysis. In contrast to the 2D case of Fig. 3, however,
unidirectional shear does not prevent rupture from occurring in a 3D film.
Nevertheless, for $B>B_{crit}$ there is a small increase in rupture time, as
shown in Fig. 6$(b)$, caused by the fact that the perturbations need to
rearrange themselves when they are suppressed in the shear direction.

Figure 7: Contour plots illustrating the height $h(x,y,t)$ of a 3D thin film,
subjected to a rotating shear force with strength $B=30$ and angular frequency
$\Omega=3\pi/32$ in the direction of the arrow in panel $(a)$, as a function
of time, with $Q=10^{-2}$ and $h_{0}=1$. The color represents the deviation
$\delta h(x,y)=h(x,y)-h_{0}$ from the initial film thickness $h_{0}$.

$(a)$

$(b)$

Figure 8: $(a)$ Time evolution of the height perturbation $h_{0}-h_{min}(t)$
for different values of the shear parameter $B$ and angular frequency
$\Omega=\pi/8$, in 3D simulations with rotating shear for $Q=10^{-2}$. The
dotted black line line indicates exponential behaviour with a growth rate
given by Eq. (9) and an arbitrary amplitude. $(b)$ Rupture time $t_{R}$ as a
function of dimensionless shear stress $B$, for various angular frequencies
$\Omega$.

Since the film rupture time is weakly affected by unidirectional shear at the
interface, we then attempted to impose a rotating shear stress, as was done
previously by Davis et al. [34]. To do so, Eq. (II.1) is again solved
numerically, but the direction of the shear force,
$\overrightarrow{e}_{\tau}$, is now varied sinusoidally in time with an
angular frequency $\Omega$. The resulting rupture process can be seen in Fig.
7, and is further analyzed in Fig. 8. During the linear regime, the
perturbations rotate as the direction of $\overrightarrow{e}_{\tau}$ changes,
but grow much as they do without shear. Once the perturbations become big
enough for nonlinear effects to emerge, the perturbations begin to be
suppressed in the direction of shear, leading to a reorientation as in the
unidirectional case. For very small frequencies, $\Omega\ll 2\pi/t_{R,B=0}$,
the change in direction is too slow to affect the rupture process
significantly. When $\Omega$ is increased, however, the shear has time to
switch directions and suppress the perturbations that have grown in the
initially cross-shear direction. This leads to periodic reorientation of the
perturbation profile with an angular frequency equal to $\Omega$, and an
overall delay of rupture. For rapidly rotating shear, the perturbations remain
essentially circular in cross-section and resemble the perturbations seen in
the shear-free case shown in Fig. 5$(a-c)$. If $B$ is sufficiently large, and
$\Omega$ has an appropriate value, the perturbation can be completely
suppressed until the maximum simulation time of 50, as is the case in Fig. 7.

For a given angular frequency, $\Omega$, the change in rupture time as $B$ is
increased, as depicted in Fig. 8, is actually quite similar to the 2D case of
section III.1. For $B<10$, shear does not significantly delay rupture, but
does so once $B\gtrapprox 15$ due to slower growth in the nonlinear regime.
The critical value of $B$ required is significantly higher than that observed
in 2D, and varies with $\Omega$. As seen in Fig. 8$(b)$, it seems that the
rupture-delaying effect of shear is strongest when $\Omega\approx 3\pi/4$. For
lower frequencies, the shear suppresses rupture in one direction, but does not
change fast enough to suppress it in the other, whereas for high frequencies,
the shear does not remain in one direction for long enough to suppress rupture
in that direction.

###  III.3 Implications for polymer processing

Let us now discuss the implications of the numerical results obtained above in
a practical context. We specifically consider a polymer processing method:
nanolayer coextrusion. This innovative process is based on a series of layer
multiplying elements (LME) which apply the baker’s transformation (successive
slicing and recombining) to a stratified polymer melt flow to achieve
multilayer systems, made of up to thousands of alternating layers of two or
more polymers, each having nanometric thicknesses [42]. Many polymer pairs
have been processed this way, whether glassy or semi-crystalline, miscible or
immiscible and in most cases stable, continuous and regular layers with
thicknesses below 50 nm have been obtained [43, 44, 45]. Let us focus on the
widely studied case of polystyrene (PS) / polymethylmetacrylate (PMMA)
multilayers, for which it has been shown that optimized extrusion conditions
can lead to continuous layers as thin as 20 nm [4]. Typical values for
extrusion temperature, surface tension, dynamic viscosity and Hamaker constant
can be set to 220°C, 1 mN/m [46], 10000 Pa.s [4] and 2.10-18 J [47, 35],
respectively. In this case, the thermal roughness at the interface has an
amplitude $\sqrt{k_{B}T/\gamma}\approx 2.5\cdot 10^{-9}$ m, and $Q\approx
10^{-1}$ which is above the values presented in this study. Nonetheless, from
Figs. 2 and 3, it can be hypothesized that though faster, rupture will occur
similarly as with lower $Q$ values. Using Eq. (13) and putting back
dimensions, the rupture time without shear for a layer having an initial
thickness of 20 nm ($h_{0}\approx 1.1$) would be in the 5 - 10 s range, hence
much less than the typical processing time, between 1 to 2 minutes [48]. The
stabilizing effect of shear could explain the stability of nanometric layers
(around 20 nm), as already discussed in [35]. Looking at the 2D case (Fig.
4(c)), a significant increase of the rupture time occurs when B is higher than
the critical value around 12.5, corresponding to a shear rate in the order of
20 - 25 s-1, that is easily reached in classical extrusion conditions [48].
However, the novel 3D simulations of the present study in uniaxial shear
suggest that the shear force, $B$, has in fact a limited influence on the
rupture time, which appears to contradict the previous conclusions drawn in
2D. The following hypotheses can be made concerning processing, especially
nanolayer coextrusion: the minimal thickness achieved after the LMEs is
actually higher than the one discussed previously, nanometric thicknesses
being only reached after the LMEs when the flow goes through the flat die. As
the rupture time scales with the thickness to the power of 5, rupture may only
occur when the polymer flow passes through this flat die, where the layer
thickness reaches values below 100 nm. Interestingly, in the flat die, there
is not only unidirectional shear but also a diverging radial flow, which may
result in a more stabilizing situation as suggested by results presented in
Figs. 7 and 8. Besides, we stress that the effect of elongation has not been
elucidated yet, and we shall also mention that the boundary conditions used
here in a single layer differ quite substantially with those occurring in a
multilayer flow [49].

Finally, we comment on the 3D surface profiles showing the appearance of
ridge-like and valley-like features aligned along the shear direction prior to
rupture. This pattern formation is reminiscent of the experimental
observations by Dmochowska et al. [36], where a PS thin film sandwiched
between two PMMA thicker layers was dewetted under shear. Though this study
focused on the dewetting kinetics, i.e. the hole growth after hole formation,
it was observed that the dewetting holes take ellipsoidal shapes aligned with
the local shear directions, contrary to the no-shear case [50, 51]. The hole
growth was accelerated in the shear direction while it remained similar to the
no-shear case (with a slight increase caused by shear-thinning of the PMMA
matrix) in the perpendicular direction.

##  IV Conclusion

In this work, we have described how thermal fluctuations and shear affect the
stability of thin nanometric films described by the stochastic thin film
equation. Finite element numerical solutions of the latter show that the role
of the fluctuations is essentially to initiate perturbations of a
characteristic size that is proportional to the thermal roughness of the
interface. At later times, thermal fluctuations do not play a significant role
in the dynamics of thin film rupture. By using a linear stability analysis, we
give a rather simple prediction for the rupture time of a thin film as a
function of initial height and temperature, in the absence of shear. When a
shear force is introduced to a 2D system, rupture can be suppressed, resulting
in the formation of a permanent travelling wave, while the initial size of the
perturbations seems unaffected in the parameter range explored by the
simulations. In the more physically realistic case of unidirectional shear in
3D films, however, rupture is not suppressed, as cross-shear perturbations
grow unimpeded. Our results may explain why the reported rupture-inhibiting
effect of shear has not been reported experimentally. This may be relevant in
physical processes such as the thin air film formed below a droplet as it
impacts on a surface, for which rupture is observed despite a high shear rate
[52]. Our simulations do indicate, however, that rupture can be delayed if the
direction of shear varies sufficiently rapidly with time – a situation of
potential practical relevance for dewetting experiments and nanocoextrusion
processes.

###### Acknowledgements.

The authors acknowledge financial support from the European Union through the
European Research Council under EMetBrown (ERCCoG-101039103) grant, from the
Research Council of Norway through the program NANO2021 (project number
301138), and from the Agence Nationale de la Recherche under EMetBrown
(ANR-21-ERCC-0010-01), Softer (ANR-21-CE06-0029) and Fricolas
(ANR-21-CE06-0039) grants.

##  Appendix A Adaptive time step simulations

As described in section II.2, we stop our simulations when the minimum height
of the film reaches a threshold value of $h^{*}=0.33$, and define the time at
which this occurs as the rupture time of the film, $t_{R}$. Our justification
for this is that the late stage of rupture is an almost instantaneous event
that follows the power law described by Zhang and Lister [14]. In order to
confirm that this is indeed a good assumption, we have performed simulations
as described in section II.2, but with an adaptive time step that decreases as
the film height decreases, as well as a reduced grid spacing of $\Delta
x=0.001$ in order to resolve the flow details in this regime. Fig. 9 shows the
results of such simulations. It is clear here that the late stage of rupture
is accurately described by the aforementioned power law. The time elapsed
between $h_{min}=0.33$ and $h_{min}=0.015$ is less than $0.1$, which is a
negligible fraction of the rupture time $t_{R}$, for all the cases simulated
in this article. The results when shear is included are similar.

$(a)$

$(b)$

$(c)$

Figure 9: Results of the simulations with an adaptive time step near rupture,
for $Q=0.0001$, $B=0$, and $h_{0}=1$. $(a)$ Film profiles at various times.
$(b)$ Minimum film height $h_{min}$ as a function of the reverse time
$t_{R}-t$, in log-log scale. The $(t_{R}-t)^{1/5}$ power law [14] with
arbitrary prefactor is indicated by the dotted black line. $(c)$ Amplitude of
the height perturbation as a function of time. The dotted black line line
indicates exponential behaviour with a growth rate given by Eq. (9) and an
arbitrary amplitude. The dotted red line indicates the $Q$ value.

.

##  Appendix B Thermal roughness of a flat interface with shear

In section III.1, we find that shear does not significantly change the initial
size of the perturbations that lead to thin film rupture when thermal
fluctuations are imposed. This is rather surprising in light of previous works
which show that imposing a shear force significantly decreases the equilibrium
thermal roughness of interfaces in other geometries [31, 30]. We thus simulate
a 2D thin film subjected to thermal fluctuations and shear but no disjoining
pressure, for which the interface fluctuates around its initial height $h_{0}$
indefinitely. We then compute the thermal roughness of the interface at each
time step as $\sigma=(\sum_{n=1}^{N}(h-h_{0})^{2}/N)^{1/2}$, where $n$
represents each of $N$ total gridpoints. In Fig. 10$(a)$ we show how $\sigma$
varies in time when there is no shear. In both cases, $\sigma$ fluctuates
around a mean value of around $Q$ as expected due to the definition of $Q$ in
section II.1.

When we impose a shear $B\neq 0$ to the interface, we do indeed observe a
decrease in the thermal roughness, as is shown in Fig. 10$(b)$ . In fact, for
$B>3$, $\sigma$ is reduced by more than half. Although consistent with what we
expect from the literature, this result seems to conflict with the results in
Fig. 3$(d)$ where the initial perturbation size is independent of $B$. Upon
deeper reflection, however, it is not obvious that the equilibrium thermal
roughness observed over time is the same as the size of instantaneous
fluctuations. In our simulations of film rupture with both fluctuations and
shear, it seems that perturbations to a flat film of size approximately $Q$
are created during the first few timesteps before the disjoining pressure or
shear have any effect. The sinusoidal perturbations that lead to rupture then
form before the shear effect on thermal roughness has had any effect. It is
thus important to distinguish between the time-averaged equilibrium effect of
fluctuations and their instantaneous impact on the film.

$(a)$

$(b)$

Figure 10: Results of the simulations of a flat film without disjoining
pressure with $Q=0.01$ and $h_{0}=1.0$. $(a)$ Instantaneous measured thermal
roughness $\sigma$ as a function of time for the case of $B=0$. The dotted
orange line shows the expected thermal roughness of $\sqrt{2}Q$. $(b)$ The
thermal roughness as a function of the dimensionless shear force $B$. The
black line shows the value of $\sqrt{2}Q$.

## References

  * Wong  _et al._ [1996] H. Wong, I. Fatt, and C. Radke, Deposition and thinning of the human tear film, Journal of Colloid and Interface Science 184, 44 (1996).
  * Chandran Suja  _et al._ [2022] V. Chandran Suja, A. Verma, E. Mossige, K. Cui, V. Xia, Y. Zhang, D. Sinha, S. Joslin, and G. Fuller, Dewetting characteristics of contact lenses coated with wetting agents, Journal of Colloid and Interface Science 614, 24 (2022).
  * Li  _et al._ [2020] Z. Li, A. Olah, and E. Baer, Micro- and nano- layered processing of new polymeric systems, Progress in Polymer Science 102, 101210 (2020).
  * Bironeau  _et al._ [2017] A. Bironeau, T. Salez, G. Miquelard-Garnier, and C. Sollogoub, Existence of a critical layer thickness in ps/pmma nanolayered films, Macromolecules 50, 4064 (2017).
  * Craster and Matar [2009] R. Craster and O. Matar, Dynamics and stability of thin liquid films, Reviews of Modern Physics 81 (2009).
  * Ferrell and Hansford [2007] N. Ferrell and D. Hansford, Fabrication of micro- and nanoscale polymer structures by soft lithography and spin dewetting, Macromolecular Rapid Communications 28, 966 (2007).
  * Dhara  _et al._ [2018] P. Dhara, N. Bhandaru, A. Das, and R. Mukherjee, Transition from spin dewetting to continuous film in spin coating of liquid crystal 5cb, Scientific Reports 8, 7169 (2018).
  * Vrij [1966] A. Vrij, Possible mechanism for the spontaneous rupture of thin, free liquid films, Discuss. Faraday Soc. 42, 23 (1966).
  * Sheludko [1967] A. Sheludko, Thin liquid films, Advances in Colloid and Interface Science 1, 391 (1967).
  * Ruckenstein and Jain [1974] E. Ruckenstein and R. K. Jain, Spontaneous rupture of thin liquid films, J. Chem. Soc., Faraday Trans. 2 70, 132 (1974).
  * Williams and Davis [1982] M. B. Williams and S. H. Davis, Nonlinear theory of film rupture, Journal of Colloid and Interface Science 90, 220 (1982).
  * Oron  _et al._ [1997] A. Oron, S. Davis, and S. Bankoff, Long-scale evolution of thin liquid films, Reviews of Modern Physics 69, 931 (1997).
  * Sharma and Ruckenstein [1986] A. Sharma and E. Ruckenstein, An analytical nonlinear theory of thin film rupture and its application to wetting films, Journal of Colloid and Interface Science 113, 456 (1986).
  * Zhang and Lister [1999] W. W. Zhang and J. R. Lister, Similarity solutions for van der waals rupture of a thin film on a solid substrate, Physics of Fluids 11, 2454 (1999).
  * Dallaston  _et al._ [2018] M. C. Dallaston, M. A. Fontelos, D. Tseluiko, and S. Kalliadasis, Discrete self-similarity in interfacial hydrodynamics and the formation of iterated structures, Phys. Rev. Lett. 120, 034505 (2018).
  * Carlson and Mahadevan [2015] A. Carlson and L. Mahadevan, Protein mediated membrane adhesion, Physics of Fluids 27, 051901 (2015).
  * Becker  _et al._ [2003] J. Becker, G. Grün, R. Seemann, H. Mantz, K. Jacobs, K. R. Mecke, and R. Blossey, Complex dewetting scenarios captured by thin-film models, Nature Materials 2, 59 (2003).
  * Jacobs  _et al._ [2008] K. Jacobs, R. Seemann, and S. Herminghaus, Stability and dewetting of thin liquid films, Polymer Thin Films 1 (2008).
  * Davidovitch  _et al._ [2005] B. Davidovitch, E. Moro, and H. A. Stone, Spreading of viscous fluid drops on a solid substrate assisted by thermal fluctuations, Phys. Rev. Lett. 95, 244505 (2005).
  * Grün  _et al._ [2006] G. Grün, K. Mecke, and M. Rauscher, Thin-film flow influenced by thermal noise, Journal of Statistical Physics 122, 1261 (2006).
  * Carlson [2018] A. Carlson, Fluctuation assisted spreading of a fluid filled elastic blister, Journal of Fluid Mechanics 846, 1076–1087 (2018).
  * Pedersen  _et al._ [2019] C. Pedersen, J. F. Niven, T. Salez, K. Dalnoki-Veress, and A. Carlson, Asymptotic regimes in elastohydrodynamic and stochastic leveling on a viscous film, Physical Review Fluids 4 (2019).
  * Marbach  _et al._ [2018] S. Marbach, D. S. Dean, and L. Bocquet, Transport and dispersion across wiggling nanopores, Nature Physics 14, 1108–1113 (2018).
  * Diez  _et al._ [2016] J. Diez, A. González, and R. Fernandez, Metallic-thin-film instability with spatially correlated thermal noise, Physical Review E 93 (2016).
  * Nesic  _et al._ [2015] S. Nesic, R. Cuerno, E. Moro, and L. Kondic, Fully nonlinear dynamics of stochastic thin-film dewetting, Phys. Rev. E 92, 061002 (2015).
  * Mecke and Rauscher [2005] K. Mecke and M. Rauscher, On thermal fluctuations in thin film flow, Journal of Physics: Condensed Matter 17, S3515 (2005).
  * Durán-Olivencia  _et al._ [2019] M. A. Durán-Olivencia, R. S. Gvalani, S. Kalliadasis, and G. A. Pavliotis, Instability, rupture and fluctuations in thin liquid films: Theory and computations, Journal of Statistical Physics 174, 579 (2019).
  * Shah  _et al._ [2019] M. Shah, V. Steijn, C. Kleijn, and M. Kreutzer, Thermal fluctuations in capillary thinning of thin liquid films, Journal of Fluid Mechanics 876, 1090 (2019).
  * Zhao  _et al._ [2022] C. Zhao, J. Liu, D. A. Lockerby, and J. E. Sprittles, Fluctuation-driven dynamics in nanoscale thin-film flows: Physical insights from numerical investigations, Phys. Rev. Fluids 7, 024203 (2022).
  * Derks  _et al._ [2006] D. Derks, D. G. A. L. Aarts, D. Bonn, H. N. W. Lekkerkerker, and A. Imhof, Suppression of thermally excited capillary waves by shear flow, Phys. Rev. Lett. 97, 038301 (2006).
  * Thiébaud and Bickel [2010] M. Thiébaud and T. Bickel, Nonequilibrium fluctuations of an interface under shear, Physical review. E, Statistical, nonlinear, and soft matter physics 81, 031602 (2010).
  * Bresson  _et al._ [2017] B. Bresson, C. Brun, X. Buet, Y. Chen, M. Ciccotti, J. Gâteau, G. Jasion, M. N. Petrovich, F. Poletti, D. J. Richardson, S. R. Sandoghchi, G. Tessier, B. Tyukodi, and D. Vandembroucq, Anisotropic superattenuation of capillary waves on driven glass interfaces, Phys. Rev. Lett. 119, 235501 (2017).
  * Kalpathy  _et al._ [2010] S. Kalpathy, L. Francis, and S. Kumar, Shear-induced suppression of rupture in two-layer thin liquid films, Journal of colloid and interface science 348, 271 (2010).
  * Davis  _et al._ [2010] M. J. Davis, M. B. Gratton, and S. H. Davis, Suppressing van der waals driven rupture through shear, Journal of Fluid Mechanics 661, 522–539 (2010).
  * Kadri  _et al._ [2021] K. Kadri, J. Peixinho, T. Salez, G. Miquelard-Garnier, and C. Sollogoub, Dewetting of a thin polymer film under shear, Polymer 235, 124283 (2021).
  * Dmochowska  _et al._ [2022] A. Dmochowska, J. Peixinho, C. Sollogoub, and G. Miquelard-Garnier, Dewetting dynamics of sheared thin polymer films: An experimental study, ACS Macro Letters 11, 422 (2022).
  * Batchelor [2000] G. K. Batchelor,  _An Introduction to Fluid Dynamics_, Cambridge Mathematical Library (Cambridge University Press, 2000).
  * Israelachvili [2011] J. N. Israelachvili,  _Intermolecular and Surface Forces_ , third edition ed. (Academic Press, San Diego, 2011) pp. 253–289.
  * Logg  _et al._ [2012] A. Logg, K.-A. Mardal, G. N. Wells, _et al._ ,  _Automated Solution of Differential Equations by the Finite Element Method_, edited by A. Logg, K.-A. Mardal, and G. N. Wells (Springer, 2012).
  * Oliphant [2006] T. Oliphant,  _Guide to NumPy_ (2006).
  * Moreno-Boza  _et al._ [2020] D. Moreno-Boza, A. Martínez-Calvo, and A. Sevilla, Stokes theory of thin-film rupture, Physical Review Fluids 5, 014002 (2020).
  * Gholami  _et al._ [2020] F. Gholami, L. Pakzad, and E. Behzadfar, Morphological, interfacial and rheological properties in multilayer polymers: A review, Polymer 208, 122950 (2020).
  * Rijal  _et al._ [2022] B. Rijal, L. Delbreilh, C. Sollogoub, E. Baer, and A. Saiter-Fourcin, Multiscale analysis of segmental relaxation in PC/PETg multilayers: Evidence of immiscible nanodroplets, Macromolecules 55, 6562 (2022).
  * Lozay  _et al._ [2021] Q. Lozay, Q. Beuguel, N. Follain, L. Lebrun, A. Guinault, G. Miquelard-Garnier, S. Tencé-Girault, C. Sollogoub, E. Dargent, and S. Marais, Structural and barrier properties of compatibilized pe/pa6 multinanolayer films, Membranes 11, 75 (2021).
  * Nassar  _et al._ [2018] S. F. Nassar, S. Domenek, A. Guinault, G. Stoclet, N. Delpouve, and C. Sollogoub, Structural and dynamic heterogeneity in the amorphous phase of poly (l, l-lactide) confined at the nanoscale by the coextrusion process, Macromolecules 51, 128 (2018).
  * Wu [1970] S. Wu, Surface and interfacial tensions of polymer melts. ii. poly (methyl methacrylate), poly (n-butyl methacrylate), and polystyrene, The Journal of Physical Chemistry 74, 632 (1970).
  * de Silva  _et al._ [2012] J. P. de Silva, F. Cousin, A. R. Wildes, M. Geoghegan, and M. Sferrazza, Symmetric and asymmetric instability of buried polymer interfaces, Physical Review E 86, 032801 (2012).
  * Bironeau [2016] A. Bironeau,  _Films multinanocouches de polymères amorphes coextrudés : élaboration, caractérisation et stabilité des nanocouches_ , Ph.D. thesis, ENSAM (2016).
  * Lenz and Kumar [2007] R. D. Lenz and S. Kumar, Instability of confined thin liquid film trilayers, Journal of Colloid and Interface Science 316, 660 (2007).
  * Zhu  _et al._ [2016] Y. Zhu, A. Bironeau, F. Restagno, C. Sollogoub, and G. Miquelard-Garnier, Kinetics of thin polymer film rupture: Model experiments for a better understanding of layer breakups in the multilayer coextrusion process, Polymer 90, 156 (2016).
  * Chebil  _et al._ [2018] M. S. Chebil, J. D. McGraw, T. Salez, C. Sollogoub, and G. Miquelard-Garnier, Influence of outer-layer finite-size effects on the dewetting dynamics of a thin polymer film embedded in an immiscible matrix, Soft matter 14, 6256 (2018).
  * Kolinski  _et al._ [2012] J. M. Kolinski, S. M. Rubinstein, S. Mandre, M. P. Brenner, D. A. Weitz, and L. Mahadevan, Skating on a film of air: Drops impacting on a surface, Phys. Rev. Lett. 108, 074503 (2012).


# A numerical cough machine

Cesar Pairetti  Sorbonne Université and CNRS, Institut Jean Le Rond
d’Alembert, UMR 7190, Paris, France  Centro Internacional Mecánica
Computacional (CONICET - UNL), Santa Fe, Argentina  Facultad de Ciencias
Exactas, Ingeniería y Agrimensura (UNR), Rosario, Argentina  Raphaël Villiers
Sorbonne Université and CNRS, Institut Jean Le Rond d’Alembert, UMR 7190,
Paris, France  Stéphane Zaleski  Sorbonne Université and CNRS, Institut Jean
Le Rond d’Alembert, UMR 7190, Paris, France  Institut Universitaire de France,
IUF, Institut Jean Le Rond d’Alembert, UMR 7190, Paris, France
<EMAIL_ADDRESS>
###### Abstract

We introduce a simplified model of physiological coughing or sneezing, in the
form of a thin liquid layer subject to a rapid (30 m/s) air stream. The setup
is simulated using the Volume-Of-Fluid method with octree mesh adaptation, the
latter allowing grid sizes small enough to capture the Kolmogorov length
scale. The results confirm the trend to an intermediate distribution between a
Log-Normal and a Pareto distribution $P(d)\propto d^{-3.3}$ for the
distribution of droplet sizes in agreement with a previous re-analysis of
experimental results by one of the authors. The mechanism of atomisation does
not differ qualitatively from the multiphase mixing layer experiments and
simulations. No mechanism for a bimodal distribution, also sometimes observed,
is evidenced in these simulations.

Coughing and sneezing are two processes by which a large number of droplets of
muco-salivary fluid are exhaled and subsequently travel large distances in the
environment Bourouiba et al. (2014); Bourouiba (2020a). These phenomena have
acquired an acute interest in the context of the so-called aerosol
transmission route of the Covid-19 pandemic Balachandar et al. (2020), but
have been studied for near a century in the context of respiratory diseases in
general Bourouiba (2020b). Early investigations by Duguid Duguid (1946) and
Loudon & Roberts Loudon and Roberts (1967) of the number and size of droplets
emitted in coughing and sneezing events have yielded rich data, incorporating
very large numbers of droplets. A comprehensive analysis of the
characteristics of these droplets, both in size and velocity, would be of
immense interest, as it is a prerequisite for the modeling of the droplet
cloud propagating further downstream from the mouth. Rapid photographic
imaging Scharfman et al. (2016) has revealed features similar to those
observed in other atomization processes, including thin liquid sheets,
ligaments and droplets. In this context, the statistical distribution of
droplet sizes is of particular interest. Although the log-normal distribution
has been frequently mentionned in connection with exhalations Nicas et al.
(2017); Wells et al. (1955) as well as other atomizing flows Ling et al.
(2017) many other distributions $N(d)$ of the diameter $d$ have been put
forward such as compound gamma distributions Villermaux and Bossa (2011) and
many others. A re-analysis of the data of Duguid and Loudon & Roberts for
sneezes has however revealed Balachandar et al. (2020) a $N(d)\sim d^{-2}$
scaling over an impressive three orders of magnitude. This scaling allows
determining the proportion of millimeter-sized droplets that travel short
distances and the proportion of much smaller droplets that can be incorporated
in a turbulent puff or particle-laden jet and travel long distances as
discussed in ref. Bourouiba et al. (2014); Balachandar et al. (2020). The
fraction of the exhaled muco-salivary fluid in each class of droplet sizes may
also be determined in this way and the probability of having viral loads in
each class can be inferred under adequate hypotheses, such as a homogeneous
distribution of the virus in volume or surface.

In order to better understand the fluid mechanics of exhalations, King, Brock
& Lundell King et al. (1985) have designed a physical model of violent
exhalation that may be nicknamed a “cough machine”. Air is flowing at high
speed (from 10 to 30 m/s) in a rectangular-section duct, with a flow rate
analog of the observed human cough. A thin layer of muco-salivary fluid is
deposited at the bottom of the duct. While King, Brock & Lundell use the cough
machine to observe non-atomizing waves on the thin layer, it can also be used
to simulate droplet formation at higher speeds and/or lower velocities. The
device would then achieve atomization through a process similar to that of
shear flow atomization of planar sheets Gorokhovski and Herrmann (2008);
Villermaux and Bossa (2011); Agbaglah et al. (2017); Ling et al. (2017);
Zandian et al. (2019); Ling et al. (2019). (See also the recent review of
numerical approaches in Constante-Amores et al. (2020).) However there are
important differences since the liquid layer is initially at rest and the
airflow is impulsive. These differences could lead to different distributions
of droplet sizes and velocities, which are the focus of the current
investigation.

We model the muco-salivary fluid and the surrounding gas as a Newtonian fluid
(The muco salivary fluid is non-newtonian but these effects only kick-in at
very small scales). Thus the flow is described by the Navier-Stokes equations,
which we solve by a simple finite volume discretization in the one-fluid
numerical approach to two-phase flowTryggvason et al. (2011), using a Volume-
Of-Fluid method for tracking fluid interfacesTryggvason et al. (2011), an
octree grid and the Bell-Collela-Glaz advection scheme on staggered grids as
described in Popinet (2003), a Crank-Nicholson method for viscosity Lagrée et
al. (2011), and a height-function method for curvature and surface tension
Popinet (2009). The octree grid is refined using a wavelet estimate of the
local error, as described in the self-documented code Basilisk
(http://basilisk.fr). We use thirteen levels of refinement at the maximum,
which results in the equivalent of $2^{39}\simeq 10^{13}$ or sixteen billion
grid cells. However, thanks to the octree nature of the simulation, only
10-100 million grid cells are used in practice.

Figure 1: A view of the air-liquid interface at t=9 ms. The interface is
colored by the velocity of the droplets, with green the smallest velocity and
dark red the largest.

The physical domain subject to the simulation is a cube of dimensions $L^{3}$
with flow inlet at $x=0$ and flow outlet at $x=L$. The flow is channeled
through a tube of length $\ell_{x}$, and rectangular cross section with height
$\ell_{y}$ and width $\ell_{z}$, with walls or plates of thickness $e$
centered on the planes $x=L/2-\ell_{z}/2$ (bottom plate), $y=\pm\ell_{y}/2$
(lateral plate). The bottom plate is covered with a liquid phase of thickness
$h$, density $\rho_{l}$ and viscosity $\mu_{l}$. The gas phase, modelling the
exhaled air, has density $\rho_{a}$, and viscosity $\mu_{a}$. It enters the
flow from the left wall at $x=0$ with uniform and constant velocity $U$
through an inflow boundary condition in the region
$L/2-\ell_{y}/2<y<L/2+\ell_{y}/2$ and $L/2-\ell_{z}/2+h<z<L+\ell_{z}/2$. This
condition ensures the inflow of air blows just above the initial position of
the liquid layer. The surface tension of the liquid is noted $\sigma$. Gravity
is neglected as it is small compared to inertial effects $g\ell_{z}/U^{2}\ll
1$.

$L$ | $\ell_{x}$ | $\ell_{y}$ | $\ell_{z}$ | $e$ | $h$ | $\rho_{l}$ | $\mu_{l}$ | $\rho_{a}$ | $\mu_{a}$ | $U$ | $\sigma$
---|---|---|---|---|---|---|---|---|---|---|---
0.15 | 0.05 | 0.01 | 0.02 | 0.002 | 0.001 | 1000 | 0.005 | 1.2 | $2\,10^{-5}$ | 30 | 0.03
Table 1: Dimensional values of the fluid and geometrical parameters

The main dimensionless parameters are the Reynolds number of the air based on
the channel height ${\rm Re}\,_{a}=\rho_{a}U\ell_{z}/\mu_{a}=36,000$, the
Weber number of the air based on the channel height ${\rm
We}\,_{a}=\rho_{a}U^{2}\ell_{z}/\sigma=720$, and the Reynolds number of the
liquid based on the height of the liquid ${\rm
Re}\,_{l}=\rho_{l}Uh/\mu_{l}=30,000$. The values of these parameters are quite
high which makes it surprising that the simulations, performed with a
Navier–Stokes code, are still classified as Direct Numerical Simulations
(DNS). Indeed the assessment of the DNS character of the simulation may be
inferred from the value of the dissipation $\epsilon$ in the similar setup of
refs. Ling et al. (2017, 2019). There the kinetic energy dissipation per unit
volume $\epsilon$ Pope (2001) was measured and it was observed that the
maximum value of $\epsilon/(\rho_{a}U^{3}/\ell_{z})$ was about 0.01 which
yields an estimate of $\epsilon$. The Kolmogorov length scale is then
$\eta=(\rho_{a}\nu_{a}^{3}/\epsilon)^{1/4}$ With the value of $\epsilon$
estimated above, we have $\eta_{1}\simeq 24.2$ microns while the size of the
smallest grid cell with the thirteen levels of refinement in the simulation
reported here is $\Delta=2^{-13}\,L\simeq 18.3$ microns. According to the DNS
resolution criterion given by Pope Pope (2001), the smallest turbulent scales
will be well resolved if $\Delta/\eta\leq 2.1$ while we have
$\Delta/\eta_{1}\simeq 0.75$ Even if one uses the more conservative estimate
$\eta_{2}=\ell_{z}{\rm Re}\,_{a}^{-3/4}$ as in Constante-Amores et al. (2020)
one gets $\Delta/\eta_{2}\simeq 2.4$. The perhaps surprising result that the
simulations may be qualified as DNS can be explained by the fact that while
the less extreme simulations of Ling et al. (2017, 2019) were limited to
$\ell_{z}/\Delta=256$ grid points in the gas jet thickness, our simulations
using octree refinement go up to the equivalent of $\ell_{z}/\Delta=1092$.

Simulations are initialized with zero velocity in the air and liquid although
the incompressibility condition results in a non-zero velocity field
everywhere in the gas immediately after time zero. The simulation is continued
for a total time T=9 ms. The liquid surface is quickly significantly perturbed
with waves present over the entire length $\ell_{x}$ of the tube, with a much
larger wave near the inlet and some secondary waves near the outlet (Fig. 1).
The wave stretches into a thin liquid layer that fills with air under the
combined effect of pressure and vorticity. Eventually, two mechanisms lead to
the formation of droplets: the formation of fingers or ligaments at the end of
the sheets and the puncturing of the sheets. The latter causes the formation
of holes in the sheets and the subsequent expansion of those holes, leading to
the formation of ligaments that eventually pinch and break into droplets. This
process has been described in other experimental Scharfman et al. (2016) and
numerical Ling et al. (2017) investigations. In the current “closed channel”
configuration similar mechanisms of droplet formation are observed. It is seen
on Fig. 1 that the small dark red droplets near the channel exit have moved
quickly since the drag law scaling for a sphere
$F_{d}\sim\rho_{a}d^{2}U^{2}/4$ (where $d$ is the droplet radius) implies that
the small droplets accelerate much faster than the large ones.

Figure 2: In a cross section of the channel, the liquid phase and the solid
are shown in black and the gas phase in white. Three different snapshots at
regularly spaced time intervals are shown. The oblique line connects the
positions, thus giving a graphical display of the wave velocity $U_{D}$,
agreing with the Dimotakis velocity discussed in the text.

A clear phenomenon that couples with droplet formation is growth of a large
sheet or wave near the inlet and its progression downstream. The velocity of
such waves is typically the Dimotakis velocity Dimotakis (1986) and it is
expressed as $U_{D}=U{\sqrt{\rho_{a}}}/{(\sqrt{\rho_{a}}+\sqrt{\rho_{l}})}$. A
similar kind of solitary wave progressing at the Dimotakis velocity has been
observed in the simpler setting of an infinite vortex sheet between air and
liquid Hoepffner et al. (2011). The liquid surface is excited by a local
perturbation of the flat vortex sheet, and in this case the inhomogeneity of
the flow at the entrance plays the role of the localised perturbation, while
less localised waves are seen further downstream. In our case the Dimotakis
velocity is $U_{D}\sim 1.004$ m/s which may be compared to a rough measurement
from the simulations (Fig. 2) of $U_{D,\rm{num}}\simeq 1.25$ m/s. The
identification of droplets or connected fluid components, and the computation
of their volumes $V_{d}$ allows to define an equivalent droplet diameter
$d=(6V_{d}/\pi)^{1/3}$.

Figure 3: Droplet counts in each bin, proportional to the Probability
Distribution Function of droplet diameters, for $t=9$ ms and for two grid
resolutions L13: $\Delta x=18\mu$m, L12: $\Delta x=36\mu$m and L11: $\Delta
x=72\mu$m.

Because of the presence of some non-spherical droplets in the sample, we
compared the filtered counts, incorporating only near-spherical droplets using
a sphericity index. We find that this filtering removes droplets with $d>2$mm.
(There are less than ten such droplets.) The resulting number frequency is
shown in Fig. 3 together with the statistical error bars, the error being
defined as one standard deviation of the binomial. It is seen that the
distribution is close to the power law $N(d)\sim d^{-3.32}$ for small sizes
and then inches down. On Fig. 3 we do not show droplet counts for $d<100$
microns and $d>400$ microns, since the small sizes are plagued by grid
resolution errors and the larger sizes by statistical errors.

An often considered number frequency (NF) model is the log-normal, which reads
$N(d)=({B}/{d})\exp\left[-{(\ln
d-\hat{\mu})^{2}}/{(2\hat{\sigma}^{2})}\right]$ where $B$ is a normalization
constant, $\hat{\mu}$ is the expected value of $\ln d$, also called the
geometric mean, and $\sigma$ is the standard deviation of $\ln d$, also called
the geometric standard deviation (GSD, see Nicas et al. (2017)). If we plot
$y=dP(d)$ versus $x=\ln d$ the Log-Normal frequency distribution appears as a
parabola. This is done in Fig. 4.

Figure 4: Droplet counts in “log-normal coordinates” in which the log-normal
NF appears as a parabola, with labels L11-13 as in Fig. 3.

That figure shows a fit with a GSD $\hat{\sigma}=0.98$. The dimensionless GSD
is of the same order of magnitude as the GSD ($\hat{\sigma}\sim 1$) in similar
experiments and simulations Herrmann (2011); Ling et al. (2017); Marty (2015),
but in the authors’ opinion, this may not capture universal underlying physics
but the similar range of scales that is within numerical or experimental reach
in the literature cited. Comparing the NF at various resolution in Herrmann
(2011); Ling et al. (2017) it is seen that as the grid size is reduced the NF
shifts to the left, with its geometric mean decreasing and its GSD increasing.
This also seen partially in Fig. 3 where the L11 resolution peaks “before” the
two others. It is thus possible that as resolution is increased the curved NF
seen on Figs 3 and 4 would progressively asymptote to a power law, that is a
Pareto, instead of a converged Log-Normal. Such a Pareto NF necessarily has a
lower bound, if only the molecular size. This lower bound is not attainable
numerically or perhaps even experimentally. A tempting hypothesis is to
associate it to the thickness of the sheets as they break or perforate, which
involves mechanisms that are not modelled in this study.

To conclude, we have demonstrated a simple physical analog of the
physiological mechanism of coughing and sneezing, that is strongly reminiscent
of planar sheet atomization processes. A Pareto $d^{-2}$ distribution was not
found but is not excluded at very small diameters. Perspectives include higher
resolutions simulations and laboratory experiments in the regime of this
numerical experiment, and using the characteristics of the numerically
estimated droplet sizes and velocity to predict the further evolution of the
droplet cloud using Lagrangian particle methods such as those of Chong et al.
Chong et al. (2020).

We acknowledge the funding of the Bec.ar programme, the ERC-ADV grant TRUFLOW,
the ANR NANODROP grant funded by the Fondation de France and the PRACE Covid
grant of computer time. We thank Lydia Bourouiba, Pallav Kant, Detlef Lohse
and Stéphane Popinet for stimulating and fruitful discussions, and for sharing
useful material, visualisations and codes.

## References

  * Bourouiba et al. (2014) L. Bourouiba, E. Dehandschoewercker, and J. W. Bush, J. Fluid Mech. 745, 537 (2014).
  * Bourouiba (2020a) L. Bourouiba, JAMA 323, 1837 (2020a).
  * Balachandar et al. (2020) S. Balachandar, S. Zaleski, A. Soldati, G. Ahmadi, and L. Bourouiba, International Journal of Multiphase Flow 132, 103439 (2020).
  * Bourouiba (2020b) L. Bourouiba, Annual Review of Fluid Mechanics 53 (2020b).
  * Duguid (1946) J. P. Duguid, Epidemiology & Infection 44, 471 (1946).
  * Loudon and Roberts (1967) R. Loudon and R. Roberts, Nature 213, 95 (1967).
  * Scharfman et al. (2016) B. Scharfman, A. Techet, J. Bush, and L. Bourouiba, Experiments in Fluids 57, 24 (2016).
  * Nicas et al. (2017) M. Nicas, W. W. Nazaroff, and A. Hubbard, Journal of Occupational and Environmental Hygiene 2, 143 (2017).
  * Wells et al. (1955) W. F. Wells et al., _Airborne Contagion and Air Hygiene. An Ecological Study of Droplet Infections._ (Cambridge: Harvard University Press (for The Commonwealth Fund), Mass., USA., 1955).
  * Ling et al. (2017) Y. Ling, D. Fuster, S. Zaleski, and G. Tryggvason, Phys. Rev. Fluids 2, 014005 (2017).
  * Villermaux and Bossa (2011) E. Villermaux and B. Bossa, J. Fluid Mech. 668, 412 (2011).
  * King et al. (1985) M. King, G. Brock, and C. Lundell, physiology.org 58, 1776 (1985).
  * Gorokhovski and Herrmann (2008) M. Gorokhovski and M. Herrmann, Annu. Rev. Fluid Mech. 40, 343 (2008).
  * Agbaglah et al. (2017) G. Agbaglah, R. Chiodi, and O. Desjardins, J. Fluid Mech. 812, 1024 (2017).
  * Zandian et al. (2019) A. Zandian, W. A. Sirignano, and F. Hussain, International Journal of Multiphase Flow 113, 117 (2019).
  * Ling et al. (2019) Y. Ling, D. Fuster, G. Tryggvason, and S. Zaleski, J. Fluid Mech. 859, 268 (2019).
  * Constante-Amores et al. (2020) C. R. Constante-Amores, L. Kahouadji, A. Batchvarov, S. Shin, J. Chergui, D. Juric, and O. K. Matar, arXiv preprint arXiv:2012.01887 (2020).
  * Tryggvason et al. (2011) G. Tryggvason, R. Scardovelli, and S. Zaleski, _Direct Numerical Simulations of Gas-Liquid Multiphase Flows_ (Cambridge University Press, 2011).
  * Popinet (2003) S. Popinet, J. Comput. Phys. 190, 572 (2003).
  * Lagrée et al. (2011) P.-Y. Lagrée, L. Staron, and S. Popinet, J. Fluid Mech. 686, 378 (2011).
  * Popinet (2009) S. Popinet, J. Comput. Phys. 228, 5838 (2009).
  * Pope (2001) S. B. Pope, _Turbulent flows_ (2001).
  * Dimotakis (1986) P. Dimotakis, AIAA J 24, 1791 (1986).
  * Hoepffner et al. (2011) J. Hoepffner, R. Blumenthal, and S. Zaleski, Phys. Rev. Lett 106, 104502 (2011).
  * Herrmann (2011) M. Herrmann, AAS 21, 283 (2011).
  * Marty (2015) S. Marty, Ph.D. thesis, Université Grenoble Alpes, Université de Grenoble (2015).
  * Chong et al. (2020) K. L. Chong, C. S. Ng, N. Hori, R. Yang, R. Verzicco, and D. Lohse, arXiv preprint arXiv:2008.01841 (2020).