University of Southern California Dissertations and Theses

Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity

(USC Thesis Other)

Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Hydrodynamic stability of two fluid columns of different
densities and viscosities subject to gravity
by
Aditya Heru Prathama
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Mechanical Engineering)
December 2021
Copyright 2021 Aditya Heru Prathama
To my wife Valencia,
for her love and encouragement
ii
Acknowledgements
First and foremost, I would like to express my gratitude to my advisor, Prof. Carlos Pantano-
Rubino, for his guidance, encouragement and patience over the past few years. Without his invalu-
able insight and enthusiastic involvement, this thesis would not have been possible. I would also
like to thank my dissertation and qualifying exam committee members, Profs. Mihailo Jovanovic,
Mitul Luhar, Ivan Bermejo-Moreno and Alejandra Uranga, for their constructive comments and
discussions to improve the dissertation and enhance my understanding of the problem.
This work is motivated by the DNS study carried out at Caltech by Dr. Ilana Gat, so I would
like to thank her for the fruitful collaboration particularly during the development of the LSA part.
In addition, this thesis would not be completed without the tremendous assistance of Dr. Xiaoyi
Lu, who guided me through the simulation part and helped me to understand the zMachDNS
code better. I cherish the friendship and opportunity to work together with the current and former
members of the High-Performance Turbulence and Combustion Modeling Lab – David Petty, Hang
Yu, Qingquan Wang, Timothy Smith, Ben Shields, Xiaoyi Lu and Sicong Wu.
A PhD journey would not be as meaningful and joyful without the support from all friends and
families in Champaign, Los Angeles, Surabaya and Manado. First of all I would like to thank my
parents for the support and encouragement from afar. Most importantly, I would like to appreciate
the love and encouragement from my wife Valencia Ang throughout my study, especially in the
moments that she had to take care of our kids solitarily. Finally, I would like to cherish my children,
Hans and Erin, for the joy and delightful interruptions they bring each day.
iii
This work is supported in part by the the Department of Energy, National Nuclear Security
Association, under award no. DE-NA0002382, and the California Institute of Technology. Com-
puting resources are provided by the USC’s Center for Advanced Research Computing. In partic-
ular, I would like to thank Cesar Sul and Marco Olguin for their help in running the code in the
clusters. Last but not least, I would like to thank Chrissy Franks and Natalie Guevara from AME
Department for their help in facilitating all academic matters during my time at USC.
iv
Contents
Dedication ii
Acknowledgements iii
Abstract xi
Chapter 1: Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Motivation and scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Chapter 2: Inviscid linear stability analysis 8
2.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Flow analysis for small perturbations . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
Chapter 3: Viscous linear stability analysis 21
3.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Flow analysis for small perturbations . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.1 Base flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.2 Initial-value problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
v
3.3.1 Two-dimensional perturbations . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.2 Three-dimensional perturbations . . . . . . . . . . . . . . . . . . . . . . . 39
3.4 Comparison with DNS initialized with random perturbation fields . . . . . . . . . 45
Chapter 4: Nonlinear simulation 47
4.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3 Simulation setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Chapter 5: Conclusions 70
References 73
Appendices 77
A Detailed perturbation equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
A.1 Linear IVP – forward equations . . . . . . . . . . . . . . . . . . . . . . . 78
A.2 Quasi-steady-state approximation . . . . . . . . . . . . . . . . . . . . . . 81
A.3 Linear IVP – adjoint equations . . . . . . . . . . . . . . . . . . . . . . . . 82
B Inviscid solution for nonzero velocity difference Δw
0
at t = t
0
> 0 . . . . . . . . . 85
vi
List of Figures
2.1 Planar sketch of flow configuration in Cartesian coordinates for the inviscid analy-
sis (y is perpendicular to the page). . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Evolution of the interface amplitude
ˆ
ζ for various values of s. . . . . . . . . . . . . 17
2.3 Magnitudes (a,c) and phases (b,d) of
ˆ
C
∗−
and
ˆ
C
∗+
for density ratio R = 4/3 for
various values of s; legend as in figure 2.2. . . . . . . . . . . . . . . . . . . . . . . 19
2.4 Real and imaginary values of
ˆ
ϕ
+
for density ratio R= 4/3 and various values of s;
legend as in figure 2.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1 Planar sketch of flow configuration in Cartesian coordinates (y is perpendicular to
the page). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Base flow fields as a function of density ratio R at uniform viscosity M = 1. . . . . 27
3.3 Base flow fields as a function of viscosity ratio M at R= 5. . . . . . . . . . . . . . 28
3.4 Amplification factors G
T
, G
K
and G
Y
with different initializations of E
Y
(τ
0
) = 1
(left), E
K
(τ
0
)= 1 (center) and E
T
(τ
0
)= 1 (right), respectively, for R= 5, M = 1,
τ
0
= 10, m= 0.05 and k= 0 (two-dimensional). . . . . . . . . . . . . . . . . . . . 34
3.5 The E
Y
-optimal (solid) versus forward (dashed) amplification factors for R = 5,
M = 1, m= 0.05 and τ
0
= 1, 5 and 20. . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6 Comparison between forward (left) and E
Y
optimal (right) initial conditions of
ˆ
Y
for R= 5, M = 1, m= 0.05, and τ
0
= 20. . . . . . . . . . . . . . . . . . . . . . . 36
vii
3.7 The E
T
-optimized amplification factors as functions of τ
0
for R = 5, M = 1 and
m= 0.01 (left) and m= 0.05 (right). . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.8 Optimal G
T
(solid) versus their corresponding inviscid G
inv
(dashed) amplification
factors for R= 5, M = 1 and m= 0.01. . . . . . . . . . . . . . . . . . . . . . . . 38
3.9 Optimal amplification factors for m= 0.05, τ
0
= 10 as a function of R with uniform
viscosity M = 1 (left) and as a function of M with density ratio R= 5 (right). . . . 39
3.10 Two-dimensional optimal perturbation fields u
1
(a,d), w
1
(b,e) and Y
1
(c,f) at initial
time τ
0
(a–c) and at time (τ = 31.64) when G
T
≈ 10 (d–f ) of typical growing
disturbance for R= 5, M = 1, m= 0.05 and τ
0
= 10. . . . . . . . . . . . . . . . . 40
3.11 Optimal amplification factors G
T
, G
K
and G
Y
with different cost functions E
Y
(left), E
K
(center) and E
T
(right) for R= 5, M = 1, τ
0
= 10, m= 0.16 and k = 0.05. 41
3.12 E
T
-optimal amplification factors for R = 5, M = 1, τ
0
= 10, and m = 0.16, and
different values of k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.13 Two-dimensional (k= 0) optimal perturbation fields u
1
(a,d), w
1
(b,e) and Y
1
(c,f) at
initial time τ
0
(a–c) and at time (τ = 15.52) when G
T,max
≈ 2.53 (d–f ) for R= 5,
M = 1, m= 0.16 and τ
0
= 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.14 Three-dimensional optimal perturbation fields u
1
, v
1
, w
1
and Y
1
at τ
0
= 10 (left)
and τ = 15.88 (right) for R= 5, M = 1, k= 0.05 and m= 0.16. . . . . . . . . . . 44
3.15 Temporal evolution of the shear-layer width (δ/δ
i
)
p
t
i
/T
DNS
including base pro-
file (left), and perturbation width excluding the base profile(δ
1
/δ
i
)
p
t
i
/T
DNS
(right)
for forward integrations with R = 5, M = 1, τ
0
= 1.12 and several values of m,
compared against DNS result taken from Gat et al. (2017). . . . . . . . . . . . . . 45
viii
4.1 Flow field initializations for two-dimensional DNS with scalar Y (a), velocities
u (b) and w (c) shown for R = 5. Single-mode scalar LOP with m = 0.05 and
ε = 0.01 is imposed at τ = τ
0
= 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2 The evolutions of energy of perturbation scalar Y
1
(a), velocities u
1
(b) and w
1
(c) obtained from DNS (solid) and LSA (dashed), as functions of τ
0
for R = 5,
m= 0.05 and ε = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3 Power spectrum P
n
(x,τ) for R = 5, m = 0.05, τ
0
= 1 and ε = 0.01 at τ− τ
0
=
0,12,19,24,30 and 50, respectively (a–f ). See text for significance of these times. . 62
4.4 The evolutions of energy of perturbation scalar Y
1
(a), velocities u
1
(b) and w
1
(c)
obtained from DNS (solid) and LSA (dashed), as functions of ε for R= 5, m= 0.05
and τ
0
= 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.5 The evolutions of energy of perturbation scalar Y
1
(a), velocities u
1
(b) and w
1
(c) obtained from DNS (solid) and LSA (dashed), as functions of R for m= 0.05,
τ
0
= 10 and ε = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.6 Two-dimensional scalar Y (a–c), velocities u (d–f ) and w (g–i) fields at τ≈ 17.17
(a,d,g), at τ = 28.71 (b,e,h) and at τ = 32.71 (c,f,i) for R = 5, m = 0.05, τ
0
= 10
and ε = 0.01. See text for the significance of these three instances of time. . . . . . 66
4.7 Two-dimensional scalar Y (a–c), velocities u (d–f ) and w (g–i) fields at τ≈ 57.02
(a,d,g), at τ = 65.98 (b,e,h) and at τ = 72.69 (c,f,i) for R = 5, m = 0.05, τ
0
= 50
and ε = 0.01. See text for the significance of these three instances of time. . . . . . 67
4.8 Mixing layer widths δ(x,τ) for R= 5, m= 0.05, ε = 0.01 as functions of τ
0
. . . . 68
4.9 The evolutions of energy of perturbation scalar Y
1
, velocities u
1
and w
1
obtained
from DNS (solid) and LSA (dashed) for R= 5, m= 0.16, τ
0
= 10 and ε = 0.01. . . 69
ix
A.1 Imaginary frequency ω
I
as: (a) function of R for M = 1 and τ
0
= 10, (b) function
of M for R= 5 and τ
0
= 10, and (c) function of τ
0
for R= 5 and M = 1. . . . . . . 82
A.2 The eigenfunctions for R= 5, M = 1, τ
0
= 10, and m= 0.05. . . . . . . . . . . . . 83
x
Abstract
The hydrodynamic stability of a vertical interface separating two fluid columns of different
densities and viscosities, subject to a gravitational acceleration field parallel to the interface, is
investigated. This variable-density flow possesses a time-dependent laminar one-dimensional ref-
erence state, where each column accelerates at different rates owing to their different densities.
This so-called accelerating Kelvin-Helmholtz configuration is studied in the framework of zero-
Mach number limit, which allows scalar and density variations but ignores acoustic effects. First,
the inviscid linear stability analysis (LSA) is carried out with immiscible interface in the presence
of surface tension. The LSA is not amenable to Fourier or Laplace solution in time. Instead, the
equations are derived by the initial value problem (IVP) method and solved analytically. The so-
lution is expressed in terms of the well-known parabolic cylinder function. Second, the viscous
LSA problem where the fluid columns are miscible is considered, while allowing for variation in
viscosity. The interface thickness grows as the square root of time (by diffusion). Numerical inte-
gration of the linear IVP is carried out and discussed in detail as a function of vertical and spanwise
wavenumbers and the flow parameters. Adjoint-based optimization is performed in order to de-
termine linearly optimal perturbations (LOPs), which are initial conditions that lead to maximum
growth of disturbances in finite time. Sensitivity of perturbation growth with respect to initial time,
density and viscosity ratios is investigated. Third, direct numerical simulation (DNS) is carried out
under the uniform viscosity assumption, in order to understand how the LOPs evolve into fully
nonlinear manner. Each simulation run is initialized with the base states being perturbed by single-
mode perturbations with finite amplitude. The evolution of perturbation energy is compared with
Abstract adapted from Prathama & Pantano (2017, 2021).
xi
that of the LSA. As disturbances grow sufficiently large, nonlinear effects become important and
deviations from LSA results are observed. Qualitative description of the flow fields in different
growth stages, as well as evolution of mixing layer growth and energy cascade from single mode
to other modes are discussed. The overall results show that there are two mechanisms at play
in this flow configuration: one involving the growing momentum (kinetic energy) in the vertical
direction, and the other involving diffusional spread (potential energy) in the horizontal direction.
xii
Chapter 1
Introduction
1.1 Background
Hydrodynamic stability theory concerns the study of fluid behaviour in the vicinity of a base state
(Foures et al., 2013). If an initial deviation from the base state decays over time, the flow is
considered to be stable; if an initial deviation increases over time, the flow is considered unstable.
Stability analysis is of central importance in fluid mechanics ever since one realizes that solution
of the governing equations is not always observed in practice (Kerswell, 2018). This theory is
connected with the problem of turbulence since one cause of the appearance of turbulence for any
motion is instability, and for several motions, when the Reynolds number exceeds a certain critical
value, the laminar flow turns into a turbulent one (Georgescu, 1985).
Rayleigh-Taylor instability (RTI) and Kelvin-Helmholtz instability (KHI) arise in natural phe-
nomena as well as in technological applications, such as inertial confinement fusion, geophysics,
astrophysics, compressible turbulence, and combustion (Sharp, 1984; Olson et al., 2011; Gat et al.,
2017), to name a few. Canonical RTI is experienced by an interface between two quiescent fluids
where the heavier fluid is on top of the lighter one, while KHI takes place when there is shear
across the interface of two fluids with different velocities and densities (Rayleigh, 1883; Taylor,
Chapter adapted from Prathama & Pantano (2017, 2021).
1
1950; Chandrasekhar, 1981). In general, classical temporal stability theory has been focusing
on quasi-parallel or parallel shear flows, either bounded (boundary layers, Poiseuille and Couette
flows) or unbounded flows (jets, wakes, and mixing layers), as well as flows in porous media
(Truzzolillo & Cipelletti, 2018). In atmospheric flows, the linear stability of a vertical interface
between two fluids with different densities has been investigated in the case of shallow convec-
tion with or without wind shear (Bretherton, 1987; Mellado, 2017; Kirshbaum & Straub, 2019).
In these analyses, the domain is bounded in the vertical dimension, using homogeneous bound-
ary conditions, and the base flow is assumed to be stationary. The linear stability analysis (LSA)
then leads to time-independent linear differential equations for the flow perturbations which can
be solved as a superposition of Fourier modes (Drazin & Reid, 2004). Long-time instability is
subsequently determined by existence of the least stable eigenmode, with exponential growth rate
indicated by the corresponding eigenvalue. However, it is well-recognized that this approach fails
to take into account the short-term dynamics between the mean flow and perturbation. For in-
stance, the plane Couette flow, which is linearly stable based on eigenvalue analysis alone, expe-
riences large transient algebraic growth prior to asymptotic decay. This is attributed to the fact
that the underlying linear operators are not self-adjoint (Reddy et al., 1993; Trefethen et al., 1993).
Their eigenfunctions are not mutually orthogonal and any initial perturbation constructed from
an eigenfunction expansion may grow substantially (Gustavsson, 1991; Butler & Farrell, 1992;
Reddy & Henningson, 1993). In addition, modal analysis assumes that the base flow is steady or
quasi-steady, commonly known as quasi-steady-state approximation (QSSA), where the base flow
is assumed to change at a lower rate than the perturbations. However, this yields erroneous re-
sults when analyzing flows with time-dependent base states that evolve at comparable time scales
as the perturbations. These problems can be analyzed formally (if not really practically) by the
initial-value problem (IVP) method; e.g., see Richtmyer (1960) or the review by Schmid (2007).
Tan & Homsy (1986) carried out the LSA of an interface between two miscible fluids with
different viscosities in porous media, by utilizing both QSSA and IVP methods. They showed that
the QSSA agrees poorly with results from the IVP at early times because the former ignores the rate
2
of change of base state, which can be comparable to that of the fluctuations. Comparison between
the growth of random noise and the eigenfunction from QSSA as initial conditions indicates that
the latter case is less stable though. Ben et al. (2002) employed the spectral method to investigate
small-amplitude miscible fingering at inception, using self-similar diffusive front as the base state.
They argued that QSSA is not valid at times where the fingers’ amplitudes are much smaller than
their wavelengths and their growth rate is small compared to the unsteady diffusive spreading
rate of the base state. In general, it is known that the growth rate of disturbances is sensitive to
the given initial condition (Tan & Homsy, 1986) in the IVP formulation, and therefore it is not
unique. Hence, it becomes of interest to determine the most amplified initial condition, usually
up to a finite time horizon. Physically, this optimal growth arises from effective energy transfer
from the mean flow to the non-mutually-orthogonal modes (Farrell, 1988). Several techniques such
as eigenfunction expansion (Butler & Farrell, 1992), expansion in periodic functions (Criminale
et al., 2003), singular value decomposition (Rapaka et al., 2008), and adjoint equations have been
employed to determine linearly optimal disturbances. The latter approach is used in this work; see
review by Luchini & Bottaro (2014) for detailed description of this method.
In geophysical flows, Farrell & Moore (1992) used the adjoint of the quasi-geostrophic dy-
namic equations and performed the so-called direct-adjoint looping technique (Kaminski et al.,
2014), an application of the power method. Corbett & Bottaro (2001) analyzed the crossflow
instability in swept boundary layers, and found that the optimal perturbations take the form of
nearly-streamwise vortices at inception and develop into streaks. Similar study on Hiemenz flow
with non-parallel base state at the leading-edge of swept wings was also carried out to determine
both optimal disturbances and its optimal control (Gu´ egan et al., 2006). Later, Doumenc et al.
(2010) and Daniel et al. (2013) investigated the transient growth of perturbations in Rayleigh–
B´ enard–Marangoni convection and density driven fingering, respectively. More recently, Arratia
et al. (2013) and Kaminski et al. (2014) studied the transient linear growth in homogeneous and
parallel stratified shear layers, respectively. They found that at small enough Richardson num-
ber, the transient optimal perturbations are three-dimensional oblique waves due to Orr and lift-up
3
mechanisms for small and large wavenumbers, respectively; the two-dimensional KHI is the dom-
inant mechanism at large times. At large Richardson number case, however, the modal instability
is suppressed, and instead the growth of internal waves outside the vorticity layer is responsible
for short-time amplification. These works were subsequently extended to variable-density mixing
layer experiencing the primary KHI (Lopez-Zazueta et al., 2016). Finally, more studies have been
focusing on transient growth in porous media flows (Hota et al., 2015).
While the linearized governing equations assume infinitesimal perturbations, they no longer
apply when amplitude of perturbations grow large enough and nonlinear effects become impor-
tant. For wavelike disturbances, this means that the Fourier components are not independent but
coupled together (Schmid & Henningson, 2001). For instance, even in strongly stratified layers
with relatively large Richardson number (Ri> 0.3), the ensuing transient growth mechanisms may
allow a transition to (inevitably transient) turbulence in a way that the transition is not mediated by
the growth of a linear instability (Kaminski et al., 2017). Full nonlinear governing equations then
need to be solved in the framework of direct numerical simulations (DNS).
Recently, there have been attempts to find nonlinear optimal perturbations (NLOP) for a given
base flow, with very similar direct-adjoint-looping technique applied to full Navier-Stokes equa-
tions (Kerswell, 2018). For instance, NLOP have been found for unstratified flows such as plane
Couette flow (Rabin et al., 2012; Duguet et al., 2013), pipe flow (Pringle & Kerswell, 2010), and
boundary layers (Cherubini et al., 2015). In general, the structure of the nonlinear optimal per-
turbations closely matches that of the linear optimal perturbations at small initial amplitude, and
becomes increasingly localized in space as the initial perturbation amplitude is increased, until
some critical amplitude at which transition to turbulence is able to occur (Kaminski et al., 2017).
NLOP have been shown to unpack and grow via a sequence of essentially linear transient growth
mechanisms which communicate via nonlinearity as described in Duguet et al. (2013) and Ker-
swell et al. (2014). However, although previous studies indicate that NLOP may attain higher
gains than their linear counterparts for unstratified flows (Schmid, 2007; Rabin et al., 2012; Lu-
chini & Bottaro, 2014), finding NLOP introduces additional computational complexity, due to the
4
dependence of the adjoint equations on the forward variables. Frequent checkpointing is typically
employed to store velocity fields at intermediate points during the forward integration stage, to
later recompute the whole flow fields at all times using interpolation (Kerswell, 2018). Hence,
a thorough parametric search (in terms of wavenumbers, disturbance amplitudes, and target time
of the optimization) is extremely computationally demanding (Pralits et al., 2015). Alternatively,
an investigation into how the linearly optimal disturbances are evolving in fully nonlinear manner
(Kaminski et al., 2014, 2017) can be pursued.
Between the linear stability problem and fully nonlinear simulations, many researchers have
considered the weakly nonlinear stability analysis, employing the techniques such as methods of
multiple scales, amplitude expansion method (Watson, 1960), or center manifold theorem (Fu-
jimura, 1989, 1991). The analysis is performed in the vicinity of neutral stability states, about
which small amplitude assumption holds. This results in the Stuart-Landau equation (Stuart, 1960),
an evolution equation for the amplitude, which is also assumed to be weakly time-dependent. This
amplitude equation is absent in linear stability theory, since the eigenfunctions obtained from LSA
are independent of the amplitude of the disturbance. Furthermore, this assumes that the slowly
varying amplitude is separate from the exponential-like growth, hence the existence of multiple
time scales. The approach has been shown to be successful in flow configurations around a neutral
state. Method of multiple scales, for example, can then be applied in the vicinity of this neu-
trally stable state, giving rise to supercritical or subcritical bifurcations (Yaglom, 2012). However,
in flows where a neutrally stable state does not exist in the sense of linear theory, such as those
with time-dependent base states, the applicability of classical weakly nonlinear analysis is not
warranted.
1.2 Motivation and scope
In the current study, a baroclinically generated shear layer formed between two vertically oriented
columns with different densities and viscosities is considered. The columns are subject to an
5
external gravitational acceleration field, resulting in initially perpendicular pressure and density
gradients. While one is tempted to think of the flow as an RTI configuration, the situation resembles
more a KHI flow of variable density where the velocities of the columns grow linearly in time. This
variable density flow is studied in the framework of zero-Mach-number limit, where the system of
equations allows for scalar and density variations, as well as interaction with the hydrodynamic
flow field. Here, the small-scale acoustic effects are filtered-out, resulting in simpler system that
is less expensive to implement numerically than the fully compressible flow equations (Majda &
Sethian, 1985). The canonical configuration has been recently investigated by Gat et al. (2016,
2017) using DNS. This motivates the stability analysis of this flow, in particular the influence of
initial conditions to the growth of perturbations and the evolution of flow structures at short and
long times.
Initially, the LSA of this configuration is carried out in the inviscid limit with the presence of
surface tension, which is discussed in chapter 2. We derive the equations analytically by the IVP
method and express the solution in terms of the well-known parabolic cylinder function. The an-
alytic solution indicates that the flow is unconditionally unstable at all wave modes, owing to the
constantly accelerating columns that provide an unlimited source of kinetic energy. The growth of
the solution is determined to be quadratic in time exponential. Subsequently, in chapter 3 we dis-
cuss the LSA for miscible fluids with viscosity and diffusion, where the interface thickness grows
as the square root of time (by diffusion), and compare it with the inviscid limit results determined
previously. The base states is discussed and the IVP approach is elaborated. Adjoint-based opti-
mization is then performed in order to determine initial perturbation that yields maximum growth
of disturbances at later times. Different choices of the objective functions are explored to clarify
their impact on the conclusions. Results are presented and the sensitivity of perturbation growth
with respect to initial time, density, and viscosity ratios are discussed. Finally, in chapter 4 we
allow the obtained linearly optimal perturbations (LOPs) with finite initial amplitudes, rather than
the random two-dimensional field with a wide spectrum implemented in DNS study of Gat et al.
6
(2017), to evolve in a fully nonlinear manner, within the framework of direct numerical simula-
tion. An additional parameter, that is the initial energy of perturbations (equivalently, the initial
amplitude of perturbations), may play important role in flow transition and alter the flow structures
as compared to linearized analysis. Results such as perturbation energy evolutions are compared
against those of LOPs in LSA. Qualitative description of the nonlinear flow evolution in different
growth stages, as well as evolution of mixing layer growth and energy cascade from single mode to
other modes are discussed. The influence of varying density ratios and initial time is investigated.
Finally, concluding remarks are made in chapter 5.
7
Chapter 2
Inviscid linear stability analysis
2.1 Problem formulation
We consider a vertical interface which separates two immiscible fluid columns with different densi-
ties. Figure 2.1 illustrates the flow configuration and the Cartesian coordinate system used through-
out the analysis. A constant gravitational acceleration field g is acting downwards (−ˆ z z z direction),
and the unperturbed interface is initially located along the z-axis at x = 0. The densities of the
left and right columns are ρ
2
and ρ
1
, respectively. We define the following parameters for our
subsequent discussions
R≡
ρ
1
ρ
2
≥ 1, ¯ ρ≡
ρ
1
+ ρ
2
2
, (2.1)
where R is the density ratio of high- to low-density fluid and ¯ ρ is the mean density.
The flow obeys the Euler equations at each side of the interface,
ρ

