Simon Schmidt, Outi Tammisola, Lutz Lesshafft, Kilian Oberleithner
Global stability and nonlinear dynamics of wake
flows with a two-fluid interface
Open Access via institutional repository of Technische Universität Berlin
Document type
Journal article | Accepted version
(i. e. final author-created version that incorporates referee comments and is the version accepted for
publication; also known as: Author’s Accepted Manuscript (AAM), Final Draft, Postprint)
This version is available at
https://doi.org/10.14279/depositonce-16265
Citation details
Schmidt, S., Tammisola, O., Lesshafft, L., & Oberleithner, K. (2021). Global stability and nonlinear dynamics of
wake flows with a two-fluid interface. In Journal of Fluid Mechanics (Vol. 915). Cambridge University Press
(CUP). https://doi.org/10.1017/jfm.2021.150.
Terms of use
This work is protected by copyright and/or related rights. You are free to use this work in any way permitted by
the copyright and related rights legislation that applies to your usage. For other uses, you must obtain
permission from the rights-holder(s).
1
Global stability and nonlinear dynamics of
wake flows with a two-fluid interface
Simon Schmidt1†, Outi Tammisola2, Lutz Lesshafft3and Kilian
Oberleithner1
1Laboratory for Flow Instability and Dynamics, Technische Universit¨at Berlin, 10623 Berlin,
Germany
2FLOW, Engineering Mechanics, KTH Royal Institute of Technology, SE-100 44 Stockholm,
Sweden
3LadHyX, CNRS/´
Ecole Polytechnique/Institut Polytechnique de Paris, 91128 Palaiseau,
France
(Received xx; revised xx; accepted xx)
A framework for the computation of linear global modes, based on time-stepping of
a linearised Navier-Stokes solver with an Eulerian interface representation, is presented.
The method is derived by linearising the nonlinear solver Basilisk, capable of computing
immiscible two-phase flows, and offers several advantages over previous, matrix-based,
multi-domain approaches to linear global stability analysis of interfacial flows. Using
our linear solver, we revisit the study of Tammisola et al. (2012), who found a counter-
intuitive, destabilising effect of surface tension in planar wakes. Since their original study
does not provide any validation, we further compute nonlinear results for the studied
flows. We show that a surface-tension-induced destabilisation of plane wakes is observable
which leads to periodic, quasiperiodic or chaotic oscillations depending on the Weber
number of the flow. The predicted frequencies of the linear global modes, computed in
the present study, are in good agreement with the nonlinear results, and the growth
rates are comparable to the disturbance growth in the nonlinear flow before saturation.
The bifurcation points of the nonlinear flow are captured accurately by the linear solver
and the present results are as well in correspondence with the study of Tammisola et al.
(2012).
Key words:
1. Introduction
Shear flows involving two (or more) immiscible fluids can exhibit distinctively different
stability characteristics than their single-phase equivalents. This can be due to the
presence of not only a momentum gradient of the adjacent fluids but also a density and
viscosity gradient. Furthermore, the presence of a surface tension force at the interface(s)
can influence the stability of the flow.
As a consequence, stability properties of multiphase flows are significantly more com-
plex and involve a larger parameter space than their single phase equivalents. The Weber
number can be used to express the ratio of momentum to surface tension force in the
flow. In general, for liquid jets with low Weber numbers, the influence of surface tension
†Email address for correspondence: s.sc[email protected]
2S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
is dominating and, depending on whether the liquid layer is planar or round, can either
stabilize (Squire 1953; Hagerty & Shea 1955) or destabilize the flow (Rayleigh 1878). At
higher Weber numbers, aerodynamic forces, through the momentum or density/viscosity
gradient (Yih 1967; Drazin & Reid 2004; Boeck & Zaleski 2005), overcome the effects
of surface tension. However, in the case of confined plane jets and wakes, a study by
Tammisola et al. (2012) found that for certain configurations, surface tension might as
well destabilise a plane flow which is stable in its absence. These somewhat surprising
results are partially confirmed in a subsequent study by Biancofiore et al. (2014) who
found an unstable global mode for some of the configurations of the former study.
However, for other configurations they found no global instability, thereby raising doubts
concerning the validity of the surface tension-induced destabilisation found by Tammisola
et al. (2012).
Flow stability is usually assessed through the framework of linear stability analysis,
which seeks solutions of the Navier-Stokes operator, linearised around steady or time-
periodic base flows (see e.g. Schmid & Henningson 2012). By choosing a Fourier ansatz
for the perturbation quantities, the linearised system is recast as an eigenvalue problem,
where the resulting eigenvalues yield information about the exponential growth or decay
of the respective eigenvectors. In a local analysis, the underlying base flow is assumed
to be parallel and therefore only dependent on one spatial coordinate, while for a global
analysis it might be inhomogeneous in all spatial coordinates.
However, for large-scale flows or three-dimensional base flows, memory requirements
for storing the operator matrix are still prohibitive. In such cases direct construction of
the Jacobian matrix can be avoided by employing iterative techniques where the high-
dimensional system is orthogonally projected onto a low-dimensional subspace which in
turn allows for direct computation of its eigenvalues. In practice, construction of the
subspace can be achieved using standard time-stepping techniques and in principle any
nonlinear DNS solver can be utilized, as demonstrated by Barkley et al. (2008).
Although global analysis has become a standard tool for analysing single-phase flows,
its application to multiphase flow has only rarely been attempted so far. There have
been numerous approaches through local absolute stability analysis (e.g. S¨oderberg 2003;
Rees & Juniper 2009; Sevilla 2011), and by analysing the growth of small perturbations
around base flows in direct numerical simulations (Tammisola et al. 2017). However, to
our knowledge, the work of Tammisola et al. (2011, 2012) is one of the very few reported
applications of linear global analysis to a two-dimensional immiscible two-phase flow. In
their approach, the linearised operator is constructed explicitly, and discretisation of the
fluid phases is done on separate meshes, resulting in a conformal interface representation
that is aligned with the mesh boundary, separating the domains. Along the interface,
velocity and stress conditions are enforced to couple the domains. The base flow is
computed using a single-phase spectral element code, and the fluid interface is extracted
a posteriori as a streamline to obtain an artificial two-phase base flow.
In the present work, we explore a different approach, by developing a framework which
allows for computation of global modes by means of time-stepping with a linearised
DNS solver with an Eulerian interface representation, capable of computing two-phase
flows. The benefits of a successful realisation of this approach are readily seen: base
flow and perturbation computations are obtained using the same numerical schemes, so
that no grid mapping or interpolation is necessary. Further, resource requirements for
perturbation computations scale similarly to nonlinear simulations, such that analysis of
three-dimensional base flows is possible. Finally, the utilisation of an available, highly
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 3
optimised, nonlinear solver obviates the need for re-implementing the required schemes
in a new solver.
For our study, we choose the open-source framework Basilisk, developed by St´ephane
Popinet (http://basilisk.fr), which offers a geometrical Volume-Of-Fluid (VOF) interface
representation, combined with a well-balanced surface tension scheme.
The aim of this article is two-fold: First, we give a detailed presentation of the
derivation of the novel method for the computation of linear global modes of two-phase
flows. Second, we revisit the wake flows investigated by Tammisola et al. (2012), this
time using nonlinear simulations and our linear solver, thereby shedding more light on
the nonlinear dynamics of the flow and providing a rigorous validation.
In the remainder of this article, we first give an overview of the schemes implemented in
Basilisk and outline the necessary modifications to the solver for computing linear global
modes. We then proceed by presenting nonlinear results of the wake flow configurations
of Tammisola et al. (2012), and compare them with the linear results obtained with
the present method. Finally, we discuss some of the differences between the studies of
Tammisola et al. (2012) and Biancofiore et al. (2014).
2. Numerical methods for two-phase flows
2.1. Governing equations for interfacial two-phase flow
As the stability computations are facilitated using the same numerical methods as
the nonlinear solver, a recapitulation of the discretisation and schemes used in this
work is given here. From a physical viewpoint, it is assumed that both fluid phases are
separated by an interface of negligible thickness. The molecular imbalance of cohesive
forces between both fluids is modelled as a surface tension force, acting on the interface.
The numerical representation of the phases and the interface can be either Lagrangian
or Eulerian. Current nonlinear solvers usually use the level set method (Sussman et al.
1994), the Volume-Of-Fluid method (VOF) (e.g. Scardovelli & Zaleski 1999) or methods
derived thereof, all of which use an Eulerian representation, resulting in a non-conformal
interface representation of finite thickness. The incompressible continuity and momentum
equations, including variable density and surface tension, are given in unified form as
∂tρ+u· ∇ρ= 0,(2.1a)
ρ(∂tu+u· ∇u) = −∇p+∇ · (µD) + σκnδ(x−xs) (2.1b)
with u= (u, v, w)Tthe velocity vector, ρthe density, µthe dynamic viscosity, p
the pressure and xsbeing the position of the interface. The deformation tensor is D=
∇u+∇Tu. Density and viscosity differ between the phases but are constant within each
phase. The rightmost term in equation (2.1b) represents the surface tension force along
the interface, and is composed of the surface tension coefficient σ, the local interface
curvature κ, the unit normal vector of the interface nand δ, the Dirac δ-function that
is non-zero on the interface. Using a Heaviside function H(x−xs), that is 1 in phase 1
and 0 in phase 2, ρand µcan be expressed as
ρ=ρ2+H(ρ1−ρ2),(2.2a)
µ=µ2+H(µ1−µ2).(2.2b)
4S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
2.2. Interface representations
Numerically, Hand δare approximated as Hand δwhere is a characteristic length
scale related to the local grid size ∆. The method of computing the surface tension term
is closely related to the immersed boundary method, introduced by Peskin (1972) where
we find a volumetric representation of the surface force as
σκnδ(x−xs) = σκ∇H(x−xs).(2.3)
While several numerical representation of Hare possible, we will focus on two variants
of the continuum surface force method (CSF) (Brackbill et al. 1992). In the CSF method,
combined with the Volume of Fluid method we set
H(x−xs) = c(x) =
0,for xin the phase 1,
1,for xin the phase 2,
0.5,for xat the interface,
(2.4)
where cis the volume fraction field and =∆, the grid size. We then find nδ=∇c.
Similarly, when combining the CSF method with a level set method, we can find a smooth
approximation
H(x−xs) = H(φ(x)) =
0,if φ < −,
1,if φ > ,
1+φ/+sin(πφ/)/π
2,otherwise.
(2.5)
where φis usually chosen as a signed distance function with respect to the interface:
φ(x) =
φ > 0,for xin the phase 1,
φ < 0,for xin the phase 2,
φ= 0,for xat the interface.
(2.6)
Here, we find nδ=∇φ/|∇φ|δwhere the smooth Delta function can be obtained as
δ= dH(φ)/dφ. Both methods usually lead to a characteristic interface thickness of
O(∆) (Popinet 2018). Since ρis directly coupled to cor φ, respectively, equation (2.1a)
is equivalent to the advection of cor φ:
∂tc+∇ · (cu) = 0 (2.7a)
∂tφ+∇ · (φu)=0.(2.7b)
In Basilisk, the CSF method in combination with a VOF interface representation is
used. However, for reasons described in §3.3, we will adopt aspects of the level set method
for the derivation and solution of the linearised code.
The discretisation of (2.3) warrants special attention as well-balanced schemes have to
be used which are able to recover equilibrium solutions of certain continuous problems,
thus avoiding the problem of so-called parasitic currents (Harvie et al. 2006). An extensive
discussion of this matter is given in Popinet (2018). An essential requirement is that the
same discrete operators are used for the gradients of the pressure pand the Heaviside
function H. The computation of the normals and curvature depends on the method used
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 5
for the interface representation. Since for the level set method, φis a continuous function,
the curvature is easily computed as
n=∇φ
|∇φ|,(2.8a)
κ=∇ · n.(2.8b)
However, for the VOF-CSF method, direct derivation of cleads to problematic cur-
vature estimates because of its discontinuity. Therefore Basilisk, uses a height-function
method that gives second-order accurate curvature estimates (Cummins et al. 2005).
For illustration, in two dimensions a close-to-horizontal interface can be described by
y=hy(x), where hyis the vertical distance to the interface at a given x. The normals
and curvature of this interface are then computed as
n=(h0
y,1)
(1 + h2
y)1/2(2.9a)
κ=h00
y
(1 + h02
y)3/2.(2.9b)
For an interface x=hx(y) the procedure is similar. An extension to three dimensions
is straight forward and described in Popinet (2009).
2.3. Discretisation in Basilisk
The equation system (2.1) is discretised on regular Cartesian grids using a staggered-
in-time discretisation, which is second-order accurate, yielding the following velocity and
volume fraction update at every time step Popinet (2003, 2009).
un+1 =un+∆t −un+1/2· ∇un+1/2− ∇pn+1/ρn+1/2
+∇ · [2(µ/ρ)n+1/2(Dn+Dn+1)] + [σ/ρκ∇c]n+1/2,(2.10a)
cn+1/2=cn−1/2−∆t ∇ · (cnun).(2.10b)
Pressure and velocity are decoupled using a standard time-splitting projection scheme
(Chorin 1969). The advection term is computed using the Bell-Colella-Glaz (BCG)
second-order unsplit upwind scheme (Bell et al. 1989), as will be described in §3.2. A
geometric VOF method is used to advect the volume fraction in equation (2.7a). The
advection of a level set function φis equivalent to the advection of a passive scalar which
is advected using
φn+1/2=φn−1/2−∆t ∇ · (φnun).(2.11)
3. Linearisation procedure
3.1. Derivation of the linearised equations
In order to obtain the linearised equations, we first non-dimensionalise the nonlinear
equations (2.1b) with respect to ρ1,µ1and a suitable reference length and velocity scale,
thus obtaining
6S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
[˜ρ+H(φ)(1 −˜ρ)](∂tu+u· ∇u) =
− ∇p+1
Re∇ · (˜µ+H(φ)(1 −˜µ))(∇u+∇Tu)+1
Weκnδ(φ),(3.1a)
where we have used the level set methodology to account for the two phases. The
Reynolds number is Re =ρ1Uref Dref /µ1, the Weber number is We =ρ1U2
ref Dref /σ,
˜ρ=ρ2/ρ1and ˜µ=µ2/µ1. Note that, from equation (3.1) onward, all quantities are
assumed to be dimensionless.
Starting point for the assessment of linear stability is the decomposition of the flow
into a basic state and an infinitesimal perturbation, u=U+ζu0and p=P+ζp0for the
velocity and pressure respectively, with ζ1. A similar linearisation is done for the level
set function φ=Φ+ζφ0. Therewith associated are a perturbed curvature κ=K+ζκ0
and a normal vector n=N+ζn0. We also make use of the base flow volume fraction C
but, as will be described in §3.3, the use of a perturbation volume fraction poses several
challenges and is thus avoided.
We insert the respective decompositions into equations (3.1, 2.7b) and retain only
leading-order terms in ζ. Dropping ζfor convenience, we arrive at
[˜ρ+H(Φ)(1 −˜ρ)](∂tu0+u0· ∇U+U· ∇u0) + 2[δ(Φ)φ0(1 −˜ρ)](U· ∇U)
=−∇p0+1
Re∇ · (˜µ+H(Φ)(1 −˜µ)) ∇u0+∇Tu0
+1
Re∇ · (δ(Φ)φ0(1 −˜µ)∇U+∇TU+Fs(Φ, φ0),(3.2a)
∂tφ0+∇ · (Φu0) + ∇ · (φ0U)=0,(3.2b)
where we have used the fact that δ(Φ) is obtained as the distributional derivative of
H(Φ), during linearisation of terms involving H(φ). We further note that, as described in
§2.2, H(Φ) can be replaced by C. Thereby, we retain an interface thickness of O(∆). Due
to the involvement of two immiscible phases, the linearised equations contain a number
of additional terms, compared to the linearised incompressible equations of a single fluid
phase: The transport of the perturbation velocity by the base flow velocity and vice versa
is scaled by the base flow density field. Additionally, a new advective term, composed only
of the base flow velocity, appears that is acted upon by the perturbation density field.
When comparing to a usual multi-domain approach of Navier-Stokes or Orr-Sommerfeld
equations, where the interface height appears as a variable, this term corresponds to
multiplication of the same base flow terms with a perturbation of the interface height.
Similarly, a scaling of the perturbation velocity diffusion by the base flow viscosity field
is introduced and a new term, representing the action of the perturbation viscosity field
on the base flow velocity diffusion, appears.
The linearised surface tension force is given as
Fs(Φ, φ0) = 1
We [κ0Nδ(Φ) + Kn0δ(Φ) + KNφ0∂Φδ(Φ)] ,(3.3)
where ∂Φin the rightmost term denotes the functional derivative with respect to Φ, and
f(∂Φδ(Φ)) = (∂Φf)δ(Φ) for some test function f. As long as Φhas signed distance, we
can make use of the fact that, at the interface, any change of Φis along the normal N,
such that ∂NΦ= 1. Hence, we have
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 7
∂Φ(KNφ0)δ(Φ) = ∂N(KNφ0)(∂ΦN)δ(Φ) = ∂N(KNφ0)δ(Φ).(3.4)
Consequently, we can write
KNφ0∂Φδ(Φ) = N· ∇(KNφ0)δ(Φ).(3.5)
The three terms representing the linearised surface tension force account for the action
of the perturbation curvature on the basic state normal at the unperturbed interface and
vice versa, as well as the action of the perturbed interface on the basic state curvature and
normal. Similar terms can be identified in the linearised multi-domain formulation, there
multiplied with the perturbation of the interface height (see Tammisola et al. (2012)).
3.2. Implementation of linearised advection and diffusion terms
The computation of the linearised advection terms is facilitated using a slightly
modified version of the numerical scheme that is used for the nonlinear advection in
Basilisk. Here, we give a general overview of the scheme and necessary modifications.
For the calculation of the nonlinear advection term un+1/2·∇un+1/2in equation (2.10a),
a prediction of the velocity on the cell faces un+1/2
fis computed with a variant of the
BCG scheme (Bell et al. 1989), using the pressure gradient and body forces at t=n.
Solenoidality of the predicted velocity is enforced in a projection step. The advection step
itself is performed by an isolated function which takes, two arguments, a scalar field, in
this case the cell centred velocity field un, as well as un+1/2
f, to compute the advection
fluxes of u.
For the linearised advection step, a prediction of the basic state velocity is not needed,
since it is either stationary or known at all time steps (in case of a time-periodic base
flow). Thus, the two linear advection terms can be computed by calling the advection
routine twice, using u0nand Un+1
f, or Un+1 and u0n+1/2
f. A small technicality is involved
when advecting Un+1. The advection routine is implemented so that it by default adds
the action of the advection to the tracer fields that are advected, i.e. Un+1. However,
for both advection terms, the action needs to be added to u0n. To ensure this, minor
modifications of the routine are necessary.
The second advection term, stemming from the density variation between the fluids, is
computed again using the same BCG scheme. The diffusion term is added explicitly, as
it only contains the basic state velocity.
3.3. Implementation of linearised interface advection
An adaptation of Basilisk’s geometric VOF scheme for linear perturbation analysis is
not straightforward. A naive linearisation of the volume fraction transport equation and
the associated interface fluxes, as well as the geometric interface reconstructions would
retain the scheme’s inherent property of computing finite amplitude waves. The reason
is rooted in the way the interface location is defined in VOF methods. As seen from
equation (2.4), a cell where 0 < c < 1 contains an interface segment. As soon as this
cell becomes either full (c= 1) or empty (c= 0), the interface segment moves to one
of the neighbouring cells. Thus, without further modifications of the schemes, a growing
disturbance c0in a cell would eventually lead to c0= 1 and the disturbed interface would
move to an adjacent cell, producing a finite movement of the perturbed interface with
respect to the basic state interface and makes a evaluation of the interface displacement
difficult. As the level set function φis not bounded between 0 <φ<1, a disturbance φ0
8S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
can grow to any value. Its interpretation is therefore more straight forward: the value of
a disturbance φ0along the basic state interface is a measure for its displacement.
Another aspect that complicates the use of a linearised VOF method is the height-
function-based curvature computation in Basilisk. Depending on the orientation of the
interface, case distinctions are needed, regarding the choice of the component of the
height functions for computing the curvature (i.e. horizontal or vertical heights in a 2D
problem). The introduction of perturbed height functions, to calculate the curvature of
the perturbed interface, would add further case distinctions and complexity. By using a
level set representation of the interface, as introduced above, these problems are avoided.
Since the VOF method is used for the computation of the base flow volume fraction
C, an accurate reconstruction of the basic state level set field Φfrom Cis needed. In
general, one would need to compute the interface segments in each interfacial cell, where
0< C < 1, in order to compute the distance to the interface in each cell throughout
the domain. However, in the present study the basic state interface is close to horizontal,
and thus is consistently described by the height function field hy(x), used to compute the
basic state curvature (see §2.2). Therefore, we can use Φs=hy(x) as an initial, shifted
level set function, which is re-distanced by solving a Hamilton-Jacobi-type equation
(∂Φ
∂t+= sgn(Φs)(1 − |∇Φ|)
Φ(x, t+= 0) = Φs
.(3.6)
The shifted level set is integrated in a pseudo time t+until a steady state is reached.
The resulting level set Φis a signed distance function, such that |∇Φ|= 1. Note that,
for the nonlinear case, φloses its signed distance properties during advection. Thus, it
needs to be re-distanced every few time steps. In principle, this carries over to φ0in the
linear computations. Thus a linearised re-distancing function
(∂φ0
∂t+=−Φ∂φ0
∂x
∂Φ
∂x +∂φ0
∂y
∂Φ
∂y
φ0(x, t+= 0) = φ0
s
(3.7)
can be used to re-distance φ0
s. The resulting φ0has zero gradient in the base flow surface-
normal direction. Since the perturbed level set field corresponds to a displacement of
the basic state level set contours, after re-distancing, all level set contours are displaced
equally. However, we have found that re-distancing has no discernible influence on the
results of the present linear computations. In the nonlinear case, there are several reasons
for re-distancing φ. One aspect is that the delta function, defining the interface, is cal-
culated from φ. Consequently, as φloses its signed distance properties, the computation
of the interface location becomes inaccurate. Here, the linearised level set is only used
to compute the perturbed normal vector and curvature, since the base flow interface
is known a priori. Both computations were not affected by the re-distancing in our
tests. Another reason is the possible degeneration of numerical gradients if φdeviates
significantly from a signed distance function. Again, we have not found this to be a source
of error in our computations.
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 9
3.4. Implementation of the linearised surface tension force
For computing the linearised surface tension force, we utilise the volume fraction Cof
the basic state interface as well as the level set functions Φand φ0. In 2D, the normal
vectors are computed using the level sets as
N=∇Φ
|∇Φ|=∂Φ
∂x ,∂Φ
∂y
q(∂Φ
∂x )2+ (∂Φ
∂y )2
,(3.8a)
ζn0=∂ζφ0
∂x ,∂ζφ0
∂y
q(∂Φ
∂x )2+ (∂Φ
∂y )2
−∂Φ
∂x ,∂Φ
∂y ∂ζφ0
∂x
∂Φ
∂x +∂ζφ0
∂y
∂Φ
∂y
h(∂Φ
∂x )2+ (∂Φ
∂y )2i3/2.(3.8b)
Note that, although |∇Φ|= 1, we include the corresponding terms in the computation
of the normal vectors, as this leads to smoother numerical approximations. A detailed
derivation of n0is given in appendix A.
The basic state curvature Kis computed similarly to κin the nonlinear case using
height functions as described in §2.2. The curvature of the perturbed interface is com-
puted from Φand φ0as
κ0=∇ · (ζn0).(3.9)
The full expansion of κ0is given in appendix A. Further, we make use of the fact that
Nδ(Φ) = ∇Cand δ(Φ) = |∇C|,(3.10)
to compute terms which involve either Nδ(Φ) or δ(Φ). Consistently with the replace-
ment of H(Φ) by Cin §3.1, this choice as well retains an interface thickness of O(∆).
4. Formulation and solution of the eigenvalue problem
4.1. Eigenvalue Problem
The linear system (3.2) is rewritten in compact form as
∂q0
∂t =L(Q)q0,(4.1)
where q0= (u0, φ0)T, and L(Q) is the linearised Navier-Stokes operator, which might
be stationary or instationary, depending on the base flow Q= (U, Φ)T. For the former
case, we then seek eigenmodes of Lof the form
q0(x, y, z, t) = exp(λjt)ˆqj(x, y, z),(4.2)
where both λjand ˆqjare generally complex-valued. The corresponding eigenvalue
problem may be written as
Lˆqj=λjˆqj.(4.3)
This eigenvalue problem could in principle be solved directly, to give all eigenvalues, or
via various time-stepping techniques, which would yield one or several of the dominant
eigenvalues of the system, i.e. those of largest magnitude. However, the eigenvalues deter-
mining the stability of the system are the leading eigenvalues, which have the largest real
10 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
part, that in turn might be close to zero near a bifurcation point. Following Tuckerman
& Barkley (2000), the dominant eigenvalue can be extracted with the exponential power
method.
The general solution to (4.1) is
q0(t+τ) = Aq0(t) = exp(Lτ)q0(t),(4.4)
where the dominant eigenvalues of Aare the leading eigenvalues of L. Thus, we obtain
an alternative eigenvalue problem
Aˆqj=µjˆqj, µj= exp(λjτ).(4.5)
For a stationary base flow, the frequency of an eigenmode is given by the imaginary
part of the eigenvalue Im(λj) whereas its growth rate is given by the real part Re(λj).
Consequently, the system is linearly unstable when any Re(λj) is positive. Since for a
stationary Q, the stability is determined by the long-time behaviour of a perturbation,
τis an arbitrary value, which however should allow for a reasonable evolution of the
perturbation (Barkley et al. 2008).
4.2. Iterative Solution
Besides the power method which is able to recover the single most dominant eigenvalue,
the Arnoldi method (Saad 2011) is a suitable choice to recover a number of dominant
eigenvalues.
The method utilises a standard orthogonal projection of Aonto a lower dimensional
Krylov subspace that allows for an approximate direct evaluation of A. To this end, we
construct a sequence
Kn(A,q0) = [q0,Aq0, ..., An−1q0] (4.6)
which spans the Krylov subspace on which Ais projected. In practice, a direct
construction of Ais not needed in order to build Kn. We only need to be able to
compute the repeated action of Aon q0which is precisely described by equation (4.4).
The Krylov sequence can thus be updated by repeated time-stepping of the linearised
solver. The methodology to construct Knand to compute the corresponding eigenpairs is
that of the standard Arnoldi method with a modified Gram-Schmidt orthogonalisation.
Following Barkley et al. (2008), the part of the initial vector q0concerning the velocities
is prescribed as random fluctuations and solenoidality is enforced through the linearised
solver. We do not impose any initial fluctuations on the level set field, instead we let the
perturbation velocity generate those fluctuations consistently through equation (3.2b).
5. Global modes of a planar wake under the influence of surface
tension
To validate the derived framework, we revisit the work of Tammisola et al. (2012). As
noted above, this is one of very few works which study linear global modes of interfacial
flows without further simplifications of the linearised equations, contrary to what has
been done for instance in Rubio-Rubio et al. (2013). There, gravitational jets are studied
based on a 1-D long-wave model, derived by Eggers & Dupont (1994).
The flows studied in Tammisola et al. (2012) are planar, coflowing shear flows which
resemble jets and wakes, depending on the velocity ratio between inner and outer flow.
All configurations are confined flows, and for all but one flow condition, the flows are
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 11
(a) (b)
Figure 1: (a) Velocity profile at the inlet and interface separating the fluid phases at
subcritical conditions (We =∞). (b) Velocity field uand interface at the same conditions.
reported to be globally stable in the absence of surface tension. Only in the presence of
rather strong surface tension, a regime of global instability occurs.
5.1. Nonlinear simulation
In the study of Tammisola et al. (2012) no numerical or experimental validation of
the studied configurations is performed. Therefore, it remains to be shown whether the
computed linear modes really are present in the nonlinear flow, and if so, whether their
shapes and frequencies are accurately predicted by linear theory. We therefore present
nonlinear results for selected configurations of their study, with a focus on wake-type
flows.
All flows are computed in a long channel with a half-channel height h1+h2= 2 where
h1=h2= 1 are the heights of the respective fluid streams in the half-channel at the inlet,
and subscripts 1,2 refer to the inner and outer fluid stream, respectively. The resulting
confinement ratio is h=h2/h1= 1. The plug flow velocity ratio of both streams at the
inlet is Λ−1= (U1+U2)/(U1−U2) = −1.4. The density and viscosity ratios are ρ2/ρ1= 1
and µ2/µ1= 1. The configuration and flow field at subcritical conditions are shown in
figure 1. Upon nondimensionalisation of the linearised equations, using U2and h2, the
Reynolds number Re and the Weber number W e are given as
Re =ρ2U2h2
µ2
, We =ρ2U2
2h2
σ.(5.1)
We fix Re = 316 and vary We to track the behaviour of the computed instabilities.
For the simulations, we use a uniformly spaced grid of Nx×Ny= 1024 ×128 mesh
points, corresponding to 10 levels of refinement in stream wise direction and spanning a
nondimensional area of Lx×Ly= 32 ×4.
The domain is initialised with u= 0 and c= 0. At the left domain boundary Ω1an
inlet velocity is imposed such that
u|Ω1=(1,if |y|>1
1+Λ
1−Λ,if |y|⩽1,(5.2a)
∂p
∂x|Ω1= 0,(5.2b)
c|Ω1=(0,if |y|>1
1,if |y|⩽1.(5.2c)
12 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
For the right boundary Ω2, a standard outflow condition
∂u
∂x|Ω2= 0,(5.3a)
p|Ω2= 0,(5.3b)
∂c
∂x|Ω2= 0 (5.3c)
is used. The top and bottom boundaries Ω3are equipped with no-slip conditions. Time
stepping is adaptive, based on a Courant condition
∆t ⩽0.5 min (∆
u,rρ1∆3
πσ )(5.4)
where ∆denotes the grid spacing. The second term inside the braces accounts for the
propagation of capillary waves along the interface. Note that, especially for high surface
tensions and small cell sizes, this term is significantly more restrictive than the first term.
We monitor the energy of the system via the Euclidean norm ||u2||2to check whether
a steady state or a stable limit cycle has been reached. In the latter case a Dynamic
Mode Decomposition (DMD) is used for a modal decomposition of the flow field (Rowley
et al. 2009; Schmid 2010). The DMD yields an approximation of the Koopman operator,
a linear infinite dimensional operator describing a nonlinear dynamical system. As a
consequence, the frequencies associated with the DMD modes correspond to the modal
frequencies of the dynamical system. Further, with the mean flow subtracted, the DMD
modes are equivalent to modes of a Discrete Fourier Transformation (DFT), but they
usually require a lot less samples (Chen et al. 2012). A sequence composed of 2000
consecutive snapshots uand cevery ∆t = 0.5 for t > 1000 is used for computation of
the Ritz values λj. The dependence of the dominant DMD mode on the grid resolution
is given exemplarily in table 2 for We = 12.5
Following, Tammisola et al. (2012), we compute the nonlinear flow for a variety of
Weber numbers ranging from We = 3.¯
3 to We = 16.¯
6. In figure 2, the spectra obtained
from the DMD and the corresponding mode shapes with the largest amplitude that
dominate the flow are displayed for several investigated We where an unsteady solution
is found. The interface amplitude is overlayed on the modes. These are the same Weber
numbers that are presented in Tammisola et al. (2012).
For Weber numbers larger than W e ≈16, all disturbances decay, and the flow settles
to a steady solution. For We = 12.5, a periodic solution is seen that gives rise to a
sinuous mode with its amplitude maximum located around x= 8. For increasing surface
tension, at We = 10, three sinuous modes are observed. One mode is the same as seen for
We = 12.5, which has moved slightly upstream (amplitude maximum x= 7.5) and has a
lower frequency and amplitude. The second and third mode show a short-wave structure
in the vicinity of the inlet, with the amplitude maximum located at around x= 1.9,
and a long-wave structure located further downstream with an amplitude maximum
around x= 13. These modes have a higher (almost similar) frequency and a generally
larger amplitude than the first mode. Further, both of these modes display opposite
slight asymmetries in the amplitude in the short wave structure in the vicinity of the
inlet. Note, that in figure 2b only one of these modes is shown. For W e = 6.¯
6 three
modes are found. The most energetic mode has a varicose structure with an amplitude
maximum located at around x= 4. The other modes are the sinuous modes, observed
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 13
for We = 12.5 and W e = 10, with their amplitude maximum shifted far upstream. At
We = 5 a single varicose mode with its amplitude maximum close to the inlet at around
x= 2.6 is found. For even higher surface tension the flow becomes globally stable again,
so that for We = 3.¯
3 a steady solution is obtained.
To give a broader overview over the complexity of the dynamics in the flow and see
in which Weber number regime the dominant modes are unstable, their frequencies
are plotted versus the respective Weber numbers in figure 3a. Additionally, in figure
3b a bifurcation diagram is shown by plotting the local minima and maxima of the
transversal velocity vat (x, y) = (5,0), within the time period after the saturation.
From the above presentation it is seen that, under the influence of surface tension, the
flow undergoes a series of bifurcations that gives rise to an increasingly complex modal
interplay that governs the flow dynamics. Starting from a stable fixed point solution, a
first Hopf bifurcation occurs between We = 16.¯
6 and We = 12.5 where the flow reaches
a stable limit cycle, governed by a periodic oscillation produced by sinuous mode 1. In
figure 3b, this is characterised by a distinct minimum and maximum of v. A second
mode, sinuous mode 2 bifurcates between We = 12 and We = 11. Both sinuous modes
have incommensurate frequencies, resulting in a quasiperiodic motion and thus the limit
cycle looses its stability to a torus. As noted above, a separate mode very similar to
sinuous mode 2, with a slightly lower frequency is observed at We = 10. Both modes,
sinuous mode 1 and 2, are still active at We = 6.¯
6 and vanish between We = 6.¯
6
and We = 5. Between We = 9 and We = 8 a third mode with an incommensurate
frequency, varicose mode 1, emerges that governs the dynamics until the flow returns to
a stable fixed point again between We = 5 and We = 3.¯
3. Additionally, as is seen in
2a for We = 10 and We = 6.¯
6, a dense set of frequency peaks around the dominating
peaks, as well as some very low frequencies are seen in the spectra. While not shown,
this behaviour is also visible in the spectra for other intermediate Weber numbers.
Furthermore, for We = 12.5,5 a more pronounced influence of the first harmonic mode
is seen as compared to the other presented We. While a thorough investigation of these
additional modes and their dynamics is beyond the scope of this work, they are likely the
result of nonlinear interactions between the described modes and/or higher harmonics.
As a result of the presence of several competing modes for 11 ⩽W e ⩽6.¯
6, in this regime,
the flow is attracted towards higher-dimensional states, characterised by quasiperiodic
or chaotic oscillations. This is characterised in figure 3b by an increasingly dense variety
of minima and maxima at the respective Weber numbers. In §5.4 we will compare linear
and nonlinear dynamics.
5.2. Base flows
For the validation of the linear solver, we choose the basic state solution obtained for
Webase =∞, as it bears closest resemblance to the results of Tammisola et al. (2012)
who obtained the base flow by computing a single phase flow without interface and
surface tension, exploiting the fact that in the absence of surface tension the flow is
globally stable. The interface was introduced a posteriori by computing a streamline
from x= (0,1)T. It was argued that, for this particular configuration, the curvature
of the unperturbed interface is generally small, except for a small region very close to
the inlet (x < 0.1), and thus surface-tension induced pressure gradients are probably
negligible even for very high surface tension.
We have verified this assumption by computing the steady state solution at Webase =
∞, where the flow is globally stable, and Webase = 12.5, where the flow is only unstable
for sinuous perturbations and thus can be obtained by imposing symmetry conditions
at y= 0. Additionally, we have computed the interface as a streamline on the single
14 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
(a) (b)
Figure 2: (a) Magnitude of the DMD modes at each frequency, extracted from the
nonlinear simulation. (b) Shapes of the dominant modes at the respective Weber numbers.
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 15
(a) (b)
Figure 3: (a) Dominant frequencies of the appearing modes extracted from the DMD in
the flow for all investigated Weber numbers. (b) Bifurcations illustrated by the min-max
values of vat (x, y) = (5,0)
phase base flow. The interface coordinates of the first two base flows are computed
by reconstructing the piecewise linear interface segments from the respective volume
fractions using the same routines that are used by the VOF advection scheme. The
interface is plotted in figure 4a and the relative errors, with respect to the single phase
flow, where a streamline is computed as interface, are given in table 1. For both, Webase =
∞and Webase = 12.5, the error is below 1 %. The remaining error between the streamline
interface and the interface for Webase =∞is probably rooted in the inaccuracy of the
advection scheme as will be addressed in the next paragraph. The relative error of the
base flow interface for increasing grid resolutions is shown in table 2. As is seen, the
relative error between the highest and second highest resolution is again, below 1 %. For
the present flow configuration it is therefore justified to use the base flow at W ebase =∞
for the linear computations.
It has to be noted, that despite the marginal differences between the interfaces of the
basic states, the convergence of the case Webase = 12.5 is worse than for Webase =∞
or for the case of the single phase flow. While for the latter two cases an error ||un+1 −
un||∞=O(10−16) is reached, for the former it remains above O(10−6). The reason
for this behaviour is probably the occurrence of spurious currents. These are generally
not considered an issue for equilibrium solutions of interfacial flows without velocity of
the background fluid, when employing well-balanced schemes for computing the surface
tension term. However, in the dynamic case where the background fluid has non-zero
velocity, the problem persists due to numerical errors of the interface advection. These
induce curvature errors which, again, induce errors in the velocity. As demonstrated
in Popinet (2009) for a translating droplet, the induced error in the velocity field is
relatively insensitive to increased mesh resolution (its L∞norm shows less than first order
convergence) and scales approximately with We−1/2. The implications are readily seen
in the curvature plot in figure 4b. Since the κinvolves the second derivative of the height
function of the interface, the numerical errors of the advection scheme produce significant
oscillations in the curvature for Webase = 12.5. These are especially pronounced close
to the inlet where the shear of the fluid streams is largest. To a lesser degree, these
oscillations are also seen for Webase =∞, but since the surface tension force is zero,
there is no two-way coupling between the curvature and the velocity. Consequently, the
disturbances are weaker and the convergence of the velocity field is not affected.
16 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
||Is−I∞||/||Is|| ||Is−I12.5|/||Is||
0.0078 0.0086
Table 1: Relative error of the basic state interface I. Errors are measured against the
streamline-constructed interface of the single phase flow Is. The Weber numbers of the
respective two-phase solutions are stated as subscripts.
Level grid points (Nx×Ny)||I∞,Level −I∞,10||/||I∞,10|| Im(λ) (DNS) λ(linear analysis)
8 256 ×32 = 8192 0.0555 0.7392 0.0039 + 0.7431i
9 512 ×64 = 32768 0.0068 0.7414 0.0092 + 0.7416i
10 1024 ×128 = 131072 0 0.7414 0.0093 + 0.7416i
Table 2: Convergence of the nonlinear and linear flow for increasing level of refinement.
Shown is the relative error of the basic state interface for Webase =∞, the convergence
of the frequency Im(λ) of the dominant DMD mode of the nonlinear simulation and of
the most unstable eigenvalue λof the linear analysis for W e = 12.5.
(a) (b)
Figure 4: (a) Interface position and (b) curvature of the respective baseflows.
In preliminary computations we have assessed the effect of using finite Weber number
base flows, as the one for Webase = 12.5, for the linear computations. However, due to the
spurious advection errors in the base flow, spurious modes occurred in our results that
hampered a successful analysis. Abadie et al. (2015) have shown that advection of the
interface using a level set method may lead to a reduction of dynamic spurious currents
as opposed to the VOF method that is used here. This could be a possible avenue towards
reducing the base flow curvature errors for finite Weber numbers in a future study.
5.3. Computation of linear global modes
In section 5.1 we have seen that the action of surface tension on the base flow leads to
complex dynamics, with often several distinct modes being visible in the flow. Here, we
investigate whether linear analysis can successfully predict the destabilisation that gives
rise to the observed nonlinear dynamics. Further, we aim to compare the present linear
results with those of Tammisola et al. (2012). Therefore, we limit the analysis to the We
investigated in their study.
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 17
The computation of linear global modes is done by constructing a low-dimensional
subspace representation of the linearised Navier-Stokes operator, using the Arnoldi
method as outlined in §4.2. The Krylov vectors spanning the subspace are obtained by
time-stepping of the linearised solver. The integration time τbetween two consecutive
snapshots is arbitrary, however, in order to satisfy the Nyquist criterion, at least two
sample points are needed to resolve the highest desired frequency. Furthermore, too
small values of τmay lead to insufficient temporal separation between consecutive
Krylov vectors, thus hindering convergence of the method within a reasonable subspace
dimension. In practice we choose τ= 1.5. We have checked that the integration time,
within the stated restrictions, does not influence the obtained results to a relevant degree.
Arnoldi iterations are performed until a residual ||Aˆqj−µjˆqj|| <10−6for all unstable
eigenvalues is reached. The grid independence is demonstrated in table 2 for the unstable
eigenvalue, computed for W e = 12.5.
The computed growth rates Re(λ) and frequencies Im(λ) are given in figure 5 for the
stated Weber numbers. Both are compared to the results of Tammisola et al. (2012).
Further, the frequencies are compared to the mode frequencies obtained in the nonlinear
simulation. The area within the dashed vertical lines corresponds to the range of Weber
numbers where the nonlinear flow is unsteady. A direct comparison of the computed
eigenvalues of the present study and those computed in Tammisola et al. (2012) is given
in table 3. Further, in figure 6a the growth rates, obtained from the linear analysis for
We = 5,6.¯
6,10,12.5, are compared to the disturbance growth in the nonlinear simulation.
In figure 6b, the shapes of the computed unstable modes are presented for each Weber
number where the corresponding amplitude of the interface is plotted as black line and
where the amplitude is chosen arbitrarily for suitable visualisation.
The frequency values of the linear analysis and their nonlinear equivalents are very
similar in all cases. They compare especially well in the upper regime of unstable Weber
numbers, for We = 10 and We = 12.5. In both cases, all unstable mode are sinuous.
In line with the nonlinear simulation, sinuous mode 1 from We = 12.5 is still unstable
for We = 10 and W e = 6.¯
6 and the appearance of sinuous mode 2, that is the most
unstable mode for We = 10 and remains unstable for We = 6.¯
6, is accurately predicted.
The predicted mode shapes at both We show good agreement compared to the DMD
modes in figure 2b, both in terms of shape and stream wise position of the amplitude
maximum. For We = 6.¯
6, the most unstable mode is varicose mode 1. In contrast to the
computed DMD modes, the linear global modes at W e = 6.¯
6 show a more regular and
symmetric structure and their frequencies are slightly over-predicted. We believe this is
not a shortcoming of the linear analysis but rather an indication that the quasiperiodicity
of the dynamics induces non-negligible nonlinear interactions in the flow that lead to a
departure of from the linear dynamics presented here. For We = 5, varicose mode 1 is
the only unstable mode predicted. Its frequency is, again, slightly over-predicted which
again can be attributed to a non-negligible nonlinearity in the flow dynamics.
Turning to the growth rates, in figure 6a a comparison of the present growth rates
with the energy growth of the nonlinear simulation is shown. Therefore, we initialise the
nonlinear flow from the steady base flow with Webase =∞and monitor the disturbance
growth via the v component of the velocity at (x, y) = (5,0) where V= 0. The plotted
time series of the nonlinear energy growth is obtained by computing the logarithm of
the envelope of the modulus of the monitored velocity signal. The lines representing the
growth rates of the linear analysis are shifted for suitable alignment with the nonlinear
time series. For We = 12.5, a sufficiently long time interval of approximately exponential
growth is seen that is in agreement with the obtained growth rates. For We = 10
18 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
(a) (b)
Figure 5: (a) Comparison of the mode frequencies and (b) growth rates at the respective
Weber numbers.
We We−1λ(present study) λ(Tammisola et al.) Im(λ) (DNS)
16.¯
6 0.06 -0.0042 + 0.7851i -0.01 + 0.788i stable
12.5 0.08 0.0093 + 0.7416i 0.005 + 0.743i 0.7414
10 0.1 0.0214 + 0.8422i 0.021 + 0.843i 0.8419
0.0194 + 0.6946i 0.015 + 0.694i 0.6912
0.0162 + 0.831i - 0.8294
6.¯
6 0.15 0.0250 + 0.6145i 0.02 + 0.63i 0.5781
0.0243 + 0.5819i - 0.5529
0.0097 + 0.5653i 0.007 + 0.567i 0.5278
5 0.2 0.0161 + 0.4413i 0.017 + 0.453i 0.4021
3.¯
3 0.3 -0.007 + 0.2896i 0.008 + 0.302i stable
Table 3: Eigenvalues λof the present linear analysis, the ones of Tammisola et al. (2012)
and frequencies Im(λ) of the dominant DMD modes of the nonlinear simulation at the
respective Weber numbers.
and We = 6.¯
6, intervals of exponential growth are less clear which may be due to
the competing modes that are observed at this Weber numbers and that prevent a
clear exponential growth of a single mode from being observed. Furthermore it could be
hypothesised that transient growth leads to short time amplification of additional modes
that further influence the observed energy growth. Nevertheless, the modal growth rates
obtained in the linear analysis agree reasonably well with the overall energy growth,
observed in the nonlinear flow at We = 10,6.¯
6. For We = 5, an interval of clear
exponential growth is seen that is in good agreement with the predicted growth rate
of the linear analysis.
Comparing the present results to the study of Tammisola et al. (2012), good agreement
is seen between the computed mode frequencies and those of Tammisola et al. (2012).
In particular, both linear analyses show similar over-predictions of the frequencies for
We = 5,6.¯
6. On the other hand, the present growth rates show notable deviations,
compared to their results. For We = 6.¯
6,10,12.5,16.¯
6 the present growth rates are always
higher, up to almost a factor of 2, whereas for We = 3.¯
3 they are notably lower. The
prediction of the bifurcation points also differs. For We = 3.¯
3, Tammisola et al. (2012)
still predict an unstable eigenvalue, while in the present study, the flow is already linearly
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 19
(a) (b)
Figure 6: (a) Comparison of the velocity disturbance growth of the nonlinear simulation v
at (x, y) = (5,0) and the exponential growth rate of the most unstable linear global modes
of the present analysis. (b) Mode shapes of the respective eigenmodes. The corresponding
eigenvalue λis given for each mode.
20 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
stable, consistent with the nonlinear simulation. Several reasons for the different growth
rates in both studies may be found. Generally, the growth rates obtained in both analyses
are rather small for all computed Weber numbers. Therefore, slight differences in the
employed numerical schemes may lead to large relative changes in the growth rates.
From a numerical perspective, the treatment of the interface and surface tension via
separate domains and coupling conditions in Tammisola et al. (2012) might be a source
for differences. Further, the use of Chebyshev polynomials in their study, while spectrally
accurate, introduces a non-uniform grid spacing which leads to a reduced resolution in
the centre of the respective domains. Another possible reason is the use of a streamline-
based reconstruction of the interface which is shown to result in a slightly smaller interface
curvature, as seen in figure 4b. These differences are especially prominent close to the
inlet, where the global mode at this Weber number is localised.
5.4. Comparison of nonlinear and linear dynamics
In §5.1 and §5.3 we have analysed the nonlinear and linear dynamics of plane wakes
under the influence of surface tension and have found that the linear analysis reliably
predicts the dominant modes, seen in the nonlinear flow, to be linearly unstable. However,
from figure 2a it is as well seen that for the intermediate Weber numbers studied there, the
nonlinear dynamics contain significantly more modes than the linear analysis predicts.
Hence, these modes are the result of nonlinear harmonic interactions resulting from the
quadratic nonlinearity of the advection term or the higher order geometric nonlinearity
of the surface tension term. It is therefore worthwhile to shed more light on the validity of
the linear dynamics to approximate the nonlinear flow. To this end, we use the dominant
DMD modes that are shown in figure 2b, and which we have seen to closely correspond
to the linear modes presented figure 6b and table 3, to construct a low dimensional
approximation of the nonlinear dynamics of the flow.
The departure of the linear dynamics towards the full nonlinear dynamics can be
analysed by employing the time-averaged mean flow as has been shown in various works
in recent times (e.g. Barkley 2006; Oberleithner et al. 2014; Boujo et al. 2018; Schmidt &
Oberleithner 2020). According to the theory established therein, the nonlinear saturation
process is driven by the interaction of higher harmonics with their complex conjugates
through the Reynolds stress divergence that modifies the base flow towards the mean
flow. As a consequence, eigenmodes that are unstable with respect to the base flow
become marginally stable on the mean flow. Simultaneously, the mean flow modification
may alter the structure and frequency of the eigenmodes through the linear perturbation
equation. The remaining difference to the full nonlinear flow is accounted for by nonlinear
interactions between different harmonics. In the present work, the frequencies of the linear
modes are in close correspondence to their equivalent nonlinear modes which indicates
that these modes are not significantly affected by the mean flow modification. Hence, the
constructed model can be seen as a representation of the linear dynamics around the mean
flow. The difference between the reduced and full dynamics are the nonlinear interactions
between different harmonics and thus, it is reasonable to assume that a comparison of
the reconstructed and full dynamics yields an appropriate qualitative picture of how far
linear and nonlinear dynamics are apart.
In figure 7, we present velocity phase trajectories of the full and reconstructed dynamics
using probes located at the approximate stream wise location of the respective amplitude
maxima of the oscillation of the nonlinear flow. For a clean visual representation, we
plot the velocity histories using a time lag of one period τ= 1/Im(λ) where Im(λ) is
chosen as the frequency of the mode with the largest amplitude for each respective Weber
number (i.e. sinuous mode 1 (We = 12.5), sinuous mode 2 (W e = 10), varicose mode
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 21
Figure 7: Qualitative representation of the velocity phase trajectories at (x, y) = (8,0)
(We = 12.5), (x, y) = (8,0) (We = 10), (x, y) = (6,0) (We = 6.¯
6) and (x, y) = (4,0)
(We = 5). The upper row in black shows the full dynamics and the lower row in blue
the reconstructed dynamics. All plots are equally scaled.
1 (We = 6.¯
6,5)). To get a better view of the cross-sections of the phase space we plot
a variant of Poincar´e sections in figure 8. The section is chosen such that the phase
trajectory crosses it every completed cycle n, established by the period τ.
As is expected from the analysis in §5.1, for We = 12.5 and W e = 5, the full
nonlinear dynamics are represented by a limit cycle. Thus a closed orbit is seen in figure
7 which, however, shows a notable deformation, caused by higher harmonic modes. The
reconstructed mono-frequent dynamics lack this higher harmonic influence which results
in a slightly smaller orbit and thus a reduced velocity amplitude. As expected for a limit-
cycle solution, both trajectories cross the section plane at a fixed point in figure 8. For
We = 10 and We = 6.¯
6, the trajectories of the full nonlinear dynamics show complex
patterns and the orbits fill a closed region in the phase space, possibly indicating chaotic
behaviour. The reconstructed dynamics each consist of three incommensurate frequencies,
thus forming a quasiperiodic cycle on a 3-torus. The reconstruction shows a significantly
different behaviour compared to the full dynamics and particularly for We = 10 results
in a overall lower amplitude oscillation. For both Weber numbers, the section plots in
figure 8 show the trajectory crossings clustered broadly around the fixed points, observed
for We = 12.5,5. While the crossings of the reconstructed dynamics show regular closed
paths, indicative of quasiperiodicity, the full dynamics do not exhibit clear paths, again
suggesting a chaotic behaviour.
5.5. Discussion and relations to other studies
The linear stability results obtained by Tammisola et al. (2012) are counter-intuitive,
in so far as they find surface tension to be destabilising in a plane two-dimensional flow,
contrary to its common stabilising tendency in plane flows (Squire 1953; Hagerty & Shea
1955). Furthermore, at the studied Reynolds number, the velocity ratio of the adjacent
fluid streams is too low to trigger the B´enard-von K´arm´an instability which is typically
present in wake flows at higher Reynolds number or velocity ratio. As a result, the flow
is globally stable in the absence of surface tension. Viscosity and density being equal
across the fluids in the present configuration, the only possible destabilizing force is
surface tension. In fact, a prospect of a destabilizing surface tension force in plane flows
is given in studies by Rees & Juniper (2009) and Biancofiore & Gallaire (2010), based on
inviscid local stability analyses, suggesting that such effect, induced by a top-hat velocity
profile, exists in the high-Re limit. It was found that for certain configurations, moderate
22 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
Figure 8: Qualitative Poincar´e sections of the trajectories presented in figure 7. Sections
are chosen such that trajectories cross every completed cycle n. The upper row in black
shows the full dynamics and the lower row in blue the reconstructed dynamics. All plots
are equally scaled.
surface tensions both, wakes and jets, are significantly more unstable than without surface
tension. To further explore this topic, we have conducted detailed nonlinear and linear
computations of this configuration.
The linear analysis is able to predict all dominant modes present in the nonlinear flow,
even for We where several competing modes govern the dynamics. As a consequence, for
the investigated Weber numbers close to the bifurcation points, i.e. W e = 12.5,5, the
linear and nonlinear dynamics are in close correspondence. However, for intermediate
Weber numbers, e.g. W e = 10,6.6, the full nonlinear dynamics have been shown to
be significantly more complex, involving significant nonlinear harmonic interactions.
Therefore, in these regimes, discrepancies between linear and nonlinear dynamics are
found.
The presently obtained linear modes are also in very good agreement with the results of
the linear stability analysis of Tammisola et al. (2012). On the other hand, the computed
growth rates show notable differences between both analyses. However, the comparison
of nonlinear and linear analysis suggests that the presently obtained growth rates are
plausible. In the study of Tammisola et al. (2012), the bifurcation occurring between
We = 16.¯
6 and We = 12.5 is well captured, as it is in the present results. The stabilisation
of the global mode at We = 3.¯
3 is only captured in the present analysis, whereas in their
work the linear analysis still shows a slightly unstable mode. Possible reasons for this
discrepancy have been given in §5.3.
Comparing the present results to the stability maps for the uniform density case of
Rees & Juniper (2009) (Figs. 7 and 8 in their work), it is seen that our configuration is in
a region of absolute instability for We−1>0.1 (varicose modes) and 0.1< We−1<10
(sinuous modes). Note, that maps for We−1<0.1 are not considered in Rees & Juniper
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 23
(2009) but the inviscid stability limits should approximately remain similar as We−1
tends to zero, as is also seen in Fig. 9 of Biancofiore & Gallaire (2010). While their
inviscid results cannot be expected to translate accurately to the present viscous flow
results, the general trend of varicose modes remaining unstable for higher surface tensions,
while sinuous modes becoming stable again, persists. The dampening effect of viscosity in
the present flow eventually leads to stabilisation of also varicose modes below We ≈3.¯
3.
Similarly, viscosity renders the flow stable above We ≈16.¯
6, while in the inviscid limit
it remains unstable.
The physical mechanism of surface tension-induced instability of plane shear flows has
been explored in Biancofiore et al. (2015, 2017) from a Kernel wave perspective. The
resulting system consists of two counter-propagating inertial and capillary waves. As
a general concept, the vorticity perturbation of each wave induces an effect on its
counter-propagating wave. Depending on their phase angle (and speed) a dampening
or amplification is observed. Without surface tension the studied system showed a single,
purely inertial, unstable mode at large wavelengths. With the inclusion of surface tension
a second mode appeared, acting at smaller wavelengths. It was further found, that
at large wavelengths influence was mainly due to an interaction of the inertial waves,
while for small wave numbers interaction was between one inertial wave and its counter-
propagating capillary wave. This was explained using a phase-locking of the two counter-
propagating waves, resulting in a continuous amplification. Contrastingly, at smaller
wavelengths, the phase speed was found too high for prolonged interaction to take place
effectively, thus limiting the capillary-inertial interaction to small wavelengths. For an
increasing surface tension, the phase speed reduces and the upper cutoff wavelength is
shifted towards larger wavelengths. Although their studied configurations do not match
the present one, these findings translate conclusively to the present results: The long-
wave sinuous mode 1 corresponds to an inertial wave interaction, likely a von K´arm´an
mode, modified by surface tension to become unstable. Similarly, the short-wave sinuous
mode 2 as well as the varicose mode 1 correspond to an inertial-capillary wave interaction
with increasing wavelength for increasing surface tension.
In light of the present results, it is important to discuss the results of another related
study by Biancofiore et al. (2014), which investigated the same flow configuration with
direct numerical simulation, using a finite element method in conjunction with a level set
method for interface capturing. Their study confirms the global stability of the flow in
the absence of surface tension and shows comparable results for some of the cases studied
in the present work. In particular, for We = 10, they detect a clear sinuous disturbance
with a velocity amplitude of 10−2which is localised between 5 <x<15. The periodicity
of this oscillation is T= 9 which corresponds to a mode frequency ω= 0.698. The
spatial appearance and frequency of this instability are close to the sinuous mode 1
computed at this Weber number in the present study and by Tammisola et al. (2012). In
contrast, for We = 12.5, they detect a weak sinous disturbance in form of an intermittent,
transient wave packet that has a disturbance amplitude of the velocity of about 10−4
and a disturbance amplitude of the interface of about 10−8. Although it is not explicitly
mentioned, the transient character of this disturbance leads to the conclusion that this
is a convective instability. The flow at this Weber number still appears to be globally
stable. For We = 5, a very weak disturbance in close vicinity of the inlet is detected,
which, however, seems to have insufficient amplitude for a reliable frequency estimation.
Moreover, it seems only detectable when axial symmetry is imposed along the channel
axis.
Biancofiore et al. (2014) argue that the differences might stem from the fact that
24 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
a single phase flow without surface tension was computed for the basic states used
in Tammisola et al. (2012), neglecting pressure modifications due to surface tension.
In particular, they point out that the incompatibility of the plug flow inlet condition
for the velocity with the wall boundaries leads to a very localised contraction of the
interface immediately behind the inlet. The resulting curvature thus induces a non-
negligible pressure gradient. While we agree with this rationale, the present nonlinear
results suggest that this pressure gradient is not the cause for the observed differences
since the modes predicted in Tammisola et al. (2012) are clearly present in our nonlinear
simulations which account for the pressure gradient. It is not entirely clear what causes
the differences in the results of Biancofiore et al. (2014). However, one might speculate
that, for instance, numerical diffusion caused by the employed numerical schemes or
insufficient mesh resolution could attenuate some of the observed modes, leaving only
the sinuous mode 1 unstable. Additionally, this could cause the bifurcation points of the
detected unstable mode to move close together such that the flow remains only unstable
around We = 10.
6. Conclusions
In this paper we have accomplished three objectives.
First, we have presented a framework for the computation of linear global modes of
interfacial flow, by means of time-stepping of a linearised Navier-Stokes solver. The effects
of surface tension, density and viscosity differences between the different phases are fully
accounted for, thereby avoiding a few simplifications that were made in previous studies of
similar interfacial flow situations. The matrix-free time-stepper-based method also avoids
excessive memory requirements, so that the approach is fairly straightforward to apply
to more complex geometries or to three-dimensional flow. As a perspective for future
improvements, robust computations of finite Weber number or non-uniform viscosity or
density base flows are needed. This includes the stabilisation of the unstable flow by
methods like Newton iterations or selective frequency damping as well as a reduction of
spurious advection errors by improved numerical schemes.
Second, we have revisited the work of Tammisola et al. (2012) and conducted nonlinear
simulations to rigorously validate the results of their global linear stability analysis.
Further, we have reproduced the computations with the developed linear solver. The
present linear analysis is in very good agreement with the results of Tammisola et al.
(2012) and both linear predictions are confirmed by the present nonlinear simulations.
The present method accurately captures the bifurcation points as deduced from the
nonlinear analysis and growth rates show good agreement with the initial disturbance
growth found in the nonlinear simulation.
Third, our validation of the configuration provides clarity that a surface-tension-induced
destabilisation of plane wakes can indeed be observed in the nonlinear flow. It is shown
that, in the studied flow regimes, surface tension can induce both, varicose and sinuous
oscillations of the flow which may be present simultaneously at several We, thus leading
to complex quasiperiodic or chaotic dynamics.
Funding. We gratefully acknowledge the support of the Elsa-Neumann-Stipendium of
the state Berlin, Germany.
Declaration of interests. The authors report no conflict of interest.
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 25
Appendix A. Derivation of the linearised normal vector and
curvature
Here we give a detailed derivation of the linearised normal vector and curvature, based
on the level set function φ. As introduced in §3, the decomposition into basic state and
infinitesimal perturbation is given as φ=Φ+ζφ0, and the nonlinear normal vector is
given as
n=∇φ
|∇φ|=∂Φ
∂x +ζ∂φ0
∂x ,∂Φ
∂y +ζ∂φ0
∂y
∂Φ
∂x +ζ∂φ0
∂x 2+∂Φ
∂y +ζ∂φ0
∂y 21/2.(A 1)
Then the basic state normal is
N=∇Φ
|∇Φ|=∂Φ
∂x ,∂Φ
∂y
∂Φ
∂x 2+∂Φ
∂y 21/2,(A 2)
and the perturbed normal vector is found by linear Taylor expansion of naround = 0
as
ζn0=ζ∂n
∂ζ ζ=0
(A 3)
=ζ∂φ0
∂x ,∂φ0
∂y "∂Φ
∂x 2
+∂Φ
∂y 2#−1/2
+
ζ∂Φ
∂x ,∂Φ
∂y ∂
∂ζ "∂Φ
∂x +ζ∂φ0
∂x 2
+∂Φ
∂y +ζ∂φ0
∂y 2#−1/2ζ=0
=∂ζφ0
∂x ,∂ζφ0
∂y "∂Φ
∂x 2
+∂Φ
∂y 2#−1/2
+
∂Φ
∂x ,∂Φ
∂y ∂ζφ0
∂x
∂Φ
∂x +∂ζφ0
∂y
∂Φ
∂y "∂Φ
∂x 2
+∂Φ
∂y 2#−3/2
.
If we introduce
f(Φ) = "∂Φ
∂x 2
+∂Φ
∂y 2#−1/2
,(A 4a)
g(Φ) = −∂ζφ0
∂x
∂Φ
∂x +∂ζφ0
∂y
∂Φ
∂y ,(A 4b)
we can write the perturbed normal vector as
ζn0=∂ζφ0
∂x f, ∂ζφ0
∂y f+gf3∂Φ
∂x , gf3∂Φ
∂y .(A 5)
26 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
The curvature is the divergence of the normal vector. Consequently, the curvature of
the basic state is given as
K=∇ · N=
∂2Φ
∂x2∂Φ
∂y 2−2∂2Φ
∂x∂y
∂Φ
∂x
∂Φ
∂y +∂2Φ
∂y2∂Φ
∂x 2
h∂Φ
∂x +∂Φ
∂y i3/2.(A 6)
Similarly, the curvature of the perturbed state is found as
κ0=∇ · ζn0(A 7)
=∂f
∂x
∂ζφ0
∂x +f∂2ζφ0
∂x2+∂f
∂y
∂ζφ0
∂y +f∂2ζφ0
∂y2
+∂g
∂xf3∂Φ
∂x +g3f2∂f
∂x
∂Φ
∂x +gf3∂2Φ
∂x2
+∂g
∂y f3∂Φ
∂y +g3f2∂f
∂y
∂Φ
∂y +gf3∂2Φ
∂y2
where
∂f
∂x =−f3∂Φ
∂x
∂2Φ
∂x2+∂Φ
∂y
∂2Φ
∂x∂y ,(A 8)
∂f
∂y =−f3∂Φ
∂x
∂2Φ
∂x∂y +∂Φ
∂y
∂2Φ
∂y2,(A 9)
∂g
∂x =−∂ζφ0
∂x
∂2Φ
∂x2+∂2ζφ0
∂x2
∂Φ
∂x +∂ζφ0
∂y
∂2Φ
∂x∂y +∂2ζφ0
∂x∂y
∂Φ
∂y ,(A 10)
∂g
∂y =−∂ζφ0
∂y
∂2Φ
∂y2+∂2ζφ0
∂y2
∂Φ
∂y +∂ζφ0
∂x
∂2Φ
∂y∂x +∂2ζφ0
∂y∂x
∂Φ
∂x .(A 11)
REFERENCES
Abadie, T., Aubin, J. & Legendre, D. 2015 On the combined effects of surface tension force
calculation and interface advection on spurious currents within volume of fluid and level
set frameworks. J. Comput. Phys. 297, 611–636.
Barkley, D. 2006 Linear analysis of the cylinder wake mean flow. EPL (Europhys. Lett.) 75 (5),
750.
Barkley, D., Blackburn, H. M. & Sherwin, S. J. 2008 Direct optimal growth analysis for
timesteppers. Intl. J. Numer. Meth. Fluids 57 (9), 1435–1458.
Bell, J. B., Colella, P. & Glaz, H. M. 1989 A second-order projection method for the
incompressible Navier–Stokes equations. J. Comput. Phys. 85 (2), 257–283.
Biancofiore, L. & Gallaire, F. 2010 Influence of confinement on temporal stability of plane
jets and wakes. Phys. Fluids 22 (1), 014106.
Biancofiore, L., Gallaire, F. & Heifetz, E. 2015 Interaction between counterpropagating
rossby waves and capillary waves in planar shear flows. Phys. Fluids 27 (4), 044104.
Biancofiore, L., Gallaire, F., Laure, P. & Hachem, E. 2014 Direct numerical simulations
of two-phase immiscible wakes. Fluid Dyn. Res. 46 (4), 041409.
Biancofiore, L., Heifetz, E., Hoepffner, J. & Gallaire, F. 2017 Understanding the
destabilizing role for surface tension in planar shear flows in terms of wave interaction.
Phys. Rev. Fluids 2(10), 103901.
Boeck, T. & Zaleski, S. 2005 Viscous versus inviscid instability of two-phase mixing layers
with continuous velocity profile. Phys. Fluids 17 (3), 032106.
Global stability and nonlinear dynamics of wake flows with a two-fluid interface 27
Boujo, E., Bauerheim, M. & Noiray, N. 2018 Saturation of a turbulent mixing layer over a
cavity: response to harmonic forcing around mean flows. J. Fluid Mech. 853, 386–418.
Brackbill, J. U., Kothe, D. B. & Zemach, C. 1992 A continuum method for modeling
surface tension. J. Comput. Phys. 100 (2), 335–354.
Chen, K. K., Tu, J. H. & Rowley, C. W. 2012 Variants of dynamic mode decomposition:
boundary condition, Koopman, and Fourier analyses. J. Nonlinear Sci. 22 (6), 887–915.
Chorin, A. J. 1969 On the convergence of discrete approximations to the Navier–Stokes
equations. Math. Comput. 23 (106), 341–353.
Cummins, S. J., Francois, M. M. & Kothe, D. B. 2005 Estimating curvature from volume
fractions. Comput. Struct. 83 (6-7), 425–434.
Drazin, P. G. & Reid, W. H. 2004 Hydrodynamic Stability. Cambridge University Press.
Eggers, J. & Dupont, T. F. 1994 Drop formation in a one-dimensional approximation of the
Navier–Stokes equation. J. Fluid Mech. 262, 205–221.
Hagerty, W. & Shea, J. F. 1955 A study of the stability of moving liquid film. ASME J.
Appl. Mech 22, 509–514.
Harvie, D. J. E., Davidson, M. R. & Rudman, M. 2006 An analysis of parasitic current
generation in volume of fluid simulations. Appl. Math. Model. 30 (10), 1056–1066.
Oberleithner, K., Rukes, L. & Soria, J. 2014 Mean flow stability analysis of oscillating jet
experiments. J. Fluid Mech. 757, 1–32.
Peskin, C. S. 1972 Flow patterns around heart valves: a numerical method. J. Comput. Phys.
10 (2), 252–271.
Popinet, S. 2003 Gerris: A tree-based adaptive solver for the incompressible Euler equations
in complex geometries. J. Comput. Phys. 190 (2), 572–600.
Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. J.
Comput. Phys. 228 (16), 5838–5866.
Popinet, S. 2018 Numerical models of surface tension. Annu. Rev. Fluid Mech. 50, 49–75.
Rayleigh, L. 1878 On the instability of jets. Proc. London Math. Soc. s1-10 (1), 4–13.
Rees, S. J. & Juniper, M. P. 2009 The effect of surface tension on the stability of unconfined
and confined planar jets and wakes. J. Fluid Mech. 633, 71.
Rowley, C. W., Mezi´
c, I., Bagheri, S., Schlatter, P. & Henningson, D. S. 2009 Spectral
analysis of nonlinear flows. J. Fluid Mech. 641, 115–127.
Rubio-Rubio, M., Sevilla, A. & Gordillo, J. M. 2013 On the thinnest steady threads
obtained by gravitational stretching of capillary jets. J. Fluid Mech. 729, 471–483.
Saad, Y. 2011 Numerical Methods for Large Eigenvalue Problems: Revised Edition. SIAM.
Scardovelli, R. & Zaleski, S. 1999 Direct numerical simulation of free-surface and interfacial
flow. Annu. Rev. Fluid Mech. 31 (1), 567–603.
Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. J. Fluid
Mech. 656, 5–28.
Schmid, P. J. & Henningson, D. S. 2012 Stability and transition in shear flows. Springer
Science & Business Media.
Schmidt, S. & Oberleithner, K. 2020 Instability of forced planar liquid jets: mean field
analysis and nonlinear simulation. J. Fluid Mech. 883.
Sevilla, A. 2011 The effect of viscous relaxation on the spatiotemporal stability of capillary
jets. J. Fluid Mech. 684, 204–226.
S¨
oderberg, L. D. 2003 Absolute and convective instability of a relaxational plane liquid jet.
J. Fluid Mech. 493, 89.
Squire, H. B. 1953 Investigation of the instability of a moving liquid film. Br. J. Appl. Phys.
4(6), 167.
Sussman, M., Smereka, P., Osher, S. & others 1994 A level set approach for computing
solutions to incompressible two-phase flow. J. Comput. Phys. 114 (1), 146 – 159.
Tammisola, O., Loiseau, J.-C. & Brandt, L. 2017 Effect of viscosity ratio on the self-
sustained instabilities in planar immiscible jets. Phys. Rev. Fluids 2, 033903.
Tammisola, O., Lundell, F. & S¨
oderberg, L. D. 2011 Effect of surface tension on global
modes of confined wake flows. Phys. Fluids 23 (1), 014108.
Tammisola, O., Lundell, F. & S¨
oderberg, L. D. 2012 Surface tension-induced global
instability of planar jets and wakes. J. Fluid Mech. 713, 632–658.
28 S. Schmidt, O. Tammisola, L. Lesshafft and K. Oberleithner
Tuckerman, L. S. & Barkley, D. 2000 Bifurcation analysis for timesteppers. In Numerical
methods for bifurcation problems and large-scale dynamical systems, pp. 453–466. Springer.
Yih, C. 1967 Instability due to viscosity stratification. J. Fluid Mech. 27 (2), 337–352.