∂u u u
∂t
+ u u u· ∇u u u

=−∇p−(ρ− ¯ ρ)gˆ z z z, (2.2)
∂Y
∂t
+ u u u· ∇Y = 0, (2.3)
Chapter adapted from Prathama & Pantano (2017).
8
F (x, t)
x
z
n
g
r
1
r
2
Y
+

(x, t) Y
-

(x, t)
r
+

(x, t)
p
+

(x, t)
u
+

(x, t) u
-

(x, t)
p
-

(x, t)
r
-

(x, t)
Figure 2.1: Planar sketch of flow configuration in Cartesian coordinates for the inviscid analysis (y
is perpendicular to the page).
∇· u u u= 0, (2.4)
where ρ is the density of the mixture, u u u= (u,v,w) is the velocity vector, p is the pressure, and Y
is the mass fraction of one of the species. There is no chemical reaction between the two streams.
We note that the flow velocities are formulated in the frame of reference where the net vertical
momentum is zero; see Gat et al. (2017). For a binary mixture in the zero Mach number limit
(Sandoval, 1995), the density is related to the local mass fraction, according to
1
ρ(x x x,t)
=

1
ρ
1
−
1
ρ
2

Y(x x x,t)+
1
ρ
2
≡ ν
2
+ Δν Y(x x x,t), (2.5)
where
Δν =

1
ρ
1
−
1
ρ
2

, ν
i
=
1
ρ
i
. (2.6)
2.2 Flow analysis for small perturbations
We now introduce the following expansion for all flow fields
Ψ Ψ Ψ(x x x,t)= Ψ Ψ Ψ
0
(x,t)+ εΨ Ψ Ψ
1
(x x x,t)+O(ε
2
), (2.7)
9
where Ψ Ψ Ψ={u,v,w,Y, p} and ε is assumed to be very small. The subscripts 0 and 1 denote the base
flow and the first order perturbation fields, respectively. The solution for the base flow is given by
u
0
= 0, v
0
= 0, p
0
= 0, (2.8)
∂Y
0
∂t
= 0, ρ
0
=(ν
2
+ Δν Y
0
)
−1
, (2.9)
∂w
0
∂t
=−

1−
¯ ρ
ρ
0

g. (2.10)
From Eq. (2.9), we infer that Y
0
(x,t) = H(x), where H(x) denotes the Heaviside function, and
subsequently from Eq. (2.5) and Eq. (2.10) we obtain
∂w
0
∂t
=−[1− ¯ ρ(ν
2
+ Δν H(x))]g. (2.11)
Denoting± as the right and left fields with respect to the interface, respectively, we have from
Eq. (2.11)
w
±
0
=±
(1− R)
2R
±
gt≡∓2a
±
t, (2.12)
where R
+
= R and R
−
= 1, and
a
±
≡
(R− 1)
4R
±
g. (2.13)
Now, the equations for the first-order perturbation fields are
∂u
1
∂t
+ w
0
(t)
∂u
1
∂z
=−
1
ρ
0
∂ p
1
∂x
, (2.14)
∂v
1
∂t
+ w
0
(t)
∂v
1
∂z
=−
1
ρ
0
∂ p
1
∂y
, (2.15)
10
∂w
1
∂t
+ w
0
(t)
∂w
1
∂z
=−
1
ρ
0
∂ p
1
∂z
+ Δν g ¯ ρY
1
. (2.16)
∂Y
1
∂t
+ w
0
(t)
∂Y
1
∂z
= 0, (2.17)
∂u
1
∂x
+
∂v
1
∂y
+
∂w
1
∂z
= 0, (2.18)
valid both on the left and right of the interface. The fields obey different equations because
w
0
(t)= w
±
0
(t) is different at each side of the interface. We now introduce Fourier modes along the
homogeneous directions, according to
(u
1
,v
1
,w
1
,Y
1
, p
1
)=

ˆ u(x,t), ˆ v(x,t), ˆ w(x,t),
ˆ
Y(x,t), ˆ p(x,t)

e
i(ky+mz)
+ c.c., (2.19)
where k and m are the wave numbers of the disturbances in the y- and z-directions, respectively.
Eqs. (2.14)-(2.18) become
∂ ˆ u
∂t
+ imw
0
(t) ˆ u=−
1
ρ
0
∂ ˆ p
∂x
, (2.20)
∂ ˆ v
∂t
+ imw
0
(t) ˆ v=−
1
ρ
0
ik ˆ p, (2.21)
∂ ˆ w
∂t
+ imw
0
(t) ˆ w=−
1
ρ
0
im ˆ p+ Δν g ¯ ρ
ˆ
Y, (2.22)
∂
ˆ
Y
∂t
+ imw
0
(t)
ˆ
Y = 0, (2.23)
∂ ˆ u
∂x
+ ik ˆ v+ im ˆ w= 0. (2.24)
11
The scalar equation, Eq. (2.23), can be integrated immediately using Eqs. (2.12)-(2.13) to
ˆ
Y
±
(x,t)=
ˆ
A
±
(x)exp

±ima
±
t
2

, (2.25)
and in order for
ˆ
Y to decay to zero at infinity, we must have
ˆ
A
±
(x) = 0. A pressure equation can
be derived by combining Eqs. (2.20)-(2.22) with Eq. (2.24), giving
−
∂
∂x

1
ρ
0
∂ ˆ p
∂x

+
k
2
+ m
2
ρ
0
ˆ p=−imΔν g ¯ ρ
ˆ
Y. (2.26)
Using
ˆ
Y
±
= 0, Eq. (2.26) becomes
−
d
2
ˆ p
±
dx
2
+(k
2
+ m
2
) ˆ p
±
= 0, (2.27)
whose solution is
ˆ p
±
=
ˆ
C
±
(t) e
∓K x
, (2.28)
where
K
2
= k
2
+ m
2
, (2.29)
in which we have used the fact that ˆ p has to remain finite as|x|→ ∞.
Now, the pressure and scalar fields are introduced into the momentum equations Eqs. (2.20)-
(2.22), showing that the general structure of the solutions is of the form
∂
ˆ
ϕ
±
∂t
∓ 2ima
±
t
ˆ
ϕ
±
=
ˆ
C
±
(t), (2.30)
where
ˆ u
±
=±
K
ρ
±
ˆ
ϕ
±
e
∓K x
, ˆ v
±
=−
ik
ρ
±
ˆ
ϕ
±
e
∓K x
, ˆ w
±
=−
im
ρ
±
ˆ
ϕ
±
e
∓K x
. (2.31)
12
Eq. (2.30) has general solution
ˆ
ϕ
±
(t)=
ˆ
B
±
e
±ima
±
t
2
+
Z
t
0
ˆ
C
±
(t
′
)e
±ima
±
(t
2
−t
′2
)
dt
′
, (2.32)
and for zero initial perturbation velocities, we have
ˆ
B
±
= 0, which results finally in
ˆ
ϕ
±
(t)=
Z
t
0
ˆ
C
±
(t
′
)e
±ima
±
(t
2
−t
′2
)
dt
′
. (2.33)
Subsequently, we consider the interface (see figure 2.1), here defined implicitly by
F(x,y,z,t)= x− x
f
(y,z,t)= 0, (2.34)
where subscript f denotes interface quantities. From Eq. (2.34) we can calculate the normal,
according to
n n n=
∇F
|∇F|
=
1
q
1+ x
2
f,y
+ x
2
f,z
(1,−x
f,y
,−x
f,z
). (2.35)
Similarly, the velocity of the interface u
f
is given by
∂F
∂t
+ u
f
n n n· ∇F =
∂F
∂t
+ u
f
|∇F|= 0, (2.36)
which gives the normal velocity
u
f
=−
∂F
∂t
|∇F|
=
x
f,t
q
1+ x
2
f,y
+ x
2
f,z
. (2.37)
The curvature of the interface is given by the standard formula
∇· n n n=
(x
2
f,y
x
f,yy
+ x
2
f,z
x
f,zz
)−(x
f,yy
+ x
f,zz
)(1+ x
2
f,y
+ x
2
f,z
)
(1+ x
2
f,y
+ x
2
f,z
)
3/2
. (2.38)
13
The jump conditions at the interface are
−u
f
[ρ]+[ρu
n
] = 0, (2.39)
−u
f
[ρY]+[ρYu
n
] = 0, (2.40)
−u
f
[ρu
n
]+[ρu
2
n
+ p] = −σ∇· n n n, (2.41)
−u
f
[ρu
t1
]+[ρu
n
u
t1
] = 0, (2.42)
−u
f
[ρu
t2
]+[ρu
n
u
t2
] = 0, (2.43)
where σ is the interfacial surface tension, u
n
, u
t1
, and u
t2
are the normal and tangential veloc-
ity components, respectively, and square brackets denote the jump across the interface. We now
introduce the expansions
x
f
= εζ + O(ε
2
), n n n=(1,−εζ
y
,−εζ
z
)+ O(ε
2
),
t t t
1 1 1
=(εζ
y
,1,0)+ O(ε
2
), t t t
2 2 2
=(εζ
z
,0,1)+ O(ε
2
)
u
f
= εζ
t
+ O(ε
2
), ∇· n n n=−ε(ζ
yy
+ ζ
zz
)+ O(ε
2
)















(2.44)
and
u
n
=(εu
1
,εv
1
,w
0
(t)+ εw
1
)·(1,−εζ
y
,−εζ
z
)= ε(u
1
− w
0
(t)ζ
z
)+ O(ε
2
), (2.45)
u
t
1
=(εu
1
,εv
1
,w
0
(t)+ εw
1
)·(εζ
y
,1,0)= εv
1
+ O(ε
2
), (2.46)
u
t
2
=(εu
1
,εv
1
,w
0
(t)+ εw
1
)·(εζ
z
,0,1)= w
0
(t)+ εw
1
+ O(ε
2
). (2.47)
Introducing Eqs. (2.44)-(2.47) into Eqs. (2.39)-(2.43), expanding all quantities around x = 0,
and gathering terms up to first order in ε, results in
[p
1
]= σ(ζ
yy
+ ζ
zz
), (2.48a)
14
u
±
1
− w
±
0
(t)ζ
z
= ζ
t
. (2.48b)
From Eq. (2.28) and Eq. (2.48a), and introducing Fourier modes, we have
ˆ
C
+
(t)=
ˆ
C
−
(t)− σK
2
ˆ
ζ, (2.49)
and similarly, Eq. (2.48b) becomes
∂
ˆ
ζ
∂t
∓ 2ima
±
t
ˆ
ζ = ˆ u
±
(x= 0)=±
K
ρ
±

Z
t
0
ˆ
C
±
(t
′
)e
±ima
±
(t
2
−t
′2
)
dt
′

. (2.50)
At this point, we have three equations with three unknowns:
ˆ
ζ and
ˆ
C
±
. Multiplying Eq. (2.50)
by e
∓ima
±
t
2
, taking derivative with respect to time, cancelling exponential terms on both sides, and
eliminating
ˆ
C
±
using Eq. (2.49), we obtain
d
2 ˆ
ζ
dt
2
+

β
2
0
−
ω
4
0
4
t
2

ˆ
ζ(t)= 0, (2.51a)
with
β
2
0
=
K
3
σ
ρ
+
+ ρ
−
, ω
4
0
= g
2
m
2
(R− 1)
2
R
.
Eq. (2.51a) can be non-dimensionalized by introducing ˜ τ ≡ ω
0
t, and defining a dimensionless
parameter
s=
β
2
0
ω
2
0
=
p
(k
2
+ m
2
)
3
σ
√
R
(ρ
+
+ ρ
−
)(R− 1)gm
, (2.52)
giving
d
2 ˆ
ζ
d ˜ τ
2
−

1
4
˜ τ
2
− s

ˆ
ζ( ˜ τ)= 0. (2.53)
The general solution of Eq. (2.53) can be expressed as (see Abramowitz & Stegun, 1965)
ˆ
ζ( ˜ τ)=
ˆ
C
1
D
s−1/2
( ˜ τ)+
ˆ
C
2
D
−s−1/2
(i ˜ τ), (2.54)
15
where D
a
(z) denotes the parabolic cylinder function of order a (note that although an imaginary
argument is used in the second term of Eq. (2.54), the complete result is real after the constants
are determined). Enforcing the initial conditions
ˆ
ζ = 1 and
ˆ
ζ
˜ τ
= 0 (obtained by setting t = 0 in
Eq. (2.50)), we get
ˆ
C
1
=−i 2
−s
Γ

−
s
2
+
1
4

Γ

s
2
+
1
4

ˆ
C
2
, (2.55)
and
ˆ
C
2
=
(1+ i)
π
2
−
7
4
−
s
2
e
−
1
2
iπs
(1+ e
2iπs
)Γ

−
s
2
+
3
4

Γ

s+
1
2

, (2.56)
where Γ(z) is the Gamma function.
We will discuss here the general behaviour of
ˆ
ζ( ˜ τ), and postpone to the next section a more
detailed consideration of the asymptotic dependence for large ˜ τ. Figure 2.2 shows
ˆ
ζ(τ) obtained
by plotting Eq. (2.54) for various values of s; the results were also verified by direct numerical
integration of Eq. (2.53). First, it is observed that the interface is unconditionally unstable, with the
most unstable behaviour for s= 0 (a result that is expected given the stabilizing effect that surface
tension has on this type of flows). Second, we observe that
ˆ
ζ decays for the particular values
s
n
= 2n+ 1/2, where n is an integer or zero, but this does not imply that the interface is stable.
For these particular values, one arrives at a situation where the interface amplitude and its velocity
both approach zero asymptotically for large ˜ τ. In this hypothetical single-wavemode arrangement,
the interface has been flattened and no further evolution is possible. Evidently, this behaviour
will not arise if there are more wavemodes present. The root of this unusual behaviour can be
traced back to the time-dependent coefficient of Eq. (2.53), which generates solutions that have no
counterpart in constant-coefficient ordinary differential equations. For example, this behaviour is
not observed in LSA where harmonic solutions exists (normal modes). Third, one observes that
ˆ
ζ is not monotonically increasing (or decreasing). The amplitude magnitude generally shows an
oscillatory behaviour (with a growing envelope). When ˜ τ
2
< 4s, the coefficient multiplying
ˆ
ζ is
negative and the solution is oscillatory like. For large times, ˜ τ
2
> 4s, the growth is exponential
like since the coefficient in the second term of the differential equation is positive. In addition, as
16
2 4 6 8
-3
-2
-1
0
1
2
3
= 0
= 1/2
= 1
= 2
= 5/2
= 4
= 9/2
˜ τ
ˆ
ζ
s
s
s
s
s
s
s
Figure 2.2: Evolution of the interface amplitude
ˆ
ζ for various values of s.
it happens in other shear flows, the most unstable growth happens for the two-dimensional mode
(k= 0), suggesting similarity to Squire’s theorem (Squire, 1933) (although here we are concerned
with a non-harmonic linear stability problem). Finally, the mismatch between the temporal phases
of the fields at both sides of the interface, m(a
+
+ a
−
)t
2
, grows with time quadratically and there
is no uniform phase-velocity translation where the interface appears fixed.
Another quantity of interest is the pressure at the interface, which to leading order is equal to
ˆ
C
±
. Figure 2.3 shows the non-dimensional
ˆ
C
∗−
=
K
ˆ
C
−
( ˜ τ)
ω
2
0
ρ
−
=
√
R
4
"

√
R ˜ τ
2
− 2i

ˆ
ζ( ˜ τ)− 4i ˜ τ
d
ˆ
ζ
d ˜ τ
#
−
d
2 ˆ
ζ
d ˜ τ
2
, (2.57)
and
ˆ
C
∗+
( ˜ τ) =
ˆ
C
∗−
( ˜ τ)− s(R+ 1)
ˆ
ζ inferred from Eq. (2.49), as a function of time for different
values of s and density ratio R= 4/3 (behaviour is similar for other values of R). The left column
shows the magnitude of the complex pressure intensities while the right column shows their phases.
It is observed that the pressure coefficient grows fastest for s= 0, consistent with what is observed
for
ˆ
ζ. The initial values can be determined analytically, giving
ˆ
C
∗−
(0)= s−i
√
R/2 and
ˆ
C
∗+
(0)=
17
−sR− i
√
R/2. Finally, for completeness, figure 2.4 shows the real and imaginary parts of
ˆ
ϕ
+
(
ˆ
ϕ
−
is similar) as a function of time for different values of s at R = 4/3. Evidently, these are also
divergent functions of time.
2.3 Discussion
As discussed earlier, this flow is unconditionally unstable at all wavemodes, even in the presence
of surface tension, because of the ever increasing (accelerating) momentum of the free streams. We
did note that the interface amplitude can decay to zero for particular values of s = s
n
= 2n+ 1/2,
for n= 0,1,2,.... First note that these particular values of s are inconsequential because for given
σ there are always other wavemodes that are unstable. To clarify that this is the case in general,
we can study the asymptotic behaviour of the solution for large ˜ τ, using the known asymptotic
expansions of the parabolic cylinder functions (see Abramowitz & Stegun, 1965), which results in
ˆ
ζ( ˜ τ) ≈
ˆ
C
1
e
−
1
4
˜ τ
2
˜ τ
s−
1
2
"
1−
(−s+
1
2
)(−s+
3
2
)
2 ˜ τ
2
+
(−s+
1
2
)(−s+
3
2
)(−s+
5
2
)(−s+
7
2
)
2× 4 ˜ τ
4
−...
#
+
ˆ
C
2
e
1
4
˜ τ
2
(i ˜ τ)
−s−
1
2
"
1+
(s+
1
2
)(s+
3
2
)
2 ˜ τ
2
+
(s+
1
2
)(s+
3
2
)(s+
5
2
)(s+
7
2
)
2× 4 ˜ τ
4
+...
#
. (2.58)
Potential growth in the solution is only possible through the quadratic exponential term in the
second line of Eq. (2.58) that is multiplied by
ˆ
C
2
. Therefore, decay of the interface amplitude is
only possible if
ˆ
C
2
= 0. From Eq. (2.56), we see that this can only happen when e
i2πs
=−1, which
explains the condition s
n
= 2n+ 1/2. We also observe from the asymptotic behaviour that the
dominant growth in time of
ˆ
ζ is of the form
ˆ
ζ∼ e
1
4
˜ τ
2
˜ τ
−s−
1
2
. (2.59)
18
0 2 4 6 8
0
2
4
6
8
10
˜ τ

ˆ
C
∗−

(a)
2 4 6 8
-2
0
2
4
6
8
˜ τ
arg(
ˆ
C
∗−
)
(b)
0 2 4 6 8
0
2
4
6
8
10
˜ τ

ˆ
C
∗+

(c)
0 2 4 6 8
-12
-10
-8
-6
-4
-2
0
˜ τ
arg(
ˆ
C
∗+
)
(d)
Figure 2.3: Magnitudes (a,c) and phases (b,d) of
ˆ
C
∗−
and
ˆ
C
∗+
for density ratio R= 4/3 for various
values of s; legend as in figure 2.2.
19
2 4 6 8
-6
-4
-2
0
2
4
6
˜ τ
Re(
ˆ
ϕ
+
)
(a)
2 4 6 8
-6
-4
-2
0
2
4
6
˜ τ
Im(
ˆ
ϕ
+
)
(b)
Figure 2.4: Real and imaginary values of
ˆ
ϕ
+
for density ratio R = 4/3 and various values of s;
legend as in figure 2.2.
20
Chapter 3
Viscous linear stability analysis
3.1 Problem formulation
We consider the linear stability of a vertical interface separating two miscible fluid columns of
different densities and viscosities under the influence of gravity. Figure 3.1 illustrates the flow
configuration, where the interface is depicted as one that experiences viscous spreading. In addition
to density-related parameters defined in Eq. (2.1), we also denote the viscosity of the right and left
fluid columns as μ
1
and μ
2
, respectively. The following constant parameters are introduced:
M≡
μ
1
μ
2
, ¯ μ≡
μ
1
+ μ
2
2
, (3.1)
where M and ¯ μ are the dynamic viscosity ratio and the average viscosity, respectively.
The equations for conservation of linear momentum, species concentration and mass are given
by
ρ

∂u u u
∂t
+ u u u· ∇u u u

=−∇ ˜ p−(ρ− ¯ ρ)gˆ z z z+ ∇·

μ

∇u u u+ ∇u u u
T

, (3.2)
ρ

∂Y
∂t
+ u u u· ∇Y

= ∇·(ρD∇Y), (3.3)
Chapter adapted from Prathama & Pantano (2021).
21
x
z
g
r
1
r
2
Y
+

(x, t) Y
-

(x, t)
r
+

(x, t)
p
+

(x, t)
u
+

(x, t) u
-

(x, t)
p
-

(x, t)
r
-

(x, t)
m
-

(x, t) m
+

(x, t)
m
2
m
1
Figure 3.1: Planar sketch of flow configuration in Cartesian coordinates (y is perpendicular to the
page).
∂ρ
∂t
+ ∇·(ρu u u)= 0, (3.4)
where ρ is the density of the mixture, u u u=(u,v,w) is the velocity vector, ˜ p= p−
2
3
μ(∇· u u u) is the
modified pressure, μ is the dynamic viscosity of the mixture, and D is the binary species diffusivity.
The other variables are as defined in §2.1, and the density-local mass fraction relationship is given
by Eq. (2.5). The dynamic viscosity is assumed to vary linearly with the local mass fraction
according to
μ(x x x,t)=(μ
1
− μ
2
)Y(x x x,t)+ μ
2
. (3.5)
Combining Eqs. (3.3)-(3.4), while using Eq. (2.5), one obtains the mass constraint equation (Lu &
Pantano, 2020), given by
∇· u u u=−∇·

D
ρ
∇ρ

. (3.6)
Now, we introduce the following non-dimensionalization:
x x x
∗
=
x x x
L
, τ
∗
=
t
T
, u u u
∗
=
u u u
U
, ˜ p
∗
=
˜ p
¯ ρU
2
, ρ
∗
=
ρ
¯ ρ
, μ
∗
=
μ
¯ μ
, (3.7a)
22
where the length, time and velocity scales of the mixing region are defined as
L=

¯ ν
2
g

1/3
, T =

¯ ν
g
2

1/3
, U =
L
T
=( ¯ νg)
1/3
, (3.7b)
with ¯ ν = ¯ μ/ ¯ ρ the average kinematic viscosity. The Reynolds number is unity with these scalings;
another parameter of relevance is the dimensionless time τ, which is elaborated in §3.2. The
dimensionless governing equations become (dropping the asterisk and tilde notations)
∂u u u
∂τ
+ u u u· ∇u u u=−
∇p
ρ
−

1−
1
ρ

ˆ z z z+
∇·

μ

∇u u u+ ∇u u u
T

ρ
, (3.8)
∂Y
∂τ
+ u u u· ∇Y =
∇·(μ∇Y)
Sc ρ
, (3.9)
∇· u u u=
1
Sc
∇·

μ∇
1
ρ

, (3.10)
where Sc=(μ/ρ)/D is the Schmidt number, i.e., the ratio of the viscous to the species diffusivity.
Furthermore, Eq. (2.5) and Eq. (3.5) give the non-dimensional density and viscosity, respectively,
according to
1
ρ
=
R+ 1
2R
[R−(R− 1)Y], (3.11a)
and
μ =
2
M+ 1
[1+(M− 1)Y], (3.11b)
respectively.
23
3.2 Flow analysis for small perturbations
We introduce the following expansion for all flow fields in Eqs. (3.8)-(3.10):
Ψ Ψ Ψ(x x x,τ)= Ψ Ψ Ψ
0
(x,τ)+ ε Ψ Ψ Ψ
1
(x x x,τ)+O(ε
2
), (3.12)
where Ψ Ψ Ψ=[u,v,w,Y, p] and ε is assumed to be very small. The subscripts 0 and 1 denote the base
flow and the first-order perturbation fields, respectively. The latter are the main subject of this
paper.
3.2.1 Base flow
The base state for this configuration is a function of x and τ only. There are no characteristic
length and time scales in the flow, as can be seen from the scalings in Eq. (3.7b), which implies
that the base flow is invariant under the stretching transformation x→ ax, τ→ a
2
τ. This lowers
the number of independent parameters from two to one, reducing the partial differential equations
into ordinary ones (Ben et al., 2002). Therefore, it is expected that the base state will not depend
on x and τ separately, since there are insufficient dimensional parameters in this problem (see the
discussion on unsteady unidirectional flow in §4.3 of Batchelor (1967)). Therefore, the base flow
admits a self-similar solution, with the self-similar variable η defined as
η≡
x
√
2Kτ
, where K≡
R+ 1
RSc
. (3.13)
Note that the self-similar solutions apply outside the singular instant corresponding to τ = 0. This
behaviour is typical in problems with similarity solution, e.g., Blasius boundary layer and unsteady
unidirectional flow in unbounded domain (Batchelor, 1967). Introducing the self-similar variable
24
into the governing equations where subscript ‘0’ denotes the base flow fields, and using the auxil-
iary functions U
0
, V
0
, W
0
, Y
0
and P
0
, according to
u
0
(x,τ) = −
(R− 1)
2
r
K
2τ
U
0
(η), (3.14)
v
0
(x,τ) = V
0
(η), (3.15)
w
0
(x,τ) = τ W
0
(η), (3.16)
Y
0
(x,τ) = Y
0
(η), (3.17)
p
0
(x,τ) = −
(R− 1)(2Sc− 1)
4Sc τ
P
0
(η), (3.18)
one obtains the following system of equations:
U
0
= μ
0
Y
′
0
, (3.19)
Z
0
(μ
0
V
′
0
)
′
+
1
Sc

(R− 1)μ
0
Y
′
0
+ 2η

V
′
0
= 0, (3.20)
Z
0
(μ
0
W
′
0
)
′
+
1
Sc

(R− 1)μ
0
Y
′
0
+ 2η

W
′
0
−
4
Sc
(W
0
+ 1)+ 2KZ
0
= 0, (3.21)
Z
0
(μ
0
Y
′
0
)
′
+[(R− 1)μ
0
Y
′
0
+ 2η]Y
′
0
= 0, (3.22)
P
0
= μ
0
(μ
0
Y
′′
0
+ βY
′
2
0
), (3.23)
with primes denoting derivatives with respect to η, and
μ
0
(η)≡ 1+ β

Y
0
(η)−
1
2

, Z
0
(η)≡ R−(R− 1)Y
0
(η), β ≡ 2

M− 1
M+ 1

. (3.24)
25
The boundary conditions for the scalar are Y
0
(η→ ∞)→ 1 and Y
0
(η→−∞)→ 0 at the regions
of heavier and lighter fluids, respectively. All other quantities, except w
0
, approach zero far from
the interface. The spatial discretization is carried out using a spectral collocation method based
on Chebyshev polynomials, which is realized using the open-source package Chebfun (Aurentz &
Trefethen, 2017; Driscoll et al., 2014). The nonlinear boundary value problems are then solved
with Newton’s method. Figure 3.2 shows the solutions of the base flow fields as a function of η for
several values of density ratio R and uniform viscosity (M = 1). One can observe the asymmetry
of the Y
0
(η) profiles, where the heavier fluid is encroaching more into the lighter fluid region as R
increases. However, the horizontal velocity U
0
and pressure P
0
are more concentrated with higher
peaks for lower R and more spread with lower peaks at higher R. Finally, the vertical velocity
W
0
profiles show larger relative velocities between the two columns as R increases. Note that both
fluids move in the direction of the external uniform acceleration field in an inertial reference frame,
while here one adopts a frame moving downwards with the average velocity of the columns (see
Gat et al. (2017) for details). Figure 3.3 shows the effect of varying the viscosity ratio M, while
keeping the density ratio constant at R = 5. Comparison with the previous case indicates that the
influence of non-uniform viscosity plays a weaker role in the base flow behaviour than that of the
density ratio.
3.2.2 Initial-value problem
We adopt the following Fourier decomposition of the first-order perturbation fields along the ho-
mogeneous directions:
(u
1
,v
1
,w
1
,Y
1
, p
1
)=

ˆ
U(η,τ),
ˆ
V(η,τ),
ˆ
W(η,τ),
ˆ
Y(η,τ),
ˆ
P(η,τ)

e
i(ky+mz)
+ c.c., (3.25)
which transforms the problem into an IVP of the form
∂
ˆ
X X X
∂τ
= L L L(Sc,R,M,k,m,τ)
ˆ
X X X, (3.26)
26
-10 -5 0 5 10
0
0.25
0.5
η
U
0
-5 0 5 10
-1
0
1
2
3
4
5
R = 2.5
R = 5.0
R = 7.5
R = 10.0
η
W
0
-10 -5 0 5 10
0
0.5
1
η
Y
0
-5 0 5 10
-0.4
-0.2
0
0.2
0.4
η
P
0
Figure 3.2: Base flow fields as a function of density ratio R at uniform viscosity M = 1.
27
-10 -5 0 5 10
0
0.2
0.4
η
U
0
-5 0 5 10
-0.5
0
0.5
1
1.5
2 M = 1/4
M = 1/2
M = 1
M = 2
M = 4
η
W
0
-10 -5 0 5 10
0
0.5
1
η
Y
0
-5 0 5 10
-0.3
-0.2
-0.1
0
0.1
0.2
η
P
0
Figure 3.3: Base flow fields as a function of viscosity ratio M at R= 5.
28
where
ˆ
X X X =

ˆ
U,
ˆ
V,
ˆ
W,
ˆ
Y

. The detailed expressions of Eq. (3.26) and the operator L L L are documented
in Appendix A. Homogeneous boundary conditions are imposed at far fields, η →±∞, as we
expect the perturbation fields to be localized in the vicinity of the mixing region. In this study,
the equations are written in the frame of the self-similar coordinate, η. From this point of view,
the base profiles look stationary, and time dependence is transferred to the governing equation,
Eq. (3.26), as an explicit dependence of the linear operator on τ. It is equally possible to formulate
the problem in the x,τ plane and write the base profiles as functions of time. The former approach
is used here.
Since the base state is time-dependent, the time at which the perturbation is introduced, τ = τ
0
,
acts as a parameter and influences the fate of the perturbation evolution. Furthermore, the initial
conditions of the perturbation fields (at τ
0
) are chosen to correspond to the eigenmodes of an
auxiliary QSSA calculation (Tan & Homsy, 1986). This helps us chose appropriate values of
the wavenumber m, for fixed R and M, in the forward integration of the equations. The initial
conditions are immaterial in the optimized solutions, since the form of the initial condition is
selected by maximization of a cost function. Otherwise, QSSA is not investigated here in any
detail since we carry out full time-dependent integration of the linearized stability equations.
For the numerical implementation, the same spatial discretization is used for both base states
and perturbation fields (as described in §3.2.1), where we specify large enough domain size such
that the perturbation fields decay sufficiently at both η =±d, in order to ensure that the bound-
ary conditions are satisfied. The implicit second-order Crank-Nicolson time-marching method is
employed for this linearized problem. However, as is typical in variable-density flows in the zero-
Mach-number limit, there is no evolution equation for pressure. Hence, we use temporal staggered
grid, where the velocity and scalar fields are defined at full time steps n and n+ 1, while the pres-
sure is assumed to be defined at the intermediate (staggered) time step n+ 1/2 (Lu & Pantano,
2020). The conservation of mass is then satisfied at the next time level n+ 1. To this end, the
Schur complement reduction is applied to solve for the pressure field first, followed by velocity
and scalar fields. Details of the implementation is described in Appendix A.
29
Previous studies (Tan & Homsy, 1986; Doumenc et al., 2010; Daniel et al., 2013) indicate that
the maximum amplification gain of perturbation at a certain time τ
1
depends on the definition of
the disturbance energy, which in general is given by
E(τ
1
)=
∞
Z
−∞

A
1

ˆ
U(η,τ
1
)
2
+
ˆ
V(η,τ
1
)
2
+
ˆ
W(η,τ
1
)
2

+ A
2
ˆ
Y(η,τ
1
)
2

dη, (3.27)
where the values of constants A
1
and A
2
control the choice of objective function, i.e., whether the
quantity to be maximized is the kinetic energy, potential energy or total (kinetic and potential)
energy. The normalized amplification factor is defined by
G(τ
0
;τ
1
)= ln

E(τ
1
)
E(τ
0
)

1/2
, (3.28)
where a subscript is attached to E and G to denote which particular cost function/disturbance
energy is used, i.e., K (A
1
= 1,A
2
= 0), Y (A
1
= 0,A
2
= 1) and T (A
1
= 1,A
2
= 1) to denote
kinetic, potential and total, respectively. Note that E
T
= E
K
+ E
Y
. The amplification gain can
be maximized using an adjoint optimization scheme, in which E(τ
1
) is maximized subject to the
constraint E(τ
0
) = 1 (Schmid & Henningson, 2001; Schmid, 2007). The augmented Lagrangian
L of E(τ
1
) is then defined by
L

λ,
ˆ
Ψ Ψ Ψ,
ˆ
Ψ Ψ Ψ
+

= E(τ
1
)− λ[E(τ
0
)− 1]−
τ
1
Z
τ
0
∞
Z
−∞
ˆ
Ψ Ψ Ψ
+
M

ˆ
Ψ Ψ Ψ

dη dτ, (3.29)
where λ and
ˆ
Ψ Ψ Ψ
+
are the Lagrange multipliers, and
ˆ
Ψ Ψ Ψ
+
=

ˆ
U
+
,
ˆ
V
+
,
ˆ
W
+
,
ˆ
Y
+
,
ˆ
P
+

, referred to here
as the adjoint variables. The constraintM

ˆ
Ψ Ψ Ψ

is given by the perturbation Eq. (3.26). In order to
determine first-order optimal conditions, the first variation of the Lagrangian,
δL

λ,
ˆ
Ψ Ψ Ψ,
ˆ
Ψ Ψ Ψ
+

= δE(τ
1
)− λ δE(τ
0
)−
τ
1
Z
τ
0
∞
Z
−∞
δ
ˆ
Ψ Ψ Ψ
+
M

ˆ
Ψ Ψ Ψ

dη dτ−
τ
1
Z
τ
0
∞
Z
−∞
ˆ
Ψ Ψ Ψ
+
δM

ˆ
Ψ Ψ Ψ

dη dτ,
(3.30)
30
must vanish. Some boundary integrals emerge when calculating the partial derivatives ofL , re-
sulting in
τ
1
Z
τ
0

Scμ
0
Z
0
4τ

2
∂
ˆ
U
∂η
ˆ
U
+
+
∂
ˆ
V
∂η
ˆ
V
+
+
∂
ˆ
W
∂η
ˆ
W
+
+
1
Sc
∂
ˆ
Y
∂η
ˆ
Y
+

−
BZ
0
√
τ
ˆ
P
ˆ
U
+
+
(R− 1)μ
0
4τ
∂
ˆ
Y
∂η
ˆ
P
+

+∞
−∞
dτ = 0, where B≡ Sc
p
K/8, (3.31)
which is satisfied by imposing homogeneous boundary conditions as η→±∞. Finally, the opti-
mality conditions are met if the following coupling conditions are satisfied:
2

A
1

ˆ
Uδ
ˆ
U +
ˆ
V δ
ˆ
V +
ˆ
Wδ
ˆ
W

+ A
2
ˆ
Y δ
ˆ
Y

τ
1
=

A
1

ˆ
U
+
δ
ˆ
U +
ˆ
V
+
δ
ˆ
V +
ˆ
W
+
δ
ˆ
W

+ A
2
ˆ
Y
+
δ
ˆ
Y

τ
1
,
(3.32a)
2λ

A
1

ˆ
Uδ
ˆ
U+
ˆ
V δ
ˆ
V +
ˆ
W δ
ˆ
W

+ A
2
ˆ
Y δ
ˆ
Y

τ
0
=

A
1

ˆ
U
+
δ
ˆ
U +
ˆ
V
+
δ
ˆ
V +
ˆ
W
+
δ
ˆ
W

+ A
2
ˆ
Y
+
δ
ˆ
Y

τ
0
.
(3.32b)
Specifically, when maximizing the total energy, the coupling conditions are given by
2
ˆ
X X X

τ
1
=
ˆ
X X X
+

τ
1
, 2λ
ˆ
X X X

τ
0
=
ˆ
X X X
+

τ
0
, (3.33)
where
ˆ
X X X
+
=

ˆ
U
+
,
ˆ
V
+
,
ˆ
W
+
,
ˆ
Y
+

. On the other hand, when the kinetic (potential) energy is to be
maximized, we have the coupling conditions only for velocities (scalar) while the scalar (velocities)
are set to zero at τ
0
and τ
1
for both direct and adjoint variables. In the end, this leads to the adjoint
optimization problem of the form
∂
ˆ
X X X
+
∂τ
+
= L L L
+

Sc,R,M,k,m,τ
+

ˆ
X X X
+
. (3.34)
The adjoint equations are to be integrated backwards in time. The detailed expressions of Eq. (3.34)
and the operator L L L
+
are documented in Appendix A. First, we perform forward (direct) integration
31
in time from τ
0
to τ
1
, utilizing the normalized QSSA eigenmodes as our initial conditions. In
order to prevent too large a round-off error, we determine that τ
1
is the time at which G≈ 7.
At τ
1
, we integrate the adjoint governing equations backwards in time from τ
1
to τ
0
with the
homogeneous boundary conditions. The initial conditions for adjoint integration depend on the
coupling conditions, Eq. (3.32). The numerical implementation is identical to that of forward
integration, but using adjoint equations and initial and boundary conditions instead. The above
procedure is iterative, and we determine convergence between two iterations as when|G
n
(τ
1
)−
G
n−1
(τ
1
)|≤ 10
−6
.
3.3 Results
The inviscid analysis in §2 shows that the interface was unconditionally unstable, even in the
presence of surface tension. The inviscid interface amplitude
˜
ζ( ˜ τ) as a function of the inviscid
time ˜ τ by
˜
ζ( ˜ τ)= C
1
D
−1/2
( ˜ τ)+C
2
D
−1/2
(i ˜ τ)∼ e
˜ τ
2
/4
˜ τ
−1/2
(for large ˜ τ), (3.35)
where D
j
(z) is the parabolic cylinder function of order j and C
1
and C
2
are constants shown in
Appendix B and correspond to an interface that starts its evolution at τ = τ
0
when Δw
0
6= 0; the
values of the constants reported in §2 correspond to the initial condition at τ = 0, which is not
relevant here. The inviscid time is related to the non-dimensionalized time τ according to
˜ τ =

m(R− 1)
√
R

1/2
τ, (3.36)
where m is the non-dimensional wavenumber defined by the scaling in Eq. (3.7b). The normalized
amplification factor of the inviscid case is given by
G
inv
( ˜ τ)= ln

E
Y
( ˜ τ)
E
Y
( ˜ τ
0
)

1/2
=
1
2
ln
h
˜
ζ( ˜ τ)
i
, (3.37)
32
with
˜
ζ( ˜ τ
0
)= 1 as the initial condition. The inviscid asymptotic growth, using Eq. (3.35), behaves
as G
inv
∼( ˜ τ− ˜ τ
0
)
2
/8=(m(R−1)/8
√
R)(τ−τ
0
)
2
. In practice, this amplification factor is so large
that usually one does not reach this large-time regime before the LSA becomes invalid. Further-
more, the inviscid analysis also showed that purely two-dimensional perturbations (k= 0) are more
unstable, for any given vertical wavenumber m, than at any other three-dimensional perturbation
that has k6= 0. In subsequent discussions, the Schmidt number is unity in all cases, which is appro-
priate for understanding mixing of two typical gaseous fluids. But it is acknowledged that this is
less so for liquids where each column may have substantially different diffusivities; a modification
of the present analysis is required for this situation.
Contrary to the inviscid case where all wave modes are unstable, the behaviour of the viscous
case is expected to depend on the value of m. The QSSA approximation, which has limited accu-
racy in this strongly time-dependent problem, is used here to determine meaningful values of m for
some interface parameters of interest. First, we choose R = 5 and M = 1 to focus on a particular
case of reasonable interest of mixing of gases with contrasting densities; e.g. an air-SF
6
interface
(commonly used in experiments) has R= 5.036 at 1 atm. Furthermore, we can interpret that τ
0
is
related to the so-called early-time Reynolds number, Re
0
, which describes the instantaneous base
state condition when the flow is perturbed by infinitesimal disturbances, and can be expressed as
Re
0
=
Δw
0
δ
0
¯ ν
≈
gt
0
√
¯ νt
0
¯ ν
=
g
√
¯ ν
(T τ
0
)
3/2
=
g
√
¯ ν

¯ ν
1/3
g
2/3
τ
0
!
3/2
= τ
3/2
0
, (3.38)
where the velocity difference between the two streams is Δw
0
∼ gt
0
and the thickness of the inter-
face is δ
0
∼
√
¯ νt
0
. Therefore, initially we choose to concentrate on results for τ
0
= 10 correspond-
ing to Re
0
∼ 32, to explore the dynamics of the perturbations and perform a subsequent parametric
study to investigate results for a wider range of τ
0
. The wavenumber of peak amplification accord-
ing to the QSSA is m
peak
= 0.625 and the neutral wavenumber is m
neutral
= 0.148 for τ
0
= 10 from
the data shown in figure A.1. Consequently, we first focus on evaluating the impact of the energy
33
10 20 30
-5
0
5
10
15
τ− τ
0
G
10 20 30
-5
0
5
10
15
τ− τ
0
G
10 20 30
-5
0
5
10
15
G
T
G
K
G
Y
τ− τ
0
G
Figure 3.4: Amplification factors G
T
, G
K
and G
Y
with different initializations of E
Y
(τ
0
)= 1 (left),
E
K
(τ
0
)= 1 (center) and E
T
(τ
0
)= 1 (right), respectively, for R= 5, M = 1, τ
0
= 10, m= 0.05 and
k= 0 (two-dimensional).
function choice, Eq. (3.27), used to measure amplification of the perturbations at a wavenumber
close to m
peak
, i.e m= 0.05 (easy to recall).
3.3.1 Two-dimensional perturbations
First, we consider moderate wavenumbers m and purely two-dimensional perturbations, i.e. k= 0,
because the numerical solutions indicate that these perturbations are more unstable than three-
dimensional perturbations, i.e. k6= 0. Three-dimensional perturbations are investigated in the
next section for conditions close to the neutral amplification regime. Figure 3.4 shows forward
integration (non-optimized) amplification factors for R = 5, M = 1, τ
0
= 10, at the wavenumber
m = 0.05. The initial conditions were chosen to be:
ˆ
U =
ˆ
W = 0 and
ˆ
Y =
ˆ
Y
QSSA
(left),
ˆ
U =
ˆ
U
QSSA
,
ˆ
W =
ˆ
W
QSSA
and
ˆ
Y = 0 (center), and both velocity and scalar eigenfunctions from QSSA
(right), respectively. It is seen that G
K
and G
T
are almost indistinguishable because the amount
of energy in the scalar Y is typically very small, in comparison with the kinetic energy. But, note
that all curves (after some initial transient) evolve similarly as time grows. It is worth mentioning
that figure 3.4(left) shows that the scalar energy E
Y
can decay initially when the flow velocity is
initialized to zero. This is caused by diffusional decay of the scalar perturbations since there is
little flow at early times. The decay of solutions with time was indeed observed in the DNS of
Gat et al. (2017), which manifests early on because pure convective processes start driving the
34
4 6 8 10
-5
0
5
10
15
t
0
= 20
t
0
= 5
t
0
= 1
(Optimal)
(Forward)
}
}
}
˜ τ− ˜ τ
0
G
Y
Figure 3.5: The E
Y
-optimal (solid) versus forward (dashed) amplification factors for R= 5, M = 1,
m= 0.05 and τ
0
= 1, 5 and 20.
instability only at later times. This can be seen by expanding Eq. (3.35) (inviscid behaviour) for
small ˜ τ, resulting in
ˆ
ζ( ˜ τ)= 1+
1
48
˜ τ
4
+ O( ˜ τ
8
), (3.39)
which indicates that there is substantial time, up to ˜ τ ∼ 48
1/4
≈ 2.6, for viscous effects to be
felt. Note that the minimum of G
Y
where perturbations start to grow again happens at ˜ τ∼ 2; see
figure 3.5, where we plot as a function of ˜ τ both non-optimized and E
Y
-optimized solutions for
R= 5, M = 1, m= 0.05 and three values of τ
0
: 1, 5 and 20.
From figure 3.5, we can see that optimal perturbations which maximize E
Y
tend to choose initial
perturbations that minimize this initial decay. This feature disappears, though, when the velocity
is also initialized to be non-zero, as shown in the center and right figures of figure 3.4. Figure 3.6
displays the optimized initial conditions as compared with those from QSSA, for τ
0
= 20. It is
evident that the optimal solutions exhibit a shortened initial diffusional dampening phase and that
the optimal scalar profiles seem to be more compressed, which are better at extracting energy from
the mean shear (Foures et al., 2014).
35
-5 0 5 10
-1
-0.5
0
0.5
1
η
ˆ
Y
-5 0 5 10
-1
-0.5
0
0.5
1
Real
Imag.
Abs
η
ˆ
Y
Figure 3.6: Comparison between forward (left) and E
Y
optimal (right) initial conditions of
ˆ
Y for
R= 5, M = 1, m= 0.05, and τ
0
= 20.
0 10 20 30 40
0
5
10
15
τ− τ
0
G
T
τ
0
= 400
τ
0
= 250
τ
0
= 100
τ
0
= 25
τ
0
= 10
0 10 20 30 40
0
5
10
15
τ− τ
0
G
T
τ
0
= 50
τ
0
= 35 τ
0
= 20
τ
0
= 10
τ
0
= 1
Figure 3.7: The E
T
-optimized amplification factors as functions of τ
0
for R = 5, M = 1 and m =
0.01 (left) and m= 0.05 (right).
36
Subsequently, we concentrate on the optimal growth rates that maximize the total energy func-
tion, E
T
. Figure 3.7 shows the optimized amplification factor G
T
for R= 5, M= 1 at wavenumbers
m = 0.01 and 0.05. At m= 0.01, figure 3.7(left), the amplification factor increases with growing
τ
0
, as expected based on the relationship between τ
0
and the early Reynolds scaling, Eq. (3.38). In-
terestingly, at the higher wavenumber m= 0.05, figure 3.7(right), the growth rate is not monotonic
with increasing τ
0
, where solutions reverse half way from growth to decay. This is not observed
in the QSSA approximation since ω
I
> 0 for all τ
0
cases considered here (see figure A.1(right)),
suggesting that the spreading of the base profile with time is sufficiently strong to modify the qual-
itative behaviour substantially. Further inspection into the perturbation equations for small τ
0
(see
Appendix A) suggests that: (i) initially, the dominant terms are the ones involving Y
0
(equivalently,
Z
0
) and its derivatives (potential energy), yet (ii) at later times, the terms involving W
0
(shear)
and viscous dissipation (approximately proportional to m
2
) are dominant, although these effects
increase slowly.
Figure 3.8 shows a comparison between the optimal amplification factors (solid lines) at small
wavenumber, m= 0.01, and their inviscid counterparts (dashed lines) for large τ
0
plotted in terms
of the inviscid time, Eq. (3.36). The inviscid curves are given by using the exact Eq. (3.35) in
Eq. (3.37). We can observe that the viscous solutions become closer to the inviscid result as
τ
0
grows, corresponding to larger Re
0
, Eq. (3.38), suggesting that we retrieve the inviscid result
asymptotically.
Finally, figure 3.9 (left) shows the optimal amplification factors as a function of density ratio
at M = 1, m = 0.05 and τ
0
= 10. It can be observed that the least stable case is R = 10 since
high density ratio relates to a stronger potential energy of the columns that feed into the instability,
also promoting growth of the mixing layer between the two columns. Figure 3.9 (right) exhibits
the optimal amplification factors as a function of M at R = 5, m = 0.05 and τ
0
= 10. The early-
time behaviours are indistinguishable for all cases, but as viscous effect becomes dominant at later
times, we can see that the flow is least stable when the displacing (heavier) fluid is less viscous
37
0 1 2 3 4
0
5
10
15
˜ τ− ˜ τ
0
G
τ
0
= 400
τ
0
= 250
τ
0
= 100
Figure 3.8: Optimal G
T
(solid) versus their corresponding inviscid G
inv
(dashed) amplification
factors for R= 5, M = 1 and m= 0.01.
than the displaced fluid, i.e., M = 1/4. This behaviour is also observed in other studies of variable-
viscosity miscible fluids (Riaz et al., 2004; Kim & Choi, 2011; Etrati et al., 2018), where significant
flow destabilization is a consequence of the peak of the eigenfunctions tending to shift towards the
less viscous fluid region. Hence the driving kinetic energy of the heavier fluid is less damped and
this unfavourable viscosity difference renders the flow more unstable.
This section concludes with a brief description of the physical features of the perturbation
fields. Figure 3.10 displays the optimal perturbation fields of a typical growing disturbance for
R= 5, M = 1, m = 0.05 and τ
0
= 10, at initial time τ
0
and at the time (τ≈ 31.6) when G
T
≈ 10.
The perturbation fields are characterized by flow patterns that are tilted against the base-flow shear,
which persists throughout the linear amplification region considered here. This misalignment fa-
cilitates efficient energy transfer from the base flow to the perturbation (Farrell, 1988; Schmid &
Henningson, 2001; Foures et al., 2014). As time increases, the peak of u
1
shifts from the heavier
(right) to the lighter (left) fluid region. Furthermore, the peaks of w
1
seem no longer concentrated
around η = 0, but have developed more local peaks adjacent at both sides of the core mixing layer.
38
10 20 30
0
5
10
15
R = 2.5
R = 5
R = 7.5
R = 10
τ− τ
0
G
T
10 20 30
-5
0
5
10
15
M = 1/4
M = 1/2
M = 1
M = 2
M = 4
τ− τ
0
G
T
Figure 3.9: Optimal amplification factors for m = 0.05, τ
0
= 10 as a function of R with uniform
viscosity M = 1 (left) and as a function of M with density ratio R= 5 (right).
Finally, both w
1
and scalar fields have changed their orientation and are now tilted opposite to their
original alignments.
3.3.2 Three-dimensional perturbations
In this section we consider conditions close to the neutral amplification wavenumber where the
perturbations are not expected to grow. The QSSA indicates that this happens at m
neutral
= 0.148
for R = 5, M = 1 and τ
0
= 10. Therefore, we choose m = 0.16, slightly higher than m
neutral
,
and first investigate the influence of different objective functions for fixed k. Figure 3.11 shows
comparisons of the different amplification factors for k = 0.05 and optimized solutions with cost
functions E
Y
, E
K
and E
T
, from left to right, respectively. Generally, we observe that G
T
(also G
K
)
is larger than G
Y
, as before, but interestingly G
T
(also G
K
) displays transient growth for all types of
cost function optimization (Butler & Farrell, 1992; Schmid & Henningson, 2001; Kaminski et al.,
2014). The peak of this transient growth happens at approximately the same time (for different
optimization functions) but the decay is substantially slower when using E
Y
, followed by E
T
and
much faster when using E
K
. Evidently, the intermediate decay displayed by the E
T
optimization
is a combination of the effects of purely optimizing on the scalar variance E
Y
(slower decay) and
39
-5 0 5
-100
-50
0
50
100
0.05
0
-0.05
η
z
(a)
-5 0 5
-100
-50
0
50
100
0.1
0
-0.1
η
z
(b)
-5 0 5
-100
-50
0
50
100
0.8
0
-0.8
η
z
(c)
-5 0 5
-100
-50
0
50
100
1
0
-1
η
z
×10
4
(d)
-5 0 5
-100
-50
0
50
100
1
0
-1
η
z
×10
4
(e)
-5 0 5
-100
-50
0
50
100
0.5
0
-0.5
η
z
×10
3
(f)
Figure 3.10: Two-dimensional optimal perturbation fields u
1
(a,d), w
1
(b,e) and Y
1
(c,f) at initial
time τ
0
(a–c) and at time (τ = 31.64) when G
T
≈ 10 (d–f ) of typical growing disturbance for R= 5,
M = 1, m= 0.05 and τ
0
= 10.
40
5 10 15 20
-6
-4
-2
0
2
4
τ− τ
0
G
5 10 15 20
-6
-4
-2
0
2
4
τ− τ
0 G
5 10 15 20
-6
-4
-2
0
2
4
G
T
G
K
G
Y
τ− τ
0 G
Figure 3.11: Optimal amplification factors G
T
, G
K
and G
Y
with different cost functions E
Y
(left),
E
K
(center) and E
T
(right) for R= 5, M = 1, τ
0
= 10, m= 0.16 and k= 0.05.
the kinetic energy of the flow E
K
(faster decay). The curves of G
Y
for all these cases show also
transient behaviour but of different nature. When optimizing for E
Y
and E
T
(left and right figures),
we observe an initial decay of the scalar variance, which was attributed earlier to diffusion. After
this decay, which is not observed when optimizing E
K
, there is transient growth towards a peak,
followed by decay. In the optimized solutions using E
Y
and E
T
, the scalar amplification factor
never exceeds the initial scalar variance.
Now, we concentrate on optimized solutions with the total energy cost function, E
T
, since it
encompasses both flow and scalar behaviour. Figure 3.12 displays the evolution of G
T
for different
values of the spanwise wavenumber k. First, it is apparent that there is short-term transient growth
at all k plotted, with G
T,max
larger than that of the purely two-dimensional mode k= 0. This can be
attributed to the lift-up effect of fluid elements in shear flows (Schmid & Henningson, 2001). The
behaviour of G
T,max
is not monotonic with k, however. Factor G
T,max
increases with increasing
k up to k≈ 0.1, followed by a decrease of the amplification factor for larger k. All perturbations
eventually decay, with the two-dimensional mode experiencing the slowest rate of decay at large
times, which indicates that three-dimensional modes are slightly more important early on during
the evolution. Note that the differences in amplification factors are not large for all values of k
shown.
Figure 3.13 shows the optimal perturbation fields for the two-dimensional case (corresponding
to the black curve in figure 3.12) of R= 5, M= 1, m= 0.16 and τ
0
= 10, at initial time τ
0
and at the
41
10 20 30 40 50
-2
-1
0
1
2
3
k = 0 (2D)
k = 0.02
k = 0.05
k = 0.1
k = 0.16
τ− τ
0
G
T
Figure 3.12: E
T
-optimal amplification factors for R = 5, M = 1, τ
0
= 10, and m = 0.16, and
different values of k.
time (τ = 15.52) when maximum amplification is reached: G
T,max
≈ 2.53. As seen before at lower
m, initially the flow patterns are tilted slightly against the base-flow shear but, as time progresses,
the disturbances eventually rotate to become more aligned with the mean shear direction and cause
transient growth of energy, also known as the Orr mechanism (Schmid & Henningson, 2001).
Overall, the behaviour is similar to that displayed by the perturbations at lower m in figure 3.10.
Finally, figure 3.14 shows the optimal perturbation fields for a three-dimensional case (corre-
sponding to the green curve in figure 3.12) of R = 5, M = 1, k = 0.05, m = 0.16, and τ
0
= 10, at
the initial time (left) and at the time (τ = 15.88) when maximum amplification is reached (right):
G
Y,max
≈ 2.69. The perturbation patterns are similar to those of figure 3.13 except that the oblique
structure of the mode is now evident. The presence of a spanwise component of perturbation v
1
is
also worth mentioning, and noticing that it has magnitude similar to that of u
1
.
42
-5 0 5
-50
-25
0
25
50
0.05
0
-0.05
η
z
(a)
-5 0 5
-50
-25
0
25
50
0.15
0
-0.15
η
z
(b)
-5 0 5
-50
-25
0
25
50
1
0
-1
η
z
(c)
-5 0 5
-50
-25
0
25
50
7
0
-7
η
z
(d)
-5 0 5
-50
-25
0
25
50
4
0
-4
η
z
(e)
-5 0 5
-50
-25
0
25
50
1
0
-1
η
z
(f)
Figure 3.13: Two-dimensional (k = 0) optimal perturbation fields u
1
(a,d), w
1
(b,e) and Y
1
(c,f) at
initial time τ
0
(a–c) and at time (τ = 15.52) when G
T,max
≈ 2.53 (d–f ) for R= 5, M = 1, m= 0.16
and τ
0
= 10.
43
η
y
z
u
1
(a)
u
1
(b)
v
1
(c)
v
1
(d)
w
1
(e)
w
1
(f)
Y
1
(g)
Y
1
(h)
Figure 3.14: Three-dimensional optimal perturbation fields u
1
, v
1
, w
1
and Y
1
at τ
0
= 10 (left) and
τ = 15.88 (right) for R= 5, M = 1, k= 0.05 and m= 0.16.
44
10
-2
10
-1
10
0
10
-1
10
0
10
1
DNS
m
peak
m
peak
/ 2
m
peak
/ 3
m
peak
/ 4
m
peak
/ 5
m
peak
/ 6
(t+t
i
)/T
DNS
δ
δ
i
q
t
i
T
DNS
10
-2
10
-1
10
0
10
-3
10
-2
10
-1
10
0
(t+t
i
)/T
DNS
δ
1
δ
i
q
t
i
T
DNS
Figure 3.15: Temporal evolution of the shear-layer width (δ/δ
i
)
p
t
i
/T
DNS
including base profile
(left), and perturbation width excluding the base profile (δ
1
/δ
i
)
p
t
i
/T
DNS
(right) for forward inte-
grations with R= 5, M = 1, τ
0
= 1.12 and several values of m, compared against DNS result taken
from Gat et al. (2017).
3.4 Comparison with DNS initialized with random perturbation
fields
The DNS study carried out by Gat et al. (2017) considered uniform viscosity, M = 1, for several
density ratios R. They analysed the behaviour of the thickness of the interface with time and
other nonlinear statistics not covered in this study. Figure 3.15 displays the evolution of the shear-
layer width, which indicates that there are two mixing-layer temporal evolution regimes: an initial
diffusion-dominant regime with growth rate∼
√
τ, followed by a turbulence-dominated regime
with growth rate ∼ τ
3
. The early-time growth corresponds to viscous spreading contained in
the base (self-similar) profiles, whereas the growing disturbances which dominate the base flows
promote much faster growth at later times, as described below.
In order to compare the DNS results with LSA, first one needs to translate the parameters of
their initial conditions to the present variables. Defining the initial shear-layer thickness δ
i
as the
width at which 0.1 < Y < 0.9 at initial time t
i
, from figure 3(b) of Gat et al. (2017) we obtain
δ
i
≈ 0.02ℓ = 0.13 for R = 5, where ℓ = 2π. Meanwhile, this corresponds to Δη ≈ 2.96 in our
45
self-similar Y
0
profile in figure 3.2. By using Eq. (3.7b) and Eq. (3.13), we can determine their
initial (start-up) time, given in our non-dimensional units by
τ
0
=

δ
i
/L
Δη

2
1
2K
=
"
δ
i
/( ¯ ν
2
g
−1
)
1/3
Δη
#
2
RSc
2(R+ 1)
=
"
0.13/(0.0044
2/3
)
2.96
#
2
5
12
= 1.12, (3.40)
for their given μ = 0.0044 and g= 1. Secondly, we need to match their perturbation wavenumber
to our m. From figure 15(a) in Gat et al. (2017), for “Pert3” on 1024
3
grid, the dimensionless
wavenumber corresponding to the peak of the perturbation spectrum is k
r
Δx≈ 0.15, which yields
k
r
≈ (0.15)(1024)/(4π)= 12.22, since their domain size is 4π. We non-dimensionalize this re-
sult to obtain our m
peak
= k
r
L = (12.22)( ¯ ν
2
g
−1
)
1/3
= 0.33. Subsequently, we carry out forward
integration for R= 5, M = 1, τ
0
= 1.12 for several wavenumbers, starting from m
peak
down to its
sixth subharmonics. Figure 3.15 shows the temporal evolutions of the shear-layer thickness δ (left
figure) and the width of the perturbation δ
1
(right figure) compared against the DNS from Gat et al.
(2017). The horizontal axis is scaled with the time scale used in their simulation, T
DNS
, given by
T
DNS
= 2π
s
ℓ
A g
= 19.29, where A =
R− 1
R+ 1
=
2
3
for R= 5. (3.41)
We observe that the transition from viscous decay to rapid growth happens at (t+t
i
)/T
DNS
≈ 0.1.
This is comparable with the transition times observed in the DNS given by (t+ t
i
)/T
DNS
≈ 0.25.
The slightly earlier transition time predicted by LSA owes its origin to the initial condition used
in the linearized study, which is obtained from the eigenmodes of the QSSA analysis, since these
modes are more unstable than the random two-dimensional field, with a wide spectrum, used in
the DNS. Despite this difference, the agreement is quite good and helps to connect the linear and
nonlinear studies. It appears that the peak wave mode chosen in the DNS was not an amplifying
one, and instead the instability was triggered by lower wave modes present in the initial condition.
46
Chapter 4
Nonlinear simulation
The zero-Mach-number Navier-Stokes solver developed by Lu & Pantano (2020), hereinafter is
referred as zMachDNS, is employed for our purpose. The following sections follow closely the
exposition in the paper, in which they also extensively assess the convergence, stability, and other
numerical properties. Uniform viscosity (M = 1) is assumed in all simulations considered in this
chapter.
4.1 Problem formulation
In zMachDNS, we define the following non-dimensionalization:
˜ x x x=
x x x
˜
L
, ˜ τ =
t
˜
T
, ˜ u u u=
u u u
˜
U
, ˜ p=
p
ρ
1
˜
U
2
, ˜ ρ =
ρ
ρ
1
, (4.1a)
where
˜
L,
˜
T , and
˜
U are the length, time and velocity scales respectively. Here we use tilde notation
to distinguish the non-dimensional variables of DNS from those of LSA. Furthermore, we have the
Reynolds number and Froude number defined as
Re=
ρ
1
˜
U
˜
L
μ
, Fr
2
=
g
˜
L
˜
U
2
, (4.1b)
47
from which we obtain
˜
L=

μReFr
ρ
1
√
g

2/3
,
˜
U =

μ gRe
ρ
1
Fr
2

1/3
,
˜
T =
˜
L
˜
U
=

μReFr
4
ρ
1
g
2

1/3
. (4.1c)
By choosing Fr= 1 and Re= ρ
1
/ ¯ ρ = 2R/(R+1), we will have the same length, time and velocity
scales between LSA and DNS. Now we can relate the two sets of non-dimensional variables in
Eq. (4.1a) and Eq. (3.7a), as
˜ x x x= x x x
∗
, ˜ τ = τ
∗
, ˜ u u u= u u u
∗
, ˜ p= p

¯ ρU
2
ρ
1
˜
U
2

=
p
∗
Re
, ˜ ρ =
ρ
∗
Re
. (4.2)
The non-dimensional equations for conservation of linear momentum and species concentration
can be written in compact form (dropping the tilde notations) as
ρ
∂
∂τ



u u u
Y



=



−∇p− ρu u u· ∇u u u+
1
Re
∇
2
u u u+ ρ

1
ρRe
− 1

ˆ z z z
1
Pe
∇
2
Y− ρu u u· ∇Y



≡



−∇p−N(u u u)+ ρ f f f
−D(Y)



,
(4.3)
where Pe= ReSc is the P´ eclet number. The operatorsN(u u u) andD(Y) are introduced for conve-
nience. Furthermore, the dimensionless density is governed by
ρ =
1
1− α(1−Y)
, (4.4)
where α≡ 1− R. By taking time-derivative of Eq. (4.4) and substituting the partial derivatives of
ρ and Y with respect to time from continuity and scalar equations, we can also express the mass
constraint equation as
C(ρu u u,Y)=
1
ρ
∇·(ρu u u)+ αD(Y)= 0. (4.5)
The mass constraint equation must be satisfied at initial time τ = τ
0
, otherwise spurious temporal
transients will appear which destroy solution convergence under time and space refinement (Lu,
2020). This is taken into consideration when we discuss the numerical scheme in the next section.
48
Now, we can obtain the pressure Poisson-like equation (PPE), by taking time-derivative of the
mass constraint equation, Eq. (4.5), and substituting the partial derivatives of u u u and Y with respect
to time from Eq. (4.3), giving
−∇·

1
ρ
∇p

= ∇·

1
ρ
N(u u u)− f f f

−
α
Pe
∇
2

1
ρ
D(Y)

. (4.6)
Nevertheless, solving Eq. (4.6) directly is costly due to variable-coefficient Laplacian operator,
and therefore an appropriate temporal numerical approximation which enforces mass constraint is
utilized instead. Our goal is then to solve Eqs. (4.3)-(4.5) numerically, which is discussed next.
4.2 Discretization
The governing equations are discretized in space using finite-differences and a pseudo-spectral
approximation in non-periodic and periodic directions, respectively. The boundary conditions for
the finite-differences are enforced with Summation-By-Parts/Simultaneous-Approximation-Terms
(SBP/SAT) method, a penalty-based approach; see reviews by Fern´ andez et al. (2014) and Sv¨ ard
& Nordstr¨ om (2014). The SBP/SAT method imposes boundary conditions in a weak sense by
introducing penalty terms into the governing equations and using finite-differencing operators that
preserve SBP properties to ensure numerical stability of the boundary treatment. The procedure
was introduced by Nordstr¨ om et al. (2007) for incompressible flows and extended for variable-
density flows by Lu & Pantano (2020). The semi-implicit second-order Adams-Bashforth/Crank-
Nicolson (AB/CN) method is employed for time-marching, in which advection (diffusion) term is
treated explicitly (implicitly) using an AB (CN) technique.
In two-dimensional case, we have the following density-weighted discrete operators
K
ρ
()≡
1
ρPe
∇
2
h
, V V V
ρ
()≡
1
ρ Re



∇
2
h
0
0 ∇
2
h


, G G G
ρ
≡
1
ρ
∇
h
, (4.7)
49
where stands as a placeholder for a given field, ∇
h
= D
x
⊗ I
y
+ I
x
⊗ D
y
is the discrete gradient
operator,⊗ and subscript h denote Kronecker product and grid spacing, respectively. The matrix
D denotes first-order differentiation matrix, where subcripts x and y denote the corresponding
coordinate directions. The density and viscosity in Eq. (4.7) are evaluated at the grid location
where each operator on the right-hand side returns its value, i.e. the finite-difference convention.
The advective part of the equations is discretized with the following nonlinear operator in skew-
symmetric form
N N N
ρ
(u u u,)≡
1
2ρ
[∇
h
·(ρu u u)+ ρu u u· ∇
h
− ∇
h
·(ρu u u)], (4.8)
which is more numerically stable than convection form (u u u· ∇
h
) for finite-difference methods
(Zang, 1991). We note that the SBP boundary stencils have been included in the finite-differencing
matrices. Hence, Eq. (4.3) can be written in semi-discrete form as
d
dτ



u u u
Y



=



−G G G
ρ
p− N N N
ρ
(u u u,u u u)+V V V
ρ
(u u u)+ f f f + S
ρ
(u u u)
K
ρ
(Y)− N
ρ
(u u u,Y)+ S
ρ
(Y)



, (4.9)
whereas the discrete mass constraint is given by
C(ρu u u,Y)=
1
ρ
∇
h
·(ρu u u)− αρ

K
ρ
(Y)− N
ρ
(u u u,Y)+ S
ρ
(Y)

= 0. (4.10)
The boundary terms are incorporated in the SAT penalty terms given by
S
ρ
(u u u)= S
a
ρ
(u u u)+ S
d
ρ
(u u u)+ S
g
ρ
(u u u), S
ρ
(Y)= S
a
ρ
(Y)+ S
d
ρ
(Y), (4.11)
where superscripts a, d, and g denote advection, diffusion, and pressure gradient, respectively.
The emergence of the SAT terms come from spatial discretization which results in finite-difference
operators that include boundary points. The definitions of these SAT terms, their derivation for
Navier-Stokes equations, and other related aspects are elaborated in Appendices A and B in Lu &
Pantano (2020). Recall that in periodic directions without boundary stencils, the SAT terms vanish.
50
Eq. (4.9) is discretized in time, assuming with constant time step Δτ for simplicity, which reads
u u u
n+1
− u u u
n
Δτ
=
3
2
h
f f f
n
− N N N
ρ
n(u u u
n
,u u u
n
)+ S
a
ρ
n(u u u
n
)
i
−
1
2
h
f f f
n−1
− N N N
ρ
n−1

u u u
n−1
,u u u
n−1

+ S
a
ρ
n−1

u u u
n−1

i
+
1
2
h
V V V
ρ
n+1
,μ
n+1

u u u
n+1

+ S
d
ρ
n+1

u u u
n+1

i
+
1
2
h
V V V
ρ
n
,μ
n(u u u
n
)+ S
d
ρ
n(u u u
n
)
i
−G
ρ
n+1/2
p
n+1/2
+ S
g
ρ
n+1/2

u u u
n+1/2

, (4.12)
Y
n+1
−Y
n
Δτ
=
3
2
h
−N N N
ρ
n(u u u
n
,Y
n
)+ S
a
ρ
n(Y
n
)
i
−
1
2
h
−N N N
ρ
n−1

u u u
n−1
,Y
n−1

+ S
a
Y
n−1

Y
n−1

i
+
1
2
h
K
ρ
n+1

Y
n+1

+ S
d
ρ
n+1

Y
n+1

i
+
1
2
h
K
ρ
n(Y
n
)+ S
d
ρ
n(Y
n
)
i
, (4.13)
where n denotes the time index and Δτ = τ
n+1
−τ
n
. The density in the pressure term of Eq. (4.12)
is defined by
1
ρ
n+1/2
=
1
2

1
ρ
n
+
1
ρ
n+1

. (4.14)
Eqs. (4.12)-(4.13) are advanced in time subject to a CFL-like time-step restriction where
CFL= Δτ
∑
|u
k
|
max
Δx
k
< constant∼ 1, (4.15)
independent of viscosity, due to implicit treatment of diffusive terms. For simplicity, the SAT
penalty terms for non-periodic directions are omitted in the subsequent discussions. The time
integration here advances the scalar field Y and then update the density by Eq. (4.4). The veloc-
ity is subsequently updated, which exhibits one of the desirable features of the AB/CN splitting.
Defining δY ≡ Y
n+1
−Y
n
, we can first write Eq. (4.13) in delta-form as

I−
Δτ
2
K
ρ
n+1

δY =
Δτ
2
h
K
ρ
n + K
ρ
n+1

(Y
n
)− 3N N N
ρ
n(u u u
n
,Y
n
)+ N N N
ρ
n−1

u u u
n−1
,Y
n−1

i
. (4.16)
The equation is nonlinear due to the K
ρ
n+1 term and requires iteration. We use a Picard iteration
to solve for δY and update the density and the next time level, ρ
n+1
, using Y
n+1
= Y
n
+ δY and
51
Eq. (4.4). Next, defining δu u u≡ u u u
n+1
− u u u
n
and ˇ p≡ Δτ p
n+1/2
, we can express Eq. (4.12) in delta-
form as

I I I−
Δτ
2
V V V
ρ
n+1
,μ
n+1

δu u u+ G G G
ρ
n+1/2
ˇ p =
Δτ
2
h
V V V
ρ
n
,μ
n+V V V
ρ
n+1
,μ
n+1

(u u u
n
)− 3

f f f
n
− N N N
ρ
n(u u u
n
,u u u
n
)

−

f f f
n−1
− N N N
ρ
n−1

u u u
n−1
,u u u
n−1

i
. (4.17)
By defining the following operators for convenience
A A A
ρ
n+1
,μ
n+1≡ I I I−
Δτ
2
V V V
ρ
n+1
,μ
n+1, (4.18)
a a a≡
Δτ
2
h
V V V
ρ
n
,μ
n+V V V
ρ
n+1
,μ
n+1

(u u u
n
)− 3

f f f
n
− N N N
ρ
n(u u u
n
,u u u
n
)

−

f f f
n−1
− N N N
ρ
n−1

u u u
n−1
,u u u
n−1

i
,
(4.19)
Eq. (4.17) can be simply expressed as
A A A
ρ
n+1
,μ
n+1δu u u+ G G G
ρ
n+1/2
ˇ p= a a a. (4.20)
The constraint of mass conservation is then enforced at the end of the step by
Q
ρ
n+1δu u u=−C

ρ
n+1
u u u
n
,Y
n+1

, (4.21)
where we define another operator
Q
ρ
n+1≡

ρ
n+1

−1
∇
h
·

ρ
n+1

− α

ρ
n+1

N
ρ
n+1

,Y
n+1

. (4.22)
The velocity-pressure coupled system, Eqs. (4.20)-(4.21), needs to be solved without assembly
or inversion of large matrices. In order to solve large-scale variable-density flows efficiently, the
pressure projection method, i.e. non-incremental projection method, is implemented. First, we
52
solve for intermediate velocity variable, δu u u
∗
, by ignoring pressure gradient term (similar to solving
advection-diffusion equation), like in a fractional step method, according to
A A A
ρ
n+1
,μ
n+1δu u u
∗
= a a a (4.23)
Second, the momentum equations are rearranged, subtracting Eq. (4.23) from Eq. (4.20) and using
u u u
∗
= u u u
n
+ δu u u
∗
, to give



A A A
ρ
n+1
,μ
n+1 G G G
ρ
n+1/2
Q
ρ
n+1 0






δz z z
ˇ p



=



0
−C

ρ
n+1
u u u
∗
,Y
n+1




, (4.24)
where δz z z≡ u u u
n+1
− u u u
∗
= δu u u− δu u u
∗
. Eq. (4.24) explicitly exposes the saddle-point structure of
the momentum-pressure coupling, highlighting here that Q
ρ
n+16=
h
G G G
ρ
n+1/2
i
T
for variable-density
flows. A good strategy to solve this problem is to form an approximate factorization (Badia &
Codina, 2008) to express the system in a form that is amenable to block iteration. However, the
subtle issue of boundary condition must be treated carefully.
Multiplying the first row of Eq. (4.24) by Q
ρ
n+1R, whereR≡

ρ
n+1

−1
h
ρ
n+1/2
i
, and sub-
tracting from the second row and using the definition of A A A
ρ
n+1
,μ
n+1 in Eq. (4.18), we can rewrite
Eq. (4.24) as



A A A
ρ
n+1
,μ
n+1 G G G
ρ
n+1/2
Q
ρ
n+1

I I I−RA A A
ρ
n+1
,μ
n+1

−Q
ρ
n+1RG G G
ρ
n+1/2






δz z z
ˇ p



=



0
−C

ρ
n+1
u u u
∗
,Y
n+1




. (4.25)
Eq. (4.25) is in general ill-posed when using a collocated grid, which is apparent from expanding
the lower diagonal term, giving
−Q
ρ
n+1RG G G
ρ
n+1/2
=

ρ
n+1

−1
∇
h
· ∇
h
+ αρ
n+1
N
ρ
n+1

ρ
n+1

−1
∇
h
,Y
n+1

. (4.26)
53
Observe that Eq. (4.26) is still not in a satisfactory form because the operator ∇
h
· ∇
h
is too wide
in a collocated mesh arrangement and it has a null space of dimension larger than one (beyond the
hydrostatic kernel mode). This is a problem that is still present in variable-density flow and can
result in an ill-posed system since the additional term N
ρ
n+1 in Eq. (4.26) does not really remove
the difficulty. As for the collocated method, the solution of the system Eq. (4.25) is obtainable only
if one applies a residual minimization method and accepts the least-square solution. It is at this
point that the approximate projection method is introduced. The wide operator ∇
h
· ∇
h
is replaced
with the narrow ∇
2
h
that does not suffer from an enlarged null space and we admit a discretization
error of the size of the spatial approximation order. Therefore, we introduce the modified operator
−

Q
ρ
n+1RG G G
ρ
n+1/2

∗
=

ρ
n+1

−1
∇
2
h
+ αρ
n+1
N
ρ
n+1

ρ
n+1

−1
∇
h
,Y
n+1

, (4.27)
such that
−Q
ρ
n+1RG G G
ρ
n+1/2
=−

Q
ρ
n+1RG G G
ρ
n+1/2

∗
+ T
ρ
n+1, (4.28)
where
T
ρ
n+1 =

ρ
n+1

−1

∇
2
h
− ∇
h
· ∇
h

. (4.29)
This procedure is commonly taken to reduce the dimension of the null space. Note that the com-
mutator in Eq. (4.29) is of the order of the spatial discretization accuracy, and introducing the ap-
proximate projection does not alter the convergence of the overall method. The modified operator
in Eq. (4.27) is discretely consistent that uses skew-symmetric discretization, which then requires
an efficient iteration scheme, discussed in Appendix C of Lu & Pantano (2020). Eq. (4.25) can
now be written as the linear algebra problemA x x x= b b b, using the modified operator, as
A =



A A A
ρ
n+1
,μ
n+1 G G G
ρ
n+1/2
Q
ρ
n+1

I I I−RA A A
ρ
n+1
,μ
n+1

−Q
ρ
n+1RG G G
ρ
n+1/2

∗


, x x x=



δz z z
ˇ p


,
b b b =



0
−C

ρ
n+1
u u u
∗
,Y
n+1




. (4.30)
54
It is well known that Eq. (4.30) cannot be uniquely solved for arbitrary b b b in those situations
involving periodic domains or with Neumann pressure boundary conditions; this is not a problem
if one disposes of Dirichlet pressure boundary conditions. This happens because the pressure
admits a single hydrostatic (constant) mode that belongs to the null space ofA ; all other spurious
pressure modes have been precluded by Eq. (4.27). Before setting the magnitude of the hydrostatic
pressure mode, one needs to ensure solvability of the system. The well-known procedure to ensure
solvability of Eq. (4.30) is described by the Fredholm alternative theorem and it involves two
parts: first, the enforcement of a solvability condition on b b b and second, the specification of the
magnitude of the hydrostatic mode (typically set to zero). These two steps ensure a unique solution
of Eq. (4.30). Despite the lack of a physical interpretation to the modification on b b b, here, linear
algebra theory is used to provide a general solution to the problem. Concretely, one needs to find
the null eigenvector w w w (since there is only one in our case) associated with the adjoint ofA , say
w w w
T
A = 0, or A
T
w w w= 0, (4.31)
and then enforce that b b b is orthogonal to w w w, i.e. b b b
T
w w w= 0. SinceA is not self-adjoint in variable-
density flows, we note that the vector w w w is not constant. We first remove the part of b b b that projects
into w w w, accomplished by
b b b→ b b b
w
= P(w w w)b b b, (4.32)
where P(w w w) =

I I I− w w ww w w
T

(we assumed||w w w||
2
= 1 where||||
2
denotes the Euclidean norm). We
then solveA x x x= b b b
w
and set the hydrostatic pressure mode to a constant (preferable zero) to obtain a
unique solution; this is the second part of the solvability requirement. We now proceed to determine
the vector w w w, by noting that the null vector ofA is the constant vector, i.e., the hydrostatic mode.
55
Let the null vector ofA
T
be denoted by w w w=[w w w
u
,w
p
]
T
, where we split into velocity and pressure
parts w w w
u
and w
p
, respectively. Then, Eq. (4.31) reads



A A A
T
ρ
n+1
,μ
n+1

I I I−RA A A
ρ
n+1
,μ
n+1

T
Q
T
ρ
n+1
G G G
T
ρ
n+1/2

−Q
ρ
n+1RG G G
ρ
n+1/2

∗T






w w w
u
w
p



= 0 0 0. (4.33)
Eq. (4.33) can be reduced to a single scalar equation in terms of w
p
, by observing that A A A
ρ
n+1
,μ
n+1
is not singular, giving
w w w
u
=−A A A
−T
ρ
n+1
,μ
n+1

I I I−RA A A
ρ
n+1
,μ
n+1

T
Q
T
ρ
n+1
w
p
. (4.34)
Then, introduce w w w
u
into the second of Eq. (4.33) and using Eq. (4.28) we get, after some rearrange-
ments,
w
T
p

Q
ρ
n+1A A A
−1
ρ
n+1
,μ
n+1
G G G
ρ
n+1/2
+ T
ρ
n+1

= 0. (4.35)
The general solution of Eq. (4.35) depends on the order of discretization, the boundary conditions
for the pressure, the grid and the time step.
In the simulation, one should find the null vector of the transpose ofA from Eq. (4.33) (at
each time step). Note that one only needs w
p
, which is needed to enforce b b b
T
w w w= 0 or equivalently
to determine b b b
w
, since w w w
u
is always zero. However it is undesirable to use Eq. (4.35) because
it requires the inverse of A A A
ρ
n+1
,μ
n+1. The preferred approach is to solve a linear system using
Eq. (4.31) directly by setting one element of w
p
to a constant and moving that element to the
right-hand side; after the solution is found, the vector can be normalized. This is possible because
we know that the dimension of the null space is exactly one. It is very important to enforce
solvability to ensure convergence of the iterative method used to solve Eq. (4.30). We observe that
the convergence rate of the iterative linear solver stalls if the solvability condition is not enforced.
In this state, the solution returned by the linear solver may be the least-square solution, which is the
same as that obtained when enforcing the solvability condition, but at a much higher computational
cost. However, if the iterative linear solver quits too soon, and the least-square solution is not
56
reached, the simulation can become quickly unstable. Most problems observed in variable-density
solvers at sufficiently high density ratios are a consequence of the attempt to solve the singular
Eq. (4.30) without paying attention to the solvability condition. Finally, note that these problems
do not arise in the pseudo-spectral method, since the spatial commutation properties can be applied
and w
p
is constant, which is consistent with the observation that spectral codes for variable-density
flows work relatively well.
4.3 Simulation setup
We carry out two-dimensional simulations where the Schmidt number is unity in all cases. First, we
initialize the base states u
0
(x,τ), w
0
(x,τ) and Y
0
(x,τ), which depend on R (see §3.2.1), at τ = τ
0
.
Then, linearly optimal perturbations (LOPs) obtained from the LSA (as described in chapter 3) are
given a finite initial amplitude ε and added to the base states. We note here that the initialization
is single-mode, where each wavenumber case uses its respective LOP, in order to isolate the in-
fluence of different wavemodes and to focus on the nonlinear cascade process. Referring to the
illustration in figure 3.1, the z-direction is periodic and the domain is set to be L
z
= 2π/m, in order
to span exactly one initial wavelength, whereas the x-direction is discretized with a second-order
finite-difference (FD) scheme and outflow boundary conditions are imposed. The grid spacing is
Δx = Δz≈ 0.08 in order to sufficiently resolve the mixing layer. Verification includes the effects
of using a finite computational domain, as well as temporal and grid convergence due to different
resolutions. A typical initialization for the simulation is shown in figure 4.1 for scalar and ve-
locity profiles in case of R = 5, where we can observe the base states are perturbed by the scalar
disturbance with single mode m= 0.05 and finite amplitude ε = 0.01 at τ = τ
0
= 1.
57
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(a)
-60 -30 0 30 60
-60
-30
0
30
60
0
-0.2
-0.4
-0.6
x
z
(b)
-60 -30 0 30 60
-60
-30
0
30
60
2
1.2
0.4
-0.4
x
z
(c)
Figure 4.1: Flow field initializations for two-dimensional DNS with scalar Y (a), velocities u (b)
and w (c) shown for R = 5. Single-mode scalar LOP with m = 0.05 and ε = 0.01 is imposed at
τ = τ
0
= 1.
4.4 Results
In order to compare the evolution of perturbation energies, first we compute E
ψ
1
(τ
1
) from DNS,
where ψ
1
=(u
1
,v
1
,w
1
,Y
1
), as follows
E
ψ
1
(τ
1
)=
2π
Z
0
∞
Z
−∞
[εψ
1
(x,z,τ
1
)]
2
dxdz, (4.36)
which can be related to that of the LSA, given by Eq. (3.27), by utilizing the well-known Parseval’s
identity
E
ψ
1
(τ
1
)=
2π
Z
0
∞
Z
−∞
|εψ
1
(x,z,τ
1
)|
2
dxdz= 2π
∞
∑
n=−∞
∞
Z
−∞
|ε ˆ ψ
n
(x,τ
1
)|
2
dx= 2π
∞
Z
−∞
|ε ˆ ψ
n
(x,τ
1
)|
2
dx,
(4.37)
since we only have one particular mode n for each case. Finally, by using Eq. (3.13), we obtain the
comparable quantity in LSA,
E
ψ
1
(τ
1
)= 2πε
2
p
2Kτ
1
∞
Z
−∞
| ˆ ψ
k
(η,τ
1
)|
2
dη = 2πε
2
p
2Kτ
1
E
ˆ ψ
(τ
1
). (4.38)
58
We now concentrate on cases initialized with only scalar perturbation Y
1
. Figure 4.2 shows the
evolutions of the scalar and velocity perturbation energies as functions of τ
0
for R = 5, m= 0.05
and ε = 0.01, where solid (dashed) lines denote the DNS (LSA) cases. Here we can see that
initially both curves follow similar trajectory in the linear phase. However once the perturbation
amplitudes grow sufficiently large, nonlinear effects become more significant as shown by devi-
ation and saturation of the DNS perturbation growth. In general, E
Y
1
saturates a few orders of
magnitude lower than the saturation levels of E
u
1
and E
w
1
, again indicating that the amount of
energy in the scalar is typically small. In addition, the onset of nonlinearity takes place earlier
with increasing τ
0
, as expected based on the relationship between τ
0
and early Reynolds scaling,
Eq. (3.38), since higher kinetic energy promotes faster growth of the disturbances. Interestingly,
once nonlinearity sets in, the DNS vertical perturbation velocity w
1
grows faster than that of LSA,
with larger deviation from LSA for increasing τ
0
. This suggests stronger energy transfer from the
growing momentum of the base flow w
0
.
In order to observe the energy cascade taking place from the initial single mode m into other
wavemodes in the nonlinear simulation, we first express the scalar perturbation Y
1
(x,z,τ) as a
complex Fourier series in z,
Y
1
(x,z,τ)=
∞
∑
n=−∞
c
n
(x,τ)e
inmz
, (4.39)
where the Fourier coefficients c
n
(x,τ) are given by
c
n
(x,τ)=
1
L
z
L
z
Z
0
Y
1
(x,z,τ)e
−inmz
dz. (4.40)
Subsequently, we compute the power spectrum P
n
(x,τ)= c
n
(x,τ)c
∗
n
(x,τ). We consider the case of
R = 5, m = 0.05, τ
0
= 1 and ε = 0.01, which corresponds to the blue curve in figure 4.2(a). The
power spectra for the fundamental and its harmonics at τ− τ
0
= 0,12,19,24,30 and 50 are shown
in figure 4.3(a–f ), which correspond to times at the beginning, linear phase, onset of nonlinearity,
peak perturbation energy, early and late saturation phases, respectively. Note the difference in
vertical scales in the plots shown here. Initially, we can only see the non-zero P
1
(fundamental
59
mode) (figure 4.3(a)), with its peak in the heavier fluid region. This single mode grows during
linear phase with no energy cascade to other modes as shown in figure 4.3(b), with the peak shifting
into lighter fluid region. At the onset of nonlinearity, figure 4.3(c), the energy cascade to the first
and second harmonics can be observed, while the nucleation of a second peak in the fundamental
mode begins to form. Note that only multiple integers of the fundamental mode, i.e. n = 1,2,...
have nonzero values, whereas there is no energy cascade to n = 3/2,5/2,... modes. Once the
perturbation energy reaches its peak, there are two peaks at the fundamental mode, while the
energy has cascaded further into the third harmonics (figure 4.3(d)). Perturbation then decays
during early nonlinear saturation phase, figure 4.3(e), with multiple peaks form at fundamental and
its harmonics. Finally, at later saturation stage the harmonics are damped, while the fundamental
mode spreads further into the lighter fluid region (figure 4.3(f )).
We are also interested on how initial finite amplitude of perturbations influences the solutions.
Figure 4.4 displays the evolutions of scalar and velocity perturbation energies as functions of ε for
R= 5, m= 0.05 and τ
0
= 10. The energies are normalized with ε
2
in order to collapse the curves.
The evolutions for DNS cases follows closely those of LSA for longer period of time at smaller
ε, since it approaches the infinitesimal limit ε≪ 1. On the other hand, larger ε triggers onset of
nonlinearity and saturation of perturbation growths earlier. As observed in the case of varying τ
0
above, w
1
grows faster than that of LSA with larger deviation for increasing ε as nonlinear effects
become more significant. Finally, we vary the density ratio R with the evolutions of perturbation
energies shown in figure 4.5 for m = 0.05, τ
0
= 10 and ε = 0.01. Earlier onset of nonlinearity
and saturation happen at higher R due to presence of larger potential energy in the flow. Here, for
relatively low R = 2.5, the DNS w
1
does not grow faster than that of LSA, whereas for R≥ 5,
again we observe w
1
grows more strongly in DNS than in LSA, with larger deviation from LSA
for increasing R.
The physical features of the flow fields are shown for R = 5, m = 0.05, τ
0
and ε = 0.01 for
τ
0
= 10 (figure 4.6) and for τ
0
= 50 (figure 4.7). We note here that the simulations were carried out
with larger domain in the x direction than shown in the figures, which ensures that the boundaries
60
0 10 20 30
10
-2
10
-1
10
0
10
1
10
2
10
3
10
4
10
5
10
6
τ− τ
0
E
Y
1
τ
0
= 50
τ
0
= 20
τ
0
= 10
τ
0
= 1
(a)
0 10 20 30
10
-2
10
0
10
2
10
4
10
6
10
8
τ− τ
0
E
u
1 τ
0
= 50
τ
0
= 20
τ
0
= 10
τ
0
= 1
(b)
0 10 20 30
10
-2
10
0
10
2
10
4
10
6
10
8
τ− τ
0
E
w
1
τ
0
= 50
τ
0
= 20
τ
0
= 10
τ
0
= 1
(c)
Figure 4.2: The evolutions of energy of perturbation scalar Y
1
(a), velocities u
1
(b) and w
1
(c)
obtained from DNS (solid) and LSA (dashed), as functions of τ
0
for R= 5, m= 0.05 and ε = 0.01.
61
0
100
2
P
4
4
x
0 P
3
6
P
2
-100
P
1
(a)
0
100
100
P
4
200
x
0 P
3
300
P
2
-100
P
1
(b)
0
100
P
4
5000
x
0 P
3
10000
P
2
-100
P
1
(c)
0
0.5
100
P
4
1
10
4
1.5
x
0 P
3
2
P
2
-100
P
1
(d)
0
100
P
4
5000
x
0 P
3
10000
P
2
-100
P
1
(e)
0
100
1000
P
4
2000
x
0 P
3
3000
P
2
-100
P
1
(f)
Figure 4.3: Power spectrum P
n
(x,τ) for R = 5, m = 0.05, τ
0
= 1 and ε = 0.01 at τ− τ
0
=
0,12,19,24,30 and 50, respectively (a–f ). See text for significance of these times.
62
0 10 20 30
10
3
10
4
10
5
10
6
10
7
10
8
10
9
τ− τ
0
E
Y
1
ε
2
LSA
ε = 0.005
ε = 0.01
ε = 0.05
(a)
0 10 20 30
10
3
10
5
10
7
10
9
10
11
10
13
τ− τ
0
E
u
1
ε
2
LSA
ε = 0.005
ε = 0.01
ε = 0.05
(b)
0 10 20 30
10
3
10
5
10
7
10
9
10
11
10
13
τ− τ
0
E
w
1
ε
2
LSA
ε = 0.005
ε = 0.01
ε = 0.05
(c)
Figure 4.4: The evolutions of energy of perturbation scalar Y
1
(a), velocities u
1
(b) and w
1
(c)
obtained from DNS (solid) and LSA (dashed), as functions of ε for R= 5, m= 0.05 and τ
0
= 10.
63
0 10 20 30 40
10
-1
10
0
10
1
10
2
10
3
10
4
τ− τ
0
E
Y
1
R= 10
R= 7.5
R= 5
R= 2.5
(a)
0 10 20 30 40
10
-2
10
0
10
2
10
4
10
6
10
8
τ− τ
0
E
u
1
R= 10
R= 7.5
R= 5
R= 2.5
(b)
0 10 20 30 40
10
-2
10
0
10
2
10
4
10
6
10
8
τ− τ
0
E
w
1
R= 10
R= 7.5
R= 5 R= 5
R= 2.5
(c)
Figure 4.5: The evolutions of energy of perturbation scalar Y
1
(a), velocities u
1
(b) and w
1
(c)
obtained from DNS (solid) and LSA (dashed), as functions of R for m= 0.05, τ
0
= 10 and ε = 0.01.
64
do not interfere with the mixing layer. The left columns exhibit the flow fields when perturbations
are still growing within the linear phase of τ− τ
0
≈ 7, while the middle columns display the fields
when scalar perturbation energies reach their maxima, according to figure 4.2(a). Finally, the right
columns show the saturation phase where LSA solutions no longer apply. Comparison of the two
sets of figures indicate that a KHI-like billow is formed where the two fluids roll. However, the
billow seems to be more swept at higher τ
0
due to stronger shear acting on the mixing layer. This
apparently results in a lower mixing layer growth at larger τ
0
case in the nonlinear phase. In order
to quantify this, we measure the width of the mixing layer, δ(τ), which is defined as the horizontal
extent at which the mean scalar
¯
Y(x,τ) falls within the range
0.1<
¯
Y(x,τ)< 0.9, (4.41)
where
¯
Y(x,τ) is averaged over the periodic direction z, given by
¯
Y(x,τ)=
1
L
z
L
z
/2
Z
−L
z
/2
Y(x,z,τ)dz. (4.42)
Figure 4.8 shows the evolution of the mixing layer widths δ(τ) for R = 5, m = 0.05, ε = 0.01
and different τ
0
cases. As expected, initially the shear layer width at τ
0
= 50 is the largest, since
the base scalar Y
0
has spread the furthest due to diffusion. In addition, from the growth rate evo-
lutions we can observe the distinction between linear, diffusive region and the nonlinear, chaotic
phases. Interestingly, the behavior of mixing layer growth is again not monotonic with τ
0
. Once
it passes the diffusion phase, the τ
0
= 10 case experiences the earliest and steepest rise of mix-
ing layer width. On the other hand, the mixing layer growth rates of τ
0
= 50 case does not show
an appreciable difference between linear and nonlinear phases, such that it grows most slowly at
later times. This shows the strong influence of the shear (driven by kinetic energy) acting in the z
direction, which inhibits the expansion of the mixing region.
65
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(a)
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(b)
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(c)
-60 -30 0 30 60
-60
-30
0
30
60
0.4
-0.1
-0.6
x
z
(d)
-60 -30 0 30 60
-60
-30
0
30
60
25
12.5
0
-12.5
-25
x
z
(e)
-60 -30 0 30 60
-60
-30
0
30
60
25
15
5
-5
-15
x
z
(f)
-60 -30 0 30 60
-60
-30
0
30
60
35
20
5
-10
x
z
(g)
-60 -30 0 30 60
-60
-30
0
30
60
65
30
-5
-40
x
z
(h)
-60 -30 0 30 60
-60
-30
0
30
60
70
35
0
-35
x
z
(i)
Figure 4.6: Two-dimensional scalar Y (a–c), velocities u (d–f ) and w (g–i) fields at τ ≈ 17.17
(a,d,g), at τ = 28.71 (b,e,h) and at τ = 32.71 (c,f,i) for R= 5, m= 0.05, τ
0
= 10 and ε = 0.01. See
text for the significance of these three instances of time.
66
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(a)
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(b)
-60 -30 0 30 60
-60
-30
0
30
60
1
0.75
0.5
0.25
0
x
z
(c)
-60 -30 0 30 60
-60
-30
0
30
60
2.5
0
-2.5
x
z
(d)
-60 -30 0 30 60
-60
-30
0
30
60
15
5
-5
-15
x
z
(e)
-60 -30 0 30 60
-60
-30
0
30
60
15
5
-5
-15
x
z
(f)
-60 -30 0 30 60
-60
-30
0
30
60
125
75
25
-25
x
z
(g)
-60 -30 0 30 60
-60
-30
0
30
60
135
95
55
15
-25
x
z
(h)
-60 -30 0 30 60
-60
-30
0
30
60
150
105
60
15
-30
x
z
(i)
Figure 4.7: Two-dimensional scalar Y (a–c), velocities u (d–f ) and w (g–i) fields at τ ≈ 57.02
(a,d,g), at τ = 65.98 (b,e,h) and at τ = 72.69 (c,f,i) for R= 5, m= 0.05, τ
0
= 50 and ε = 0.01. See
text for the significance of these three instances of time.
67
0 10 20 30
0
30
60
90
τ− τ
0
δ(τ)
τ
0
= 50
τ
0
= 20
τ
0
= 10
τ
0
= 1
Figure 4.8: Mixing layer widths δ(x,τ) for R= 5, m= 0.05, ε = 0.01 as functions of τ
0
.
Subsequently, we investigate the condition close to the neutral amplification wavenumber,
where we choose m= 0.16 for R= 5 (following the parameter chosen in the discussion in §3.3.2).
Figure 4.9 displays the evolution of the perturbation energies for R = 5, τ
0
= 10 and ε = 0.01.
As reported by Kaminski et al. (2017), the nonlinearity in DNS results in lower peak energies,
while the decays of perturbation energies are slower at large times as compared to those of LSA.
Nevertheless, they emphasize that it is the linear growth mechanism which results in disturbance
energy growth, rather than some nonlinear response to a large disturbance, even when nonlinear
saturation is significant here.
68
0 10 20
10
-4
10
-3
10
-2
10
-1
10
0
10
1
τ− τ
0
E
E
u
1
E
w
1
E
Y
1
Figure 4.9: The evolutions of energy of perturbation scalar Y
1
, velocities u
1
and w
1
obtained from
DNS (solid) and LSA (dashed) for R= 5, m= 0.16, τ
0
= 10 and ε = 0.01.
69
Chapter 5
Conclusions
The hydrodynamic stability of two fluid columns with different densities and viscosities in a
gravitational field has been investigated in this dissertation. In this flow configuration, a uniform
gravity field acts in the direction transversal to the density stratification, resulting in opposing
constant acceleration of the free streams. The study is conducted in three coherent steps. First, the
inviscid LSA of an immiscible interface between the two columns is carried out in the presence of
surface tension. The analysis is developed in the framework of an initial value problem, since the
problem does not admit harmonic normal modes, and solved analytically. The results demonstrate
that the flow is unconditionally unstable, owing to the growing momentum of the free streams with
time. Furthermore, the growth of the solution is quadratic in time exponential.
Second, LSA of the interface separating two miscible fluid columns, now allowing for different
viscosities, is investigated and compared with previous inviscid analysis. We perform numerical
integration of the linear time-dependent IVP as a function of wavenumber. Furthermore, we carry
out adjoint optimization in order to determine the initial profiles that result in maximum amplifi-
cation at a finite time horizon. Within the strongly unstable range of wavenumbers, we confirm
that the evolution approaches the inviscid limit as the initial Reynolds number of the flow in-
creases. In this regime, the two-dimensional mode is the most unstable (dominant). At higher wave
modes, around the region of neutral stability, it is found that transient growth is present and that
Chapter adapted from Prathama & Pantano (2017, 2021).
70
three-dimensional modes (with a spanwise non-zero wavenumber) are more unstable than the two-
dimensional mode. The difference in amplification factors between two- and three-dimensional
modes is not very large, but it is noticeable.
We also investigate the effect of different optimization cost functions and concluded that using
the total energy of the flow, kinetic and potential, gives an overall picture of the evolution. Using
a cost function based only on scalar variance (potential energy) highlights initial transients of the
mixture fraction that are controlled by diffusion. Optimization of the initial scalar profiles can
largely remove this diffusional decay from the dynamics but hides flow effects (velocity) that are
important overall.
Third, the DNS of the flow configuration is performed under the uniform viscosity assump-
tion with the zero-Mach-number variable density solver developed by Lu & Pantano (2020). The
previously obtained linearly optimal perturbations are given finite amplitude, imposed on the flow
at initial time, and let to evolve into fully nonlinear manner. The nonlinear study shows that the
energy of the perturbations follow more closely the LSA for smaller initial amplitude of pertur-
bations, lower density ratio and starting time, as expected. As perturbations grow sufficiently
large, nonlinear effects become more significant as exhibited by saturation of the DNS perturba-
tion growth, followed by long time decay of perturbations. Energy cascade from the initial (funda-
mental) mode into higher modes (harmonics) is observed, starting from the onset of nonlinearity.
Physical features of the flow fields show formation of billows where the two fluids roll in the non-
linear phase, as in the Kelvin–Helmholtz instability. Measurement of mixing layer width indicates
that the evolution of mixing layer growth and the transition from the diffusion-dominant regime
into a chaotic-dominated regime are sensitive to the initial time when perturbations are imposed
into the flow. We found that there is an intermediate initial time, where the transition happens
earliest and the growth rate of mixing layer in the chaotic regime is the fastest. At low initial time,
the kinetic energy is still too low to drive the instability, while the potential energy is high due to
density jump over a thin interface thickness. The converse case is true at large initial time, where
the base flow has spread further but the vertical momentum has grown so much that the mixing
71
region is more swept, inhibiting the expansion of the mixing layer in the horizontal direction. The
overall results show that there are two mechanisms at play in this flow configuration: one involv-
ing the ever growing momentum (related to kinetic energy) in the vertical direction, and the other
involving diffusional spread (related to potential energy) in the horizontal direction.
72
Bibliography
ABRAMOWITZ, M. & STEGUN, I. A. 1965 Handbook of Mathematical Functions. Dover.
ARRATIA, C., CAULFIELD, C. P. & CHOMAZ, J.-M. 2013 Transient perturbation growth in time-
dependent mixing layers. J. Fluid Mech. 717, 90–133.
AURENTZ, J. L. & TREFETHEN, L. N. 2017 Block operators and spectral discretizations. SIAM
Rev. 59, 423–446.
BADIA, S. & CODINA, R. 2008 Algebraic pressure segregation methods for the incompressible
Navier–Stokes equations. Arch. Comput. Methods Eng. 15, 343–369.
BATCHELOR, G. K. 1967 An Introduction to Fluid Dynamics. Cambridge University Press.
BEN, Y., DEMEKHIN, E. A. & CHANG, H.-C. 2002 A spectral theory for small-amplitude mis-
cible fingering. Phys. Fluids 14, 999–1010.
BRETHERTON, C. S. 1987 A theory for nonprecipitating moist convection between two parallel
plates. Part I: Thermodynamics and “linear” solutions. J. Atmos. Sci. 44, 1809–1827.
BUTLER, K. M. & FARRELL, B. F. 1992 Three-dimensional optimal perturbations in viscous
shear flow. Phys. Fluids 4, 1637–1650.
CHANDRASEKHAR, S. 1981 Hydrodynamic and Hydromagnetic Stability. Dover.
CHERUBINI, S., PALMA, P. DE & ROBINET, J.-CH. 2015 Nonlinear optimals in the asymptotic
suction boundary layer: Transition thresholds and symmetry breaking. Phys. Fluids 27, 034108.
CORBETT, P. & BOTTARO, A. 2000 Optimal perturbations for boundary layers subject to stream-
wise pressure gradient. Phys. Fluids 12, 120–130.
CORBETT, P. & BOTTARO, A. 2001 Optimal linear growth in swept boundary layers. J. Fluid
Mech. 435, 1–23.
CRIMINALE, W. O., JACKSON, T. L. & JOSLIN, R. D. 2003 Theory and Computation of Hydro-
dynamic Stability. Cambridge University Press.
DANIEL, D., TILTON, N. & RIAZ, A. 2013 Optimal perturbations of gravitationally unstable,
transient boundary layers in porous media. J. Fluid Mech. 727, 456–487.
73
DOUMENC, F., BOECK, T., GUERRIER, B. & ROSSI, M. 2010 Transient Rayleigh–B´ enard–
Marangoni convection due to evaporation: a linear non-normal stability analysis. J. Fluid Mech.
648, 521–539.
DRAZIN, P. G. & REID, W. H. 2004 Hydrodynamic Stability. Cambridge University Press.
DRISCOLL, T. A., HALE, N., & TREFETHEN, L. N. 2014 Chebfun Guide. Pafnuty Publications.
DUGUET, Y., MONOKROUSOS, A., BRANDT, L. & HENNINGSON, D. S. 2013 Minimal transi-
tion thresholds in plane Couette flow. Phys. Fluids 25, 084103.
ETRATI, A., ALBA, K. & FRIGAARD, I. A. 2018 Two-layer displacement flow of miscible fluids
with viscosity ratio: Experiments. Phys. Fluids 30, 052103.
FARRELL, B. F. 1988 Optimal excitation of perturbations in viscous shear flow. Phys. Fluids 31,
2093–2102.
FARRELL, B. F. & MOORE, A. M. 1992 An adjoint method for obtaining the most rapidly grow-
ing perturbation to oceanic flows. J. Phys. Oceanogr. 22, 338–349.
FERN ´ ANDEZ, D. C. D. R., HICKEN, J. E. & ZINGG, D. W. 2014 Review of summation-by-
parts operators with simultaneous approximation terms for the numerical solution of partial
differential equations. Comput. Fluids 95, 171–196.
FOURES, D. P. G., CAULFIELD, C. P. & SCHMID, P. J. 2013 Localization of flow structures
using -norm optimization. J. Fluid Mech. 729, 672–701.
FOURES, D. P. G., CAULFIELD, C. P. & SCHMID, P. J. 2014 Optimal mixing in two-dimensional
plane Poiseuille flow at finite P´ eclet number. J. Fluid Mech. 748, 241–277.
FUJIMURA, K. 1989 The equivalence between two perturbation methods in weakly nonlinear sta-
bility theory for parallel shear flows. Proc. R. Soc. Lond. A 424, 373–392.
FUJIMURA, K. 1991 Methods of centre manifold and multiple scales in the theory of weakly
nonlinear stability for fluid motions. Proc. R. Soc. Lond. A 434, 719–733.
GAT, I., MATHEOU, G., CHUNG, D. & DIMOTAKIS, P.E. 2016 Acceleration-driven variable-
density turbulent flow. In Proceedings of the VIIIth International Symposium on Stratified Flows
(ISSF). http://escholarship.org/uc/item/61d722q6, San Diego, CA.
GAT, I., MATHEOU, G., CHUNG, D. & DIMOTAKIS, P.E. 2017 Incompressible variable-density
turbulence in an external acceleration field. J. Fluid Mech. 827, 506–535.
GEORGESCU, A. 1985 Hydrodynamic stability theory. Springer Netherlands.
GU ´ EGAN, A., SCHMID, P. J. & HUERRE, P. 2006 Optimal energy growth and optimal control in
swept Hiemenz flow. J. Fluid Mech. 566, 11–45.
GUSTAVSSON, L. H. 1991 Energy growth of three-dimensional disturbances in plane Poiseuille
flow. J. Fluid Mech. 224, 241–260.
74
HOTA, T. K., PRAMANIK, S. & MISHRA, M. 2015 Nonmodal linear stability analysis of miscible
viscous fingering in porous media. Phys. Rev. E 92, 053007.
KAMINSKI, A. K., CAULFIELD, C. P. & TAYLOR, J. R. 2014 Transient growth in strongly strat-
ified shear layers. J. Fluid Mech. 758, R4.
KAMINSKI, A. K., CAULFIELD, C. P. & TAYLOR, J. R. 2017 Nonlinear evolution of linear
optimal perturbations of strongly stratified shear layers. J. Fluid Mech. 825, 213–244.
KERSWELL, R.R. 2018 Nonlinear nonmodal stability theory. Annu. Rev. Fluid Mech. 50, 319–345.
KERSWELL, R. R., PRINGLE, C. C. T. & WILLIS, A. P. 2014 An optimization approach for
analysing nonlinear stability with transition to turbulence in fluids as an exemplar. Rep. Prog.
Phys. 77 (8), 085901.
KIM, M. C. & CHOI, C. K. 2011 The stability of miscible displacement in porous media: Non-
monotonic viscosity profiles. Phys. Fluids 23, 084105.
KIRSHBAUM, D. J. & STRAUB, D. N. 2019 Linear theory of shallow convection in deep, verti-
cally sheared atmospheres. Q. J. R. Meteorol. Soc. 145, 3129–3147.
LOPEZ-ZAZUETA, A., FONTANE, J. & JOLY, L. 2016 Optimal perturbations in time-dependent
variable-density Kelvin–Helmholtz billows. J. Fluid Mech. 803, 466–501.
LU, X. 2020 Hydrodynamics stability of a premixed flame subjected to transverse shear. PhD
thesis, University of Illinois at Urbana-Champaign.
LU, X. & PANTANO, C. 2020 On mass conservation and solvability of the discretized variable-
density zero-Mach Navier–Stokes equations. J. Comput. Phys 404, 109132.
LUCHINI, P. & BOTTARO, A. 2014 Adjoint equations in stability analysis. Annu. Rev. Fluid Mech.
46, 493–517.
MAJDA, A. & SETHIAN, J. 1985 The derivation and numerical solution of the equations for zero
mach number combustion. Comb. Sci. Tech. 42, 185–205.
MELLADO, J. P. 2017 Cloud-top entrainment in stratocumulus clouds. Annu. Rev. Fluid Mech. 49,
145–169.
NORDSTR ¨ OM, J., MATTSSON, K. & SWANSON, C. 2007 Boundary conditions for a divergence
free velocity-pressure formulation of the Navier–Stokes equations. J. Comput. Physics 225, 874–
890.
OLSON, B. J., LARSSON, J., LELE, S. K. & COOK, A. W. 2011 Nonlinear effects in the com-
bined Rayleigh–Taylor/Kelvin–Helmholtz instability. Phys. Fluids 23, 114107.
PRALITS, J. O., BOTTARO, A. & CHERUBINI, S. 2015 Weakly nonlinear optimal perturbations.
J. Fluid Mech. 785, 135–151.
75
PRATHAMA, A. H. & PANTANO, C. 2017 Inviscid linear stability analysis of two vertical columns
of different densities in a gravitational acceleration field. J. Fluid Mech. 826, R4.
PRATHAMA, A. H. & PANTANO, C. 2021 Linear stability analysis of two fluid columns of differ-
ent densities and viscosities in a gravity field. J. Fluid Mech. 920, A26.
PRINGLE, C. C. T. & KERSWELL, R. R. 2010 Using nonlinear transient growth to construct the
minimal seed for shear flow turbulence. Phys. Rev. Lett. 105 (15).
RABIN, S. M. E., CAULFIELD, C. P. & KERSWELL, R. R. 2012 Triggering turbulence efficiently
in plane Couette flow. J. Fluid Mech. 712, 244–272.
RAPAKA, S., CHEN, S., PAWAR, R. J., STAUFFER, P. H. & ZHANG, D. 2008 Non-modal growth
of perturbations in density-driven convection in porous media. J. Fluid Mech. 609, 285–303.
RAYLEIGH, LORD 1883 Investigation of the character of the equilibrium of an incompressible
heavy fluid of variable density. Proc. Lond. Math. Soc. 14, 170–177.
REDDY, S. C. & HENNINGSON, D. S. 1993 Energy growth in viscous channel flows. J. Fluid
Mech. 252, 209–238.
REDDY, S. C., SCHMID, P. J. & HENNINGSON, D. S. 1993 Pseudospectra of the Orr–Sommerfeld
operator. SIAM J. Appl. Math. 53, 15–47.
RIAZ, A., PANKIEWITZ, C. & MEIBURG, E. 2004 Linear stability of radial displacements in
porous media: Influence of velocity-induced dispersion and concentration-dependent diffusion.
Phys. Fluids 16, 3592–3598.
RICHTMYER, R.D. 1960 Taylor instability in shock acceleration of compressible fluids. Commu-
nications in Pure and Applied Mathematics 13 (2), 297–319.
SANDOVAL, D. L. 1995 The dynamics of variable-density turbulence. PhD thesis, University of
Washington.
SCHMID, P. J. 2007 Nonmodal stability theory. Annu. Rev. Fluid Mech. 39, 129–162.
SCHMID, P. J. & HENNINGSON, D. S. 2001 Stability and Transition in Shear Flows. Springer.
SHARP, D. H. 1984 An overview of Rayleigh–Taylor instability. Physica D 12, 3–18.
SQUIRE, H.B. 1933 On the stability of three-dimensional disturbances of viscous flow between
parallel walls. Proc. R. Soc. Lond. A 142 (621-8), 129–155.
STUART, J. T. 1960 On the non-linear mechanics of wave disturbances in stable and unstable
parallel flows Part 1. The basic behaviour in plane Poiseuille flow. J. Fluid Mech. 9, 353–370.
SV ¨ ARD, M. & NORDSTR ¨ OM, J. 2014 Review of summation-by-parts schemes for
initial–boundary-value problems. J. Comput. Phys. 268, 17–38.
TAN, C. T. & HOMSY, G. M. 1986 Stability of miscible displacements in porous media: Recti-
linear flow. Phys. Fluids 29, 3549–3556.
76
TAYLOR, G. 1950 The instability of liquid surfaces when accelerated in a direction perpendicular
to their planes. I. Proc. R. Soc. Lond. A 201, 192–196.
TREFETHEN, L. N., TREFETHEN, A. E., REDDY, S. C. & DRISCOLL, T. A. 1993 Hydrodynamic
stability without eigenvalues. Science 261, 578–584.
TRUZZOLILLO, D. & CIPELLETTI, L. 2018 Hydrodynamic instabilities in miscible fluids. J. Phys.
Condens. Matter 30, 033001.
WATSON, J. 1960 On the non-linear mechanics of wave disturbances in stable and unstable parallel
flows Part 2. The development of a solution for plane Poiseuille flow and for plane Couette flow.
J. Fluid Mech. 9, 371–389.
YAGLOM, A. M. 2012 Hydrodynamic Instability and Transition to Turbulence. Springer Nether-
lands.
ZANG, T. A. 1991 On the rotation and skew-symmetric forms for incompressible flow simulations.
Appl. Num. Math. 7, 27–40.
77
Appendices
A Detailed perturbation equations
A.1 Linear IVP – forward equations
The forward equations for the first-order perturbation fields are
∂
ˆ
U
∂τ
=
Scμ
0
Z
0
2τ
D
2
ˆ
U+

[(R− 1)μ
0
+ 2βScZ
0
]Y
′
0
+ 2η
4τ

D
ˆ
U (A.1)
−


i(kV
0
+ mW
0
τ)+

k
2
+ m
2

ScKμ
0
Z
0
2
−
(R− 1)

μ
0
Y
′′
0
+ βY
′
2
0

4τ

 ˆ
U
+
ikBμ
0
Z
0
√
τ
D
ˆ
V +
imBμ
0
Z
0
√
τ
D
ˆ
W−
β(R− 1)BZ
0

μ
0
Y
′′
0
+ βY
′
2
0

2
√
τ
3
D
ˆ
Y
+


i

kV
′
0
+ mW
′
0
τ

βBZ
0
√
τ
+
(R− 1)
2
B

μ
2
0
Y
′′′
0
+ 4β μ
0
Y
′
0
Y
′′
0
+ β
2
Y
′
3
0

4Sc
√
τ
3
−
β(R− 1)BZ
0

μ
0
Y
′′′
0
+ 3βY
′
0
Y
′′
0

2
√
τ
3
#
ˆ
Y−
BZ
0
√
τ
D
ˆ
P,
78
∂
ˆ
V
∂τ
=
ikBμ
0
Z
0
√
τ
D
ˆ
U+

ikβBZ
0
Y
′
0
√
τ
−
V
′
0
√
2Kτ

ˆ
U+
Scμ
0
Z
0
4τ
D
2
ˆ
V (A.2)
+

[(R− 1)μ
0
+ βScZ
0
]Y
′
0
+ 2η
4τ

D
ˆ
V−
"
i(kV
0
+ mW
0
τ)+

2k
2
+ m
2

ScKμ
0
Z
0
2
#
ˆ
V
−
kmScKμ
0
Z
0
2
ˆ
W +
βScZ
0
V
′
0
4τ
D
ˆ
Y +
Sc

βZ
0
V
′′
0
−(R− 1)

μ
0
V
′′
0
+ βV
′
0
Y
′
0

4τ
ˆ
Y
−
ikScKZ
0
2
ˆ
P,
∂
ˆ
W
∂τ
=
imBμ
0
Z
0
√
τ
D
ˆ
U+

imβBZ
0
Y
′
0
√
τ
−
r
τ
2K
W
′
0

ˆ
U (A.3)
−
kmScKμ
0
Z
0
2
ˆ
V +
Scμ
0
Z
0
4τ
D
2
ˆ
W +

[(R− 1)μ
0
+ βScZ
0
]Y
′
0
+ 2η
4τ

D
ˆ
W
−
"
i(kV
0
+ mW
0
τ)+

k
2
+ 2m
2

ScKμ
0
Z
0
2
#
ˆ
W +
βScZ
0
W
′
0
4
D
ˆ
Y
+
(
Sc

βZ
0
W
′′
0
−(R− 1)

μ
0
W
′′
0
+ βW
′
0
Y
′
0

4
−
R
2
− 1
2R
)
ˆ
Y−
imScKZ
0
2
ˆ
P,
∂
ˆ
Y
∂τ
=−
Y
′
0
√
2Kτ
ˆ
U+
μ
0
Z
0
4τ
D
2
ˆ
Y +

[(R− 1)μ
0
+ 2βZ
0
]Y
′
0
+ 2η
4τ

D
ˆ
Y (A.4)
−


i(kV
0
+ mW
0
τ)+

k
2
+ m
2

Kμ
0
Z
0
2
+
(R− 1)

μ
0
Y
′′
0
+ βY
′
2
0

− βZ
0
Y
′′
0
4τ

 ˆ
Y,
1
√
2Kτ
D
ˆ
U+ ik
ˆ
V + im
ˆ
W +(R− 1)
(
μ
0
4τ
D
2
ˆ
Y +
βY
′
0
2τ
D
ˆ
Y +
"
βY
′′
0
4τ
−

k
2
+ m
2

Kμ
0
2
#
ˆ
Y
)
= 0,
(A.5)
79
where the symbols D and primes denote η-derivatives. The equations above can be expressed in
compact form as
∂
ˆ
X X X
∂τ
= A A A
ˆ
X X X+ B B B
ˆ
P, (A.6a)
C C C
ˆ
X X X = 0, (A.6b)
where
ˆ
X X X =

ˆ
U,
ˆ
V,
ˆ
W,
ˆ
Y

. Taking derivative of Eq. (A.6b) with respect to time, we have
∂C C C
∂τ
ˆ
X X X+C C C
∂
ˆ
X X X
∂τ
= 0, (A.7)
and substituting Eq. (A.6a) yields the equation for pressure
ˆ
P=−(C C CB B B)
−1

∂C C C
∂τ
+C C CA A A

ˆ
X X X. (A.8)
The pressure can be substituted back into Eq. (A.6a) to obtain
∂
ˆ
X X X
∂τ
=

A A A− B B B(C C CB B B)
−1

∂C C C
∂τ
+C C CA A A

ˆ
X X X≡ L L L
ˆ
X X X. (A.9)
In the numerical implementation, after performing the spatial discretization (as described in
§3.2.1), the implicit second-order Crank-Nicolson time-marching method has the form of
ˆ
X X X
n+1
−
ˆ
X X X
n
Δτ
=
1
2
h
A
n
ˆ
X X X
n
+A
n+1
ˆ
X X X
n+1
i
+B
n+1/2
ˆ
P
n+1/2
, (A.10a)
C
n+1
ˆ
X X X
n+1
= 0. (A.10b)
80
where the block matrices A, B, and C are the discrete counterparts of those in Eq. (A.6). Re-
arranging terms, we have



E F
C 0






ˆ
X X X
n+1
ˆ
P
n+1/2



=



G
ˆ
X X X
n
0



, (A.11)
where E =

I− ΔτA
n+1
/2

, F =

−ΔτB
n+1/2

, G = (I+ ΔτA
n
/2), andI is the identity block
matrix. We apply Schur complement reduction to obtain



E F
0 −CE
−1
F






ˆ
X X X
n+1
ˆ
P
n+1/2



=



G
ˆ
X X X
n
−CE
−1
G
ˆ
X X X
n



, (A.12)
from which we can solve for pressure and subsequently velocity and scalar fields.
A.2 Quasi-steady-state approximation
In the QSSA, the base profiles are assumed to be frozen at a time τ
0
> 0, when they are then
perturbed, and the disturbance quantities have the following forms
(u
1
,v
1
,w
1
,Y
1
, p
1
)=

ˆ
U(η,τ
0
),
ˆ
V(η,τ
0
),
ˆ
W(η,τ
0
),
ˆ
Y(η,τ
0
),
ˆ
P(η,τ
0
)

e
i(ky+mz−ωτ)
+ c.c.,
(A.13)
where ω = ω
R
+ iω
I
, the imaginary part of which indicates temporal perturbation growth rate.
This results in a generalized eigenvalue problem with τ = τ
0
as a parameter (Tan & Homsy, 1986),
which is then solved with the QZ algorithm (Corbett & Bottaro, 2000). The results presented below
are for two-dimensional cases (k= 0) since ω
I,2D
> ω
I,3D
for all cases. Figure A.1(a) displays the
imaginary amplification frequency ω
I
(R,m) for uniform viscosity ratio M = 1 and τ
0
= 10, where
ω
I
increases with increasing R. Next, figure A.1(b) shows ω
I
(M,m) for density ratio R = 5 and
τ
0
= 10, where the largest ω
I
is achieved at the lowest M. Furthermore, figure A.1(c) shows
ω
I
(τ
0
,m) for density ratio R = 5 and uniform viscosity M = 1 as a function of τ
0
. For the largest
81
0.05 0.1 0.15 0.2
-0.1
0
0.1
0.2
0.3
0.4 R = 2.5
R = 5.0
R = 7.5
R = 10.0
m
ω
I
(a)
0.05 0.1 0.15 0.2 0.25 0.3
-0.1
0
0.1
0.2
0.3
0.4 M = 1/4
M = 1/2
M = 1
M = 2
M = 4
m
ω
I
(b)
0.05 0.1 0.15 0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6 τ
0
= 5
τ
0
= 10
τ
0
= 25
τ
0
= 40
τ
0
= 50
m
ω
I
(c)
Figure A.1: Imaginary frequency ω
I
as: (a) function of R for M = 1 and τ
0
= 10, (b) function of
M for R= 5 and τ
0
= 10, and (c) function of τ
0
for R= 5 and M = 1.
τ
0
= 50 case considered here, the peak growth rate is the highest among all, but the range of
wavenumber m at which ω
I
> 0 is the narrowest. The opposite trend is observed at smallest τ
0
,
where the positive growth rate is seen across a wider range of m, but the peak ω
I
is the lowest.
Finally, figure A.2 displays the typical eigenfunctions obtained from QSSA for R = 5, M = 1,
τ
0
= 10, and m= 0.05, which are utilized as initial conditions for the forward integration.
A.3 Linear IVP – adjoint equations
The corresponding adjoint equations are
D

BZ
0
√
τ
+
ˆ
U
+

+
ikScKZ
0
2
V
+
+
imScKZ
0
2
ˆ
W
+
= 0, (A.14)
∂
ˆ
U
+
∂τ
+
= D
2

Scμ
0
Z
0
2τ
+
ˆ
U
+

− D

[(R− 1)μ
0
+ 2βScZ
0
]Y
′
0
+ 2η
4τ
+
ˆ
U
+

(A.15)
+


i(kV
0
+ mW
0
τ
+
)−

k
2
+ m
2

ScKμ
0
Z
0
2
+
(R− 1)

μ
0
Y
′′
0
+ βY
′
2
0

4τ
+

 ˆ
U
+
+D

ikBμ
0
Z
0
√
τ
+
ˆ
V
+

−

ikβBZ
0
Y
′
0
√
τ
+
+
V
′
0
√
2Kτ
+

V
+
+ D

imBμ
0
Z
0
√
τ
+
ˆ
W
+

−

imβBZ
0
Y
′
0
√
τ
+
+
r
τ
+
2K
W
′
0
!
ˆ
W
+
−
Y
′
0
√
2Kτ
+
ˆ
Y
+
−
1
√
2Kτ
+
D
ˆ
P
+
,
82
-10 0 10 20
-0.1
-0.05
0
0.05
0.1
η
ˆ
U
-10 0 10 20
-0.1
-0.05
0
0.05
0.1
0.15
Real
Imag
Abs
η
ˆ
W
-10 0 10 20
-0.01
0
0.01
η
ˆ
Y
-20 -10 0 10 20
0
0.5
1
η
ˆ
P
Figure A.2: The eigenfunctions for R= 5, M = 1, τ
0
= 10, and m= 0.05.
83
∂V
+
∂τ
+
= D

ikBμ
0
Z
0
√
τ
+
ˆ
U
+

+ D
2

Scμ
0
Z
0
4τ
+
V
+

− D

[(R− 1)μ
0
+ βScZ
0
]Y
′
0
+ 2η
4τ
+
V
+

+
"
i(kV
0
+ mW
0
τ
+
)−

2k
2
+ m
2

ScKμ
0
Z
0
2
#
V
+
−
kmScKμ
0
Z
0
2
ˆ
W
+
− ik
ˆ
P
+
, (A.16)
∂
ˆ
W
+
∂τ
+
= D

imBμ
0
Z
0
√
τ
+
ˆ
U
+

−
kmScKμ
0
Z
0
2
V
+
+ D
2

Scμ
0
Z
0
4τ
+
ˆ
W
+

(A.17)
−D

[(R− 1)μ
0
+ βScZ
0
]Y
′
0
+ 2η
4τ
+
W
+

+
"
i(kV
0
+ mW
0
τ
+
)−

k
2
+ 2m
2

ScKμ
0
Z
0
2
#
W
+
− im
ˆ
P
+
,
∂
ˆ
Y
+
∂τ
+
= D


β(R− 1)BZ
0

μ
0
Y
′′
0
+ βY
′
2
0

2
p
(τ
+
)
3
ˆ
U
+


−
"
i

kV
′
0
+ mW
′
0
τ
+

βBZ
0
√
τ
+
(A.18)
−
(R− 1)
2
B

μ
2
0
Y
′′′
0
+ 4β μ
0
Y
′
0
Y
′′
0
+ β
2
Y
′
3
0

4Sc
p
(τ
+
)
3
+
β(R− 1)BZ
0

μ
0
Y
′′′
0
+ 3βY
′
0
Y
′′
0

2
p
(τ
+
)
3

 ˆ
U
+
−D

βScZ
0
V
′
0
4τ
+
V
+

+
Sc

βZ
0
V
′′
0
−(R− 1)

μ
0
V
′′
0
+ βV
′
0
Y
′
0

4τ
+
V
+
−D

βScZ
0
W
′
0
4
ˆ
W
+

+
(
Sc

βZ
0
W
′′
0
−(R− 1)

μ
0
W
′′
0
+ βW
′
0
Y
′
0

4
−
R
2
− 1
2R
)
ˆ
W
+
+D
2

μ
0
Z
0
4τ
+
ˆ
Y
+

− D

[(R− 1)μ
0
+ 2βZ
0
]Y
′
0
+ 2η
4τ
+
ˆ
Y
+

+


i(kV
0
+ mW
0
τ
+
)−

k
2
+ m
2

Kμ
0
Z
0
2
+
βZ
0
Y
′′
0
−(R− 1)

μ
0
Y
′′
0
+ βY
′
2
0

4τ
+

 ˆ
Y
+
+(R− 1)
(
D
2

μ
0
4τ
+
ˆ
P
+

− D

βY
′
0
2τ
+
ˆ
P
+

+
"
βY
′′
0
4τ
+
−

k
2
+ m
2

Kμ
0
2
#
ˆ
P
+
)
.
84
The equations above can be expressed in compact form as
∂
ˆ
X X X
+
∂τ
+
= A A A
+
ˆ
X X X
+
+ B B B
+
ˆ
P
+
, (A.19a)
C C C
+
ˆ
X X X
+
= 0, (A.19b)
where
ˆ
X X X
+
=

ˆ
U
+
,
ˆ
V
+
,
ˆ
W
+
,
ˆ
Y
+

. Taking derivative of Eq. (A.19b) with respect to time, we have
∂C C C
+
∂τ
+
ˆ
X X X
+
+C C C
+
∂
ˆ
X X X
+
∂τ
+
= 0, (A.20)
and substituting Eq. (A.19a) yields the equation for pressure
ˆ
P
+
=−

C C C
+
B B B
+

−1

∂C C C
+
∂τ
+
+C C C
+
A A A
+

ˆ
X X X
+
. (A.21)
The pressure can be substituted back into Eq. (A.19a) to obtain
∂
ˆ
X X X
+
∂τ
+
=

A A A
+
− B B B
+

C C C
+
B B B
+

−1

∂C C C
+
∂τ
+
+C C C
+
A A A
+

ˆ
X X X
+
≡ L L L
+
ˆ
X X X
+
. (A.22)
B Inviscid solution for nonzero velocity difference Δw
0
at t =
t
0
> 0
We generalize the inviscid formulation in §2 to account for a nonzero vertical velocity difference
between the two streams in the z direction, Δw
0
, at t = t
0
> 0, when the infinitesimal perturbation
is initially imposed on the base flow. Prior to that, i.e., when t < t
0
, it is assumed that the base flow
is not perturbed. Denoting± as the right and left fields with respect to the interface, respectively,
the general structure of the equations for perturbation velocities is of the form
∂
ˆ
ϕ
±
∂t
± imw
±
0
(t)
ˆ
ϕ
±
= C
±
0
(t), (B.1)
85
where
w
±
0
(t)=∓2a
±
t, a
±
=
R− 1
4R
±
g, R
+
= R, R
−
= 1, (B.2)
and C
±
0
(t) is the coefficient for the perturbation pressure, whose solution is
ˆ p
±
= C
±
0
(t)e
∓K x
, (B.3)
whereK
2
= k
2
+ m
2
. The relations between
ˆ
ϕ
±
and the perturbation velocities are given by
ˆ u
±
=±
K
ρ
±
ˆ
ϕ
±
e
∓K x
, ˆ v
±
=−
ik
ρ
±
ˆ
ϕ
±
e
∓K x
, ˆ w
±
=−
im
ρ
±
ˆ
ϕ
±
e
∓K x
. (B.4)
Solving Eq. (B.1) with initial time t
0
> 0, we obtain
ˆ
ϕ
±
(t)=
ˆ
ϕ
±
0
e
±ima
±
(t
2
−t
2
0
)
+
Z
t
t
0
C
±
0
(t
′
)e
±ima
±
(t
2
−t
′2
)
dt
′
. (B.5)
Following the formulation in §2, the differential equation for the interface amplitude
ˆ
ζ is given
by
∂
ˆ
ζ
∂t
∓ 2ima
±
t
ˆ
ζ = ˆ u
±
(x= 0)=±
K
ρ
±

ˆ
ϕ
±
0
e
±ima
±
(t
2
−t
2
0
)
+
Z
t
t
0
C
±
0
(t
′
)e
±ima
±
(t
2
−t
′2
)
dt
′

, (B.6)
and for zero surface tension, C
+
0
(t)= C
−
0
(t). Multiplying Eq. (B.6) by e
∓ima
±
t
2
, taking derivative
with respect to time, canceling exponential terms on both sides, and eliminating C
±
0
, we obtain
d
2 ˆ
ζ
dt
2
−
ω
4
0
4
t
2
ˆ
ζ(t)= 0, (B.7a)
with
ω
4
0
= g
2
m
2
(R− 1)
2
R
.
86
Eq. (B.7a) can be non-dimensionalized by introducing ˜ τ≡ ω
0
t, giving
d
2 ˜
ζ
d ˜ τ
2
−
˜ τ
2
4
˜
ζ = 0. (B.8)
The general solution of Eq. (B.8) can be expressed as
˜
ζ( ˜ τ)= C
1
D
−1/2
( ˜ τ)+C
2
D
−1/2
(i ˜ τ), (B.9)
where D
j
(z) denotes the parabolic cylinder function of order j. We require two initial conditions,
the first of which is unity interface amplitude
˜
ζ( ˜ τ = ˜ τ
0
)= 1. The other condition is
˜
ζ
˜ τ
( ˜ τ = ˜ τ
0
)= 0,
obtained from setting t = t
0
in Eq. (B.6) and choosing
ˆ
ϕ
±
0
such that the following condition is
satisfied
∓2ima
±
t
0
=±
K
ρ
±
ˆ
ϕ
±
0
. (B.10)
Finally, the constants are given by
C
1
=
˜ τ
0
D
−1/2
(i ˜ τ
0
)+ 2iD
1/2
(i ˜ τ
0
)
2C
3
, (B.11)
and
C
2
=
˜ τ
0
D
−1/2
( ˜ τ
0
)− 2D
1/2
( ˜ τ
0
)
2C
3
, (B.12)
where
C
3
= iD
−1/2
( ˜ τ
0
)D
1/2
(i ˜ τ
0
)+ D
−1/2
(i ˜ τ
0
)

˜ τ
0
D
−1/2
( ˜ τ
0
)− D
1/2
( ˜ τ
0
)

. (B.13)
87

Abstract (if available)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Pattern generation in stratified wakes

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

PDF

Passive rolling and flapping dynamics of a heaving Λ ﬂyer

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Synchronization in carpets of model cilia: numerical simulations and continuum theory

PDF

Understanding the formation and evolution of boundaries and interfaces in nanostructured metallic alloys

PDF

Near wake characteristics of towed bodies in a stably stratified fluid

PDF

Multiscale modeling of cilia mechanics and function

PDF

Numerical study of focusing effects generated by the coalescence of multiple shock waves

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Development and characterization of hierarchical cellular structures

PDF

The development of a hydraulic-control wave-maker (HCW) for the study of non-linear surface waves

PDF

Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation

PDF

Modeling and simulation of complex recovery processes

Asset Metadata

Creator Prathama, Aditya Heru (author)

Core Title Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Degree Conferral Date 2021-12

Publication Date 09/21/2021

Defense Date 08/27/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag adjoint-based optimization,baroclinic instability,OAI-PMH Harvest,optimal perturbations,shear layers,two fluid columns,variable density

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Pantano-Rubino, Carlos (committee chair), Bermejo-Moreno, Ivan (committee member), Jovanovic, Mihailo R. (committee member), Luhar, Mitul (committee member)

Creator Email aditya.h.prathama@gmail.com,prathama@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC15921173

Unique identifier UC15921173

Legacy Identifier etd-PrathamaAd-10082

Document Type Dissertation

Format application/pdf (imt)

Rights Prathama, Aditya Heru

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu