Mr Kim

Motoyashi Kalland

Sagar Patankar

Palas Kumar Farsoiya

Manpreet Singh (MTech, IIT Bombay) and Aishwarya Nair (Summer intern, IIT Bombay

Table of Contents Introduction

Surface waves "...that are easily seen by everyone and which are usually used as an example of waves in elementary courses, are water waves. As we shall soon see, they are the worst possible example, because they are in no respects like sound and light; they have all the complications that waves can have."

-Feynman et al. (1964)

Literature Review

The motion of an interface between two immiscible fluids is ubiquitous. The most common of these viz. the air-water interface in oceans spans 71% of the surface area of the earth. Breaking of capillary-gravity waves on the ocean surface play an essential role in the generation of spray from the ocean [START_REF] Veron | Ocean spray[END_REF] Interplay between surface waves and a turbulent velocity field (in air as well as water) is an an important component in fundamental studies of wind-wave interactions [START_REF] Janssen | The interaction of ocean waves and wind[END_REF]. Interfacial oscillations can also be found in raindrops where close to spherical water droplets, oscillate as they fall (Szakáll et al., 2010).

Cylindrical falling jets from faucets undergo oscillations or the Rayleigh-Plateau instability during their motion. Such oscillations have been of interest since the early 19th century. An important analysis for interfacial waves was carried out by Cauchy and Poisson (Cauchy, 1827;[START_REF] Poisson | M {é} moire sur la th {é} orie des ondes[END_REF]. They independently solved the linearised, initial-value problem for inviscid, irrotational waves with arbitrary (Fourier transformable) initial conditions. Their solution comprised of superposition of standing waves excited due to initial conditions [START_REF] Craik | The origins of water wave theory[END_REF].

Waves in the real world are viscous with boundary layers at the interface and additionally can be affected by the presence of a vorticity field (e.g. due to turbulence). In the near two hundred years that have elapsed since the landmark studies on waves by Cauchy and Poisson, the subject of interfacial waves has expanded with myriad applications in environmental sciences and engineering.

In addition to free oscillations, the subject of parametrically forced oscillations and Faraday instability has also seen intense, focussed activity. Of particular interest has been engineering applications where atomisation and spray generation for medical utility has been proposed using the Faraday mechanism (Tsai et al., 2017). From a fundamental perspective, recent results have demonstrated the possibility of generating Faraday waves on a sphere (Adou and Tuckerman, 2016). In the literature survey that follows, we provide an overview of the free and parametrically forced oscillation literature. While we study capillary-gravity waves in this thesis, note that the inverse problem of a capillary-gravity wave is the Rayleigh-Taylor instability when the heavier fluid lies above the lighter one. The dispersion relation is obtained from the same equation 1.1 with densities are interchanged (See Table 1.1h-1.1i).

Free oscillations of the interface

Laplace (1776) posed the general initial value problem (IVP hereafter) for waves on an interface. Given an initial surface disturbance on a deep pool, how does the interface evolve in time? The French Academy of Sciences announced this as a prize problem in December, 1813. Augustin-Louis Cauchy submitted his linearised analysis (Cauchy, 1827). One of the judges, Simeéon D. Poisson also independently solved the linearised, problem and submitted an entry [START_REF] Poisson | M {é} moire sur la th {é} orie des ondes[END_REF]. The solutions of Cauchy and Poisson were both correct but restricted to the linearised, inviscid, irrotational assumption. They obtained solutions in cylindrical axisymmetric as well as two-dimensional Cartesian coordinates (See Table 1.1a-1.1d ). Their solution comprised of linear superposition of modes excited by initial conditions (viz. the Fourier transform, before its usage became popular [START_REF] Craik | The origins of water wave theory[END_REF])) expressing the shape of the interface at any time t as a Fourier integral. Since the seminal work of Cauchy & Poisson it has become customary in the literature to refer to the solution to an IVP for interfacial/free surface waves as a Cauchy-Poisson problem (Debnath, 1994;[START_REF] Ellingsen | Initial surface disturbance on a shear current: The cauchypoisson problem with a twist[END_REF].

Despite the insight of Cauchy & Poisson in their landmark analysis, the restrictive nature of the solution is clear. For instance, if one wishes to understand the effect of a finite region of vorticity on wave generation on the interface, it is clear that Cauchy & Poisson's solution is not of much interest (see [START_REF] Prosperetti | Small-amplitude waves produced by a submerged vorticity distribution on the surface of a viscous liquid[END_REF] for linearised analysis). From a linearised perspective, solving the inviscid Euler's equation with the vortical initial condition is also of not much utility as the linearised equation for perturbation vorticity is ∂w ∂t = 0 which predicts that the perturbation vorticity field does not evolve in time. Waves in the real world are viscous and are affected by (and in turn effect) vortical regions around them. The simplest non-trivial model which describes this is a linearised viscous model (although interactions governed by nonlinear effects are also crucially important). One of the first viscous analysis was by Lamb (1932) who obtained the viscous dispersion relation in his textbook. Interestingly, Lamb also anticipated that the temporal spectrum of the viscous problem (in deep water) also contained the continuous spectrum in addition to the discrete spectrum (see discussion around equation 26 in [START_REF] Lamb | Hydrodynamics[END_REF], chapter XI, article 349). The first analysis of the Cauchy-Poisson linearised initial-value problem for standing waves on deep water (with zero vorticity initial conditions) was by [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF]. In his seminal study, Prosperetti (1976) solved the IVP for the interface deformed as a single Fourier mode. While the connection to normal mode analysis was not explicitly discussed, Prosperetti (1976) showed that the amplitude a(t) for a linearised standing wave could be written as a sum of discrete exponentials and a set of other terms (error functions). The former arose from the discrete spectrum while the latter was the continuous spectrum contribution. For axisymmetric perturbations on a radial geometry, Miles (1968) solved the viscous Cauchy-Poisson problem for free surface waves. The viscous dispersion relation at the interface of two fluids was obtained by Chandrasekhar (1961a) but did not discuss the continuous spectrum. In Chapter 3 we show that the continuous spectrum also exists for a cylindrical shaped deep pool. In particular, we also show that the cylindrical pool is entirely analogous to the rectangular pool in the linearised framework. Knowledge of the temporal spectrum for this problem constitutes an important first step in modeling the effect of vortical regions of disturbances on the propagation of waves on the deep ocean.

Similar to the study of waves on a deep pool of water, study of perturbations on a cylindrical fluid filament constitutes a classical problem in interfacial wave theory. For a cylindrical filament of radius R 0 (taken to be infinite in length), the Rayleigh-Plateau instability is well known (Plateau, 1873)(Table 1.1e-1.1g). From inviscid, irrotational linearised analysis it can be shown that axisymmetric modes with kR 0 < 1 are unstable whereas any three dimensional or axisymmetric perturbation with kR 0 > 1 are stable and produce oscillations. For perturbations on a spherical free surface (droplet or bubble) Kelvin (1890) and Rayleigh (1896) used inviscid irrotational theory to obtain the dispersion relation with the restoring force being either a radial body force and surface tension respectively(See Table 1.2a-1.2b) Purely azimuthal perturbations on a cylindrical filament were studied by Rayleigh and (Bohr, 1909) )with the aim of using the dispersion relation to estimate surface tension. Here too, viscous effects may be significant. The dispersion relation governing small amplitude, axisymmetric perturbations on a viscous liquid cylinder was first obtained by Rayleigh (1892b) (also see [START_REF] Bauer | Natural damped frequencies of an infinitely long column of immiscible viscous liquids[END_REF]; [START_REF] Liang | Linear stability analysis of capillary instabilities for concentric cylindrical shells[END_REF]; [START_REF] Stone | Note on the capillary thread instability for fluids of equal viscosities[END_REF]; [START_REF] Tomotika | On the instability of a cylindrical thread of a viscous liquid surrounded by another viscous fluid[END_REF] for further results). Many microfluidic applications involve a gas inside the filament whose density and viscosity may be negligible compared to the liquid outside. We call this a "hollow filament". For axisymmetric perturbations on a hollow filament surrounded by a viscous liquid, the dispersion relation has been obtained in Chandrasekhar (1981) (also see [START_REF] Bauer | Natural damped frequencies of an infinitely long column of immiscible viscous liquids[END_REF]). For purely azimuthal perturbations, the viscous dispersion relation is not known in the literature and is derived in this thesis. Analogous to the deep water wave problem discussed in the earlier paragraph, we show in Chapter 4 that the temporal spectrum for a hollow filament as obtained from discrete normal modes is not complete. The analysis here is rendered complicated by the circular geometry which causes special functions (modified Bessel functions) to appear. However, the fundamental idea still remains the same. Discrete normal mode eigenfunctions arising from the dispersion relation decay exponentially in the far field (i.e. at r ! •) and cannot represent vortical disturbances there. A continuum of eigenfunctions (Bessel J and Y) are also found to exist which are necessary for adequately representing such vortical initial conditions. The solution to the IVP validates what we anticipate from a careful normal mode analysis. We also refer the reader to the studies by [START_REF] García | Normal-mode linear analysis and initial conditions of capillary jets[END_REF] and [START_REF] Berger | Initial-value stability analysis of a liquid jet[END_REF] for some conclusions on the related problem of studying the temporal spectrum of a liquid cylinder. Table 1.2 lists some of the important analysis relevant to our study highlighting the novel contributions of this thesis.

Forced oscillations of the interface

The study of parametrically forced oscillations started with Faraday (1831b) who reported observing wave patterns on mercury, ink, water, alcohol, turpentine, milk and white of egg covering a vibrating horizontal plate. He remarked that these waves have frequency half of that of excitation. These waves are often referred to as Faraday waves in the literature. [START_REF] Miles | Parametrically forced surface waves[END_REF]. [START_REF] Benjamin | The stability of the plane free surface of a liquid in vertical periodic motion[END_REF] reported the first theoretical analysis under a linearised, irrotational inviscid approximation and showed that the amplitude of a standing wave is governed by the Mathieu equation. (h) Geoffrey Ingram Taylor (1886Taylor ( -1975) ) Source: Biographical Memoirs (Batchelor, 1976) (Reproduced under fair use policy for educational use)

n 2 = K(g + g 1 ) r 1 r 2 r 1 + r 2 (1.1)
F i = A cos ( t s g a i ✓ 1 3 2i + 1 ◆ E )
Equation 99 "Oscillations of a liquid sphere" Kelvin (1890), where i and g are spherical harmonic and body force in radial direction respectively.

p 2 = n(n 1)(n + 2) T ra 3 Eq. 9 on page 373 in "The theory of sound" Rayleigh (1896), where n is spherical harmonic and T is surface tension coefficient.

(c) Niels Bohr (1885Bohr ( -1962) ) Source: Niels Bohr -Facts. NobelPrize.org. ©Nobel Media AB 2018. accessed: 27 May 2019 (Reproduced with special permission for educational use)

(d) azimuthally perturbed jet cross section

q 2 = T r 1 + r 2 n 3 n a 3
where q is angular frequency of oscillation as in page 297 in "Determination of the Surface-Tension of Water by the Method of Jet Vibration" Bohr (1909). Normal mode analysis and solution to IVP including inertia of both fluids and zero initial vorticity, explaining the presence of the continuous spectrum and its contribution, interfacial Cauchy-Poisson problem Prosperetti (1976, 1981), Prosperetti (1980a,b) Free surface, as well as interfacial initial value problem in Cartesian geometry. Spherical free surface as well as interfacial normal mode analysis and IVP.

Prosperetti and Cortelezzi (1982)

Initial value problem for non zero initial vorticity distribution for Cartesian free surface of viscous fluid [START_REF] Berger | Initial-value stability analysis of a liquid jet[END_REF] Initial value problem for axisymmetric perturbations on a cylindrical fluid filament solving for fluid inside only.

García and González (2008)

Normal mode analysis claiming absence of continuous spectrum for a viscous cylindrical filament solving for fluid inside only.

Present Thesis and Farsoiya et al. (2019)

Normal mode analysis for non-axisymmetric two dimensional perturbations for a cylindrical filament solving for outer fluid only. Shows and explains the presence of the continuous spectrum from normal mode analysis as well as IVP approach with zero initial vorticity. They compared their linearised theory with experimental results which confirmed their stability predictions (see figure 1.1). The small discrepancy between the theory and experiments was attributed to wetting effects.

Kumar and Tuckerman (1994) have studied the effect of viscosity on these oscillations concluding that the equation governing the amplitude cannot be simply reduced to Mathieu equation with linear damping. Recent studies by this group have extended the analysis showing that similar Faraday waves can arise on a sphere. Their study is supplemented with linear theory (Adou and Tuckerman, 2016) and numerical simulations (Ebo- [START_REF] Ebo-Adou | Faraday instability on a sphere: numerical simulation[END_REF]. An extensive literature exists on analysing Faraday waves using linear and weakly nonlinear theory (see reviews by Miles and Henderson (1990) and Perlin and Schultz (2000)). We note in particular recent work by Tsai et al. (2017) which has focussed on applications of Faraday waves related to drug delivery. In Chapter 5, we show the existence of Faraday waves on a cylindrical filament subjecting to a radially pulsating body force. Our theory is linearised, inviscid and irrotational and shows excellent agreement with DNS for small forcing amplitude. We also suggest potential applications of these Faraday waves.

Thesis overview

The aim of this thesis is to study a set of problems concerning viscous and inviscid, free and forced interfacial waves. We demonstrate that the oscillations of interfacial waves in various geometrical configurations (See figure 1.2a-1.2c) have commonalities between them and are connected through a common thread. The author has also developed a two dimensional Cartesian and cylindrical axisymmetric, incompressible Navier-Stokes solver in C, based on the Volume of Fluid (VoF) method [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. This solver has been used for validating theoretical results presented in this thesis. This thesis is organised as follows:

In chapter 2, we briefly discuss the algorithms used for developing the Navier-Stokes-VoF solver and present results for benchmarking the same. The VoF module was developed earlier by Dasgupta (2010) and is benchmarked, extended and coupled to a Navier-Stokes solver developed by the author. The axisymmetric part of the solver is implemented together with Ms. Aishwarya Nair (IIT Bombay Research Intern from BITS Pilani, Hyderabad Campus, 2017). Chapter 3 investigates capillary-gravity standing, free oscillations. Here we show the presence of the continuous spectrum for linearised, capillary waves on a cylindrical pool of infinite depth. The discrete spectrum coming from the dispersion relation, is found to provide only two eigenfunctions. We solve the initial value problem for an interfacial perturbation, showing that this choice of initial condition has projections on the discrete as well as the continuous spectrum eigenfunctions. The consequence of this is that, the time evolution of the interface cannot be represented as a sum of discrete exponentials (in time) alone. Our analytical results show very good agreement with Direct Numerical Simulations and extend earlier classical results of Prosperetti (Prosperetti, 1980a[START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF] in Cartesian and spherical geometries to interfacial oscillations occurring on a cylindrical, radially unbounded pool of infinite depth. In Chapter 4, we investigate the continuous spectrum for azimuthal perturbations on a gaseous filament surrounded by a liquid. Using a normal mode approach, we show the existence of the discrete and continuous spectrum for this problem. An initial value problem is solved and shows good agreement with DNS. This analysis extends unites spherical, Cartesian and cylindrical geometries into a single framework. Chapter 5 investigates and shows the possibility of Faraday waves on an inviscid cylindrical fluid filament surrounded by another inviscid fluid. A Mathieu equation governing the amplitude of parametrically forced standing waves on a cylindrical filament, is derived in the linearised approximation. Predictions from this equation show good agreement with DNS. The possibility of stabilising the Rayleigh-Plateau instability is examined and compared against DNS. This study was done in collaboration with Mr. Sagar Patankar (PhD student, Dept. Chemical Engineering, IIT Bombay, Patankar et al. ( 2018)). Chapter 6 investigates the effect of density ratio on inviscid, irrotational interfacial oscillations identifying the regimes where a free surface approximation provides acceptable accuracy. We discuss linearised analytical results for waves occurring on base states described by Cartesian and plane polar geometries, comparing these with results obtained from two in-house developed codes i.e one fluid solver (OFS) and free surface solver (FSS). This study was done in collaboration with Mr. Manpreet Singh (MTech, Dept. Chemical Engineering, IIT Bombay,

Singh et al. (2019a) ).

Chapter 7 concludes the thesis with a summary of results obtained and future directions of investigation.

Two-phase volume of fluid based Navier-Stokes solver

We acknowledge contributions from Ms. A. Nair (IIT Bombay Research intern) in writing some parts of the axisymmetric solver.

In this chapter we present the development of an in-house code for simulating two-phase flows in two-dimensional axisymmetric and planar geometry. Engineers, mathematicians and physicists have studied the dynamics of multiphase flows for many years now. Such flows have myriad applications in the industry and are also interesting from a fundamental science perspective. In addition to the governing equations being nonlinear, the position of the interface needs to be found out as a part of solution. Linearised analytical solutions exist for problems where the base state velocity is zero such as oscillations of bubbles and droplets, linear waves or in the limit of low Reynolds numbers such as steady-state motion of bubbles and droplets in Stokes flow [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. In many practical applications involving such multiphase flows, numerical solutions are necessary to access regimes where linearised predictions may not be sufficient. A need for robust algorithms which include surface tension and can handle strong deformation of the interface including breakup and merger has been felt in the multiphase research community since late 50's and 60's. Four class of methods broadly exist in the current literature [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF].

1. Front capturing -Marker and Cell (MAC), Volume Of Fluid (VOF), Level Set (LS), Constrained Interface Profile (CIP) and Phase Field (PF)

Disadvantages:-Volume of fluid method has a geometric advection algorithm 11 which is fairly complex to implement (especially in three-dimensions) while Level Set methods have traditionally suffered from mass conservation issues.

2. Boundary fitted grids-uses separate, boundary-fitted grids for each phase.

e.g. The axisymmetric bubble simulation by Ryskin and [START_REF] Ryskin | Numerical solution of free-boundary problems in fluid mechanics. part 1. the finite-difference technique[END_REF] and three-dimensional bubble by Takagi and Matsumoto (1994)

Disadvantages:

-Increased intricacies in code for complex geometries. During the evolution of the interface, especially near breakup scenarios, the grid may become very skewed and fresh grid generation will be required. Also, numerical clipping of the interface is an artificial method that would be required to capture topological changes. Disadvantages:-Complex approach in axisymmetric and other non-planar geometries.

4.

Front tracking -where separate front is represented by connecting marker points which are directly advected.

Disadvantages:-Changes in topology, where fluid regions coalesce or break up are not performed implicitly. It can be done but at the expense of increased code complexity. Quite complex to maintain connectivity of the marker particles (or the grid point for the lagrangian representation of the interface). Also, implementation in parallel is not efficient since the entire interface needs to be dealt with in a single processor.

In 1970's Noh and Woodward (1976) came up with the idea of Simple Line Interface Calculation (SLIC) to approximate the free boundaries in numerical simulations as straight lines parallel to one of the coordinate axes. [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] improved upon this introducing the Volume of Fluid (VOF) method based on the concept of fractional volume of fluid. This algorithm is now widely used due to its flexibility and efficiency compared to other methods. The VOF method conserves mass up to machine accuracy and is one of the strengths of the algorithm. There are some hybrid methods which use the strengths of level set and VoF methods such as coupled level set VoF (CLSVOF) [START_REF] Gerlach | Comparison of volumeof-fluid methods for surface tension-dominant two-phase flows[END_REF].

All two dimensional simulations reported in this thesis have been carried out using an in-house developed code based on the finite volume method for solving the incompressible Navier-Stokes equations. The interface is captured using the VOF algorithm [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF][START_REF] Noh | Slic (simple line interface calculation)[END_REF]. At the start of PhD, a VOF based module developed earlier by Dasgupta (2010) was inherited. This uses the LVIRA (Puckett et al., 1997) reconstruction algorithm. Following the same structure as the open source code Gerris [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF], the solver is written in a manner such that it can be changed from cylindrical axisymmetric coordinates to two dimensional Cartesian coordinates, using a single flag. We follow the one fluid approach in VOF [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF], wherein the code necessarily solves the equations of motion in two fluid phases, in contrast to the free surface approach (Popinet and Zaleski, 2002; Tan, 2016) . Surface tension is added as a body force in the momentum equation following [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. We explain the VOF methodology next.

Volume of Fluid method

The VOF method starts with defining a quantity called volume fraction f (also sometimes called color function), for each cell in the computational domain. For a two-phase system there is a dark fluid and light fluid with f defined such that it is the ratio of volume of the dark fluid in a computational cell, to the volume of the cell itself. Cells with only light fluid thus have f = 0, those with only dark fluid f = 1 and interfacial cells having 0 < f < 1. The volume of fluid when advected does not change with respect to the fluid parcel (Hirt and Nichols, 1981). Hence, it satisfies the following relation,

∂ f ∂t + (u • r) f = 0 (2.1)
For incompressible flow, the Eulerian form of 2.1 would be,

∂ f ∂t + r• (u f ) = 0 (2.2)
where u u u is underlying velocity field. The algorithm consists of two basic steps:-

1. Interface Reconstruction 2.
Advection of f field

Interface Reconstruction

One of the reasons for the high accuracy of the VOF method is due to its usage of geometrical techniques to reconstruct and advect the interface. We follow the piecewise-linear approach (PLIC) [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Youngs | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF], where a straight line consistent with the value of f is plotted inside each computational cell where 0 < f < 1. Note that a line may be reconstructed in many ways due to the fact that there are infinite choices of straight lines for a prescribed (local) value of f . In order to choose the line uniquely, further constraints are necessary. A logical choice is to use information of the discrete f field from the neighboring cells. We choose the Least Squares Volume of Fluid Reconstruction algorithm proposed by (Pilliod Jr and Puckett, 2004) which is considered one of the most accurate algorithms available in the literature. Please note that the extension of this method to three dimensions is computationally expensive and impractical. The steps of this method are described below:

Least Squares Volume of Fluid Interface Reconstruction Algorithm (LVIRA) (Puckett

et al., 1997)

As mentioned above, our first task in VOF is to reconstruct the interface from a given discrete f field by fitting a straight line inside (every) computational cell labeled (r, c) (r, c 2 Z + ), satisfying 0 < f r,c < 1 (See Fig. 2.1). The choice of the fitted line in LVIRA is guided by the additional constraint that when this line is extrapolated into a 3 ⇥ 3 stencil centered on the (r, c) cell, the difference between the given f field (for the 3 ⇥ 3 stencil) and the f field obtained from linear extrapolation, is minimum for all choices of straight lines. The detailed sub-steps in interface reconstruction in LVIRA are described below:

Step I : Obtain F field

The first step is obtaining an f -field by knowing the initial interface and obtaining from it, the f for each computational cell, see figure 2.1a.

Step II : Initial guess of slope

In order to fit a straight line into a computational cell, we require a guess for the slope of the line. A guess for slope of the normal to the interface (defined as pointing away from the dark fluid towards the lighter fluid) is obtained using (Green-Gauss gradient (Gerlach et al., 2006)) given by

N R = 1 Dx [ f r+1,c+1 + 2 f r,c+1 + f r 1,c+1 f r+1,c 1 2 f r,c 1 f r 1,c 1 ] (2.3) N z = 1 Dy [ f r+1,c+1 + 2 f r+1,c + f r+1,c 1 f r 1,c+1 2 f r 1,c f r 1,c 1 ] (2.4)
where N R and N z are the horizontal and vertical components of the normal. Step III: Quadrant Identification

To get the interface as a line in two-dimensional, we have to determine its orientation inside the cell. The signs of N R and N z given by 2.3 and 2.4 determine the quadrant (first, second, third or fourth) into which the normal points.

Step IV: Angle (q) calculation

The angle, q is the absolute acute angle made by the interface with x-axis. All other quadrants are treated as the 1 st quadrant by suitable rotation. For the second and the fourth quadrant q has to be redefined. Refer to Figure 2.2 for definition of q.

Step VI: Perpendicular distance (p) calculation

Given a volume fraction and the slope of the interface in a cell the perpendicular distance from the left bottom corner of the cell can be calculated using the reverse relations (Scardovelli and Zaleski, 2000) as explained in the next section.

Step VII : Line extrapolation After, reorientation of the line to the first quadrant. The equation of the line with angle q and perpendicular distance p is, (see figure 2.2)

x sin q + y cos q = p (2.5) Step VIII: Norm Calculation

Let the new volume fractions in the extrapolated cells be e f r,c . A norm is defined as,

L 2 (q) = 1 Â k,l= 1 ( e f r+k,c+l f r+k,c+l ) 2 (2.6)
where e f are the new volume fractions (due to extrapolation) while f is the (given) volume fraction with e f r,c = f r,c .

Step IX : Norm Minimization Note that L 2 is a function of the angle q (see figure 2.3). The LVIRA algorithm now repeats all the above steps for an incrementally different q + dq until a local minimum of L 2 is found. The q corresponding to this minimum is accepted. A sample circle reconstructed with the LVIRA is shown in Fig. 2.1b

Advection

One of the advantages of the VOF technique is good mass conservation. This is primarily due to geometrical techniques used for reconstruction and advection. Before we described the VOF advection algorithm, we show the demerits involved in solving 2.7 equation using a regular central difference scheme. Given an f field, equation 2.7 is solved by central differencing and the resulting f field is plotted at every instant using the VOF reconstruction algorithm outlined above. As seen in figure 2.4, we initialise an annulus with f = 1 on the annulus and 0 otherwise (Rudman, 1997). As seen in the figure, the annulus quickly loses the sharpness of the interface present at t = 0. We will see that the VOF method (with its geometric technique (See Fig. 2.5) for solving equation 2.7), performs significantly better in this test case.

∂ f ∂t + r• (u f ) = 0 (2.7)

Geometric flux calculation

In order to calculate fluxes geometrically we solve forward and reverse problem following the nomenclature by [START_REF] Scardovelli | Analytical relations connecting linear interfaces and volume fractions in rectangular grids[END_REF]. Forward problem is defined as calculation of area under a line given the perpendicular distance and orientation of line. Reverse problem is to calculate perpendicular distance given the area under the line and its orientation (Scardovelli and Zaleski, 2000).

We described the details below following the algorithm provided in (Scardovelli and Zaleski, 2000):

As shown in the Fig. 2.6, given a rectangle with sides c 1 and c 2 and a line whose equation is

m 1 x + m 2 y = a 0 (2.8)
The area under the line is (Scardovelli and Zaleski, 2000)

A = 1 2m 1 m 2 h a 02 H(a 0 m 1 c 1 ) a 0 m 1 c 1 2 H(a 0 m 2 c 2 ) a 0 m 2 c 2 2 i
(2.9) where H is the Heaviside function. The above equation can be thought of as:

Area[ADGFBA] = Area[HAE] Area[BFE] Area[HDG]
We differentiate equation 2.9 with respect to a 0 and equate to zero to obtain,

2a 0 2(a 0 m 1 c 1 ) 2(a 0 m 2 c 2 )) 2m 1 m 2 = 0
which gives us the maximum value of a 0 ,

a 0 max = m 1 c 1 + m 2 c 2
The equation of line in 2.8 can therefore be normalised by using this value to result in the following.

m 1 (m 1 c 1 + m 2 c 2 ) x + m 2 (m 1 c 1 + m 2 c 2 ) y = a 0 (m 1 c 1 + m 2 c 2 ) m 1 (m 1 c 1 + m 2 c 2 ) x + m 2 (m 1 c 1 + m 2 c 2 ) y = a 0 a 0 max ) m c 1 x + (1 m) c 2 y = a (2.10)
where m ⌘ m 1 c 1 m 1 c 1 +m 2 c 2 and a ⌘ a 0

(m 1 c 1 +m 2 c 2 ) . Note that m 1 and m 2 are positive nos. and 0  m  1. Similarly because 0  a 0  (m 1 c 1 + m 2 c 2 ), hence 0  a  1. In terms of a and m, the area equation can be written as Forward and reverse relations over a range from 0 to D

A = c 1 c 2 2m(1 m) [a 2 H(a m)(a m) 2 H(a (1 m))(a (1 m)) 2 ]

Forward relations -Case 1:

For case 1, 0  a  0.5 and 0  m  0.5, the following two possibilities arise Case a:

0  a  m  0.5  1 m A = (c 1 c 2 ) a 2 2m(1 m) Case b: 0  m  a  0.5  1 m A = (c 1 c 2 ) 2a m 2(1 m)

Forward relations -Case 2:

For case 2, 0.5 < a  1.0 and 0  m  0.5, the following two possibilities arise

Case a:

0  m  0.5 < a  1 m A = (c 1 c 2 ) 2a m 2(1 m) Case b: 0  m  0.5  1 m  a A = (c 1 c 2 ) ✓ 1 (a 1) 2 2m(1 m)

◆

Forward relations -Case 3:

For case 3, 0  a  0.5 and 0.5 < m  1.0, the following two possibilities arise Case a:

0  a  1 m < 0.5 < m A = (c 1 c 2 ) a 2 2m(1 m) Case b: 0  1 m  a  0.5 < m A = (c 1 c 2 ) 2a 1 + m 2m
Forward relations -Case 4:

For case 4, 0.5 < a  1.0 and 0.5 < m  1.0, the following two possibilities arise Case a:

0  1 m < 0.5 < a  m A = (c 1 c 2 ) 2a 1 + m 2m
Case b:

0  1 m < 0.5 < m  a A = (c 1 c 2 ) ✓ 1 (a 1) 2 2m(1 m)

◆

For calculating perpendicular distance as required in reconstruction step-IV, we need the reverse relations. Note that for writing the reverse relations we assume

c 1 = c 2 = D.

Reverse relations Case 1

Case a:

0  m  0.5  1 m and 0  A  (D 2 ) m 2(1 m) a = p 2Am(1 m) D Case b: 0  m  0.5  1 m and (D 2 ) m 2(1 m) < A  0.5D 2 a = A(1 m) D 2 + m 2 Case 2
Case a:

0  m  0.5  1 m and 0.5D 2 < A  (D 2 ) 2 3m 2(1 m) a = A(1 m) D 2 + m 2
Case b:

0  m  0.5  1 m and (D 2 ) 2 3m 2(1 m) < A  D 2 a = 1 s 2m ✓ 1 A D 2 ◆ ( 1 m)

!

Since a  1, so other root is not considered.

Reverse relations -Case 3

Case a:

0  1 m < 0.5 < m and 0  A  (D 2 ) 1 m 2m a = p 2Am(1 m) D Case b: 0  1 m < 0.5 < m and (D 2 ) 1 m 2m < A  0.5D 2 a = Am D 2 + 1 m 2
Reverse relations -Case 4

Case a:

0  1 m < 0.5 < m and 0.5D 2 < A  (D 2 ) 3m 1 2m a = Am D 2 + 1 m 2
Case b:

0  1 m < 0.5 < m and (D 2 ) 3m 1 2m < A  D 2 a = 1 s 2m ✓ 1 A D 2 ◆ (1 m)

!

Since a  1, other root is not considered. The reverse relations are used to calculate the perpendicular distance. Area is calculated using forward relations [START_REF] Scardovelli | Analytical relations connecting linear interfaces and volume fractions in rectangular grids[END_REF].

(

Volume Fraction Calculation

The fluxes calculated in the previous section are used to update the volume fractions in the cells for the next time step. For this we use Direction Split

Young's method (DSY). [START_REF] Youngs | Fully multidimensional flux-corrected transport algorithms for fluids[END_REF]. In this method the order of directions is interchanged after each time step in order to avoid systematic errors. To achieve mass conservation, for this the basic condition is velocity divergence. In this method, the mass cannot be conserved until all directions are taken into consideration. Hence the values of volume fraction larger than unity and less than zero may arise after first sweep which violates the restriction, 0 6 F 6 1 for the next sweep. This difficulty can be solved by introducing effective volume of the cells ((Rudman, 1997)). But this induces risks of small under or overshoots. Undershoot occurs when the all the volume fraction has to be fluxed out from the cell but it does not and cell cannot be emptied, overshoot occurs when the cell in fluxed with the volume fraction to get filled but it does not. Both occurs because of use of effective volume in calculation of fluxes through wall. This can be resolved by using flux correction in the second sweep (Lörstad and Fuchs, 2004). The following algorithm is used to calculate the volume fractions for the 2D flow.

Effective volume(non-dimensional, characteristic area D 2 ) for the I st sweep,

dV I r,c = 1 Dt V out R out V in R in R c D (2.12)
where, dV I r,c is the effective volume after I st sweep, R out , R in and R c are the radial metric. The fluxes J can be calculated for all the directions using forward relations. Net Flux out DJ out , (non-dimensional) in or out in a cell is given by,

DJ out = (J out R out J in R in ) R c D 2 (2.13)
From which Volume Fraction after I st sweep of the cell is given by, 2.14) where, F 0 and F I are the volume fractions before I st sweep and after I st sweep respectively.

F I r,c = F 0 DJ out dV 1 r,c (
For I I nd sweep,

U ⇤ out = Dt D u out (2.15) z = 8 < : z 1 = J out U ⇤ out if|z 1 1 2 | > |z 2 1 2 | z 2 = F r,c J out 1 U ⇤ out if|z 1 1 2 | < |z 2 1 2 | (2.16) Corrected volume, dV corr = dV I r,c + (1 dV I r,c )z (2.17)
Corrected outgoing flux,

J corr out = dV I r,c (J out U ⇤ out ) + dV corr U ⇤ out (2.18)
Corrected Total flux,

DJ corr T = J corr out r out J in r in r c (2.19)
Final volume fraction after last sweep is given by,

F I I r,c = F I r,c dV I r,c DJ corr T (2.

Verification

The VOF solution methodology for advancing 2.7 is solved employing the above forward and reverse relations.

Circle in translational flow

In this test the given velocity field has zero gradient and is the simplest among all the tests. An annulus (with f = 1 inside and f = 0 outside is advected with this field solving equation 2.7). The f field is advected by two velocity fields

u = 1, v = 0 and u = 2, v = 1.
The number of grid points on the domain [0,4] x [0,4] is 200 ⇥ 200. The Courant number is 0.25 and advection proceeds for 500 steps.

Advection test for solid body rotation

The test has rate of strain tensor zero and a uniform vorticity tensor (constant angular velocity). A slotted circle Zalesak (1979), with center at (2.0, 2.75) and unit diameter is initialised. The length and width of slot are 0.5 and 0.12 respectively, see figures 2.9 and 2.10. The domain is [0,4] x [0,4] with 200x200 grid. The velocity field is given by u = W(y y 0 ), v = W(x x 0 ), where axis of rotation (normal to the plane of the paper) passes through the (x 0 , y 0 ). The angular velocity W = 0.5 and (x 0 , y 0 ) is (2,2) with Dt = 0.005.

Calculation of error

In order to quantify error in the above results, we define where F n is the solution of volume fraction field after n time steps of computation, F e is the exact solution, and F 0 is the initial solution. The initial solution can be calculated by the initial volume fraction field, exact solution of field for translational velocity fields can be easily calculated by recreating the circle at the center which has moved with the velocity field. For solid body rotation the exact solution is equal to the initial solution after one full rotation. For the shear test after 1000 backward steps the final solution should also be equal to initial solution. The errors calculated for various tests are shown in Table 2.1 using data obtained from Rudman (1997), including data from in-house developed code in the last two columns. Note that D f ⌘ 100 ⇥

E = Â |F n r,c F e r,c | Â F 0 r,c
f f f i f i
, is a measure of gain/loss of volume fractions.

Numerical solution of Navier-Stokes equations

We now describe, implementation of an in-house code which solves the incompressible Navier-Stokes equations and is integrated with the VOF solver for equation 2.7.

∂u ∂t + r • (uu) = 1 r rp + 1 r ⇣ r • n µ ⇣ ru + ru T ⌘o⌘ + F b (2.23) F b = g + T kd s n (2.24)
where u is the velocity field, p is the pressure field, g is acceleration due to gravity, T ⌘ T/r where T is the coefficient of surface tension and r the density, k is inverse of the local radius of curvature at the interface, n is the local unit normal to the interface, d s is the surface delta function, g is the acceleration due to gravity and f is the volume fraction field. The density r and viscosity µ are defined as

r = f r U + (1 f )r L and µ = f µ U + (1 f )µ L . Equations 2.22-2.
24 are integrated over a computational cell and the resulting equations are discretised and solved on a staggered grid. Note that the equation 2.23 is written for an approximation that r is constant over velocity control volumes. We define average velocity field over a control volume as,

ū(t) ⌘ 1 V Z V u(x, t) dV (2.25)
The projection method (Chorin, 1968) is used to calculate an intermediate velocity field ū? which is not divergence free.

ū? ūn

Dt = A n + 1 2 (D ? + D n ) + S n (2.26) ūn+1 ū? Dt = P n+1 (2.27) r • ūn+1 = 0 (2.28) P ⌘ 1 V Z V 1 r rp dV, A ⌘ 1 V Z A (uu) • ê dA, D ⌘ 1 rV Z A µ ⇣ ru + ru T ⌘ • ê dA and S ⌘ 1 V Z V F b dV. (2.29)
Note that in the integral for D, r is outside the integral as we make the approximation that r is constant over velocity control volumes. In contrast for P, r is retained inside the integral [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF]. The advection and source terms (in integral form) are treated explicitly in time, the former being discretised using a 5th order accurate WENO algorithm (Shu, 1998), described below.

Weighted Essentially Non-Oscillatory (WENO) scheme

In order to evaluate the advection term A in equation 2.29, we require the velocities at the faces of the control volume whereas these are only available at the center, as shown in figure 2.14. As seen in figure 2.15, the Weighted Essentially Non-Oscillatory(WENO) (Shu, 1998) algorithm obtains an approximation to u c+1/2 from neighboring values of u defined at the (velocity) cell centers [START_REF] Shu | Essentially non-oscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws[END_REF].

u (0) c+1/ 2 = 1 3 u c 2 7 6 u c 1 + 11 6 u c (2.30) u (1) c+1/ 2 = 1 6 u c 1 + 5 6 u c + 1 3 u c+1 (2.31) u (2) c+1/ 2 = 1 3 u c + 5 6 u c+1 1 6 u c+2 (2.32)
Once we get u k c+1/ 2 from the three sub-stencils, we can get, where w k is given by,

u c+1/ 2 = 2 Â k=0 w k u (k) c+1/ 2 (2.33) c+1/2 c c-1 c-2 c+1 c+2 S1 S0 S2
w k = wk  2 k=0 wk , wk = g k (e + b k ) 2 g 0 = 1 10 , g 1 = 3 5 , g 2 = 3 10 and e = 10 6
Expressions for b k are (Shu, 1998)

b 0 = 13 12 (u c 2 2u c 1 + u c ) 2 + 1 4 (u c 2 4u c 1 + 3u c ) 2 b 1 = 13 12 (u c 1 2u c + u c+1 ) 2 + 1 4 (u c 1 u c+1 ) 2 b 2 = 13 12 (u c 2u c+1 + u c+2 ) 2 + 1 4 (3u c 4u c+1 + u c+2 ) 2
From equation 2.29 we have

A = 1 V Z S (uu • ê)dS = 1 V  Z right (uu • ê)dz Z le f t (uu • ê)dz + Z top (uu • ê)dR Z bottom (uu • ê)dR (2.34)
The explicit expressions for the horizontal (radial with superscript R) and vertical (axial with superscript z) are,

A R r,c 1/2 = 1 2pRDRDz 2 4 Z Dq 0 Z Dz 0 u(u • ê)Rdqdz right Z Dq 0 Z Dz 0 u(u • ê)Rdqdz le f t + Z Dq 0 Z DR 0 u(u • ê)RdqdR top Z Dq 0 Z DR 0 u(u • ê)RdqdR bottom 3 5 = F r,c R r,c F r,c 1 R r,c 1 R r,c 1/2 DR + F r+1/2,c 1/2 F r 1/2,c 1/2 Dz (2.35)
Similarly,

A z r 1/2,c = F r 1,c+1/2 R r 1,c+1/2 F r 1,c 1/2 R r 1,c 1/2 R r 1,c DR + F r,c F r 1,c Dz (2.36)
where F are the fluxes calculated by WENO scheme at respective faces e.g.

F r,c+1/ 2 = u r,c+1/ 2 u r,c . For a 2D Cartesian coordinate system, equations corresponding to 2.35 and equation 2.36, become

A x r,c 1/2 = F r,c F r,c 1 Dx + F r+1/2,c 1/2 F r 1/2,c 1/2 Dy (2.37) A y r 1/2,c = F r 1,c+1/2 F r 1,c 1/2 Dx + F r,c F r 1,c Dy (2.38)

Diffusion terms

Diffusion terms are discretised using second order central difference scheme in space and Crank-Nicolson (Patankar, 1980) implicit scheme in time which removes the time step restriction due to the diffusion term. Gradient of velocity for flow in cylindrical coordinate system is given by, ru = From equation 2.29, we obtain

D = 1 rV Z V r • (µ(ru + ru T ))dV = 1 rV Z A µ(ru + ru T ) • ê dA (2.39)
Like the advection term, for a cylindrical coordinate system we obtain the discretised version of 2.39

D R r,c 1/2 = 1 r r,c 1/2  2 ✓ µ r,c (u r,c+1/2 u r,c 1/2 )R r,c µ r,c 1 (u r,c 1/2 u r,c 3/2 )R r,c 1 (DR) 2 R r,c 1/2 ◆ + ✓ µ r+1/2,c 1/2 ✓ u r+1,c 1/2 u r,c 1/2 (Dz) 2 + v r+1/2,c v r+1/2,c 1 DRDz ◆ µ r 1/2,c 1/2 ✓ u r,c 1/2 u r 1,c 1/2 (Dz) 2 + v r 1/2,c v r 1/2,c 1 DRDz ◆ 2µ r,c 1/2 u r,c 1/2 R 2 r,c 1/2 !# (2.40) and D z r 1/2,c = 1 r r 1/2,c  2 ✓ µ r,c (v r+1/2,c v r 1/2,c ) µ r 1,c (v r 1/2,c v r 3/2,c ) (Dz) 2 ◆ + ✓ µ r 1/2,c+1/2 ✓ v r 1/2,c+1 v r 1/2,c (DR) 2 + u r,c+1/2 u r 1,c+1/2 DRDz ◆ R r 1/2,c+1/2 R r 1/2,c µ r 1/2,c 1/2 ✓ v r 1/2,c v r 1/2,c 1 (DR) 2 + u r,c 1/2 u r 1,c 1/2 DRDz ◆ R r 1/2,c 1/2 R r 1/2,c ◆ (2.41)
For two-dimensional Cartesian coordinate system, equation 2.40 and equation 2.41, reduces to,

D x r,c 1/2 = 1 r r,c 1/2  2 ✓ µ r,c (u r,c+1/2 u r,c 1/2 ) µ r,c 1 (u r,c 1/2 u r,c 3/2 ) (Dx) 2 ◆ + ✓ µ r+1/2,c 1/2 ✓ u r+1,c 1/2 u r,c 1/2 (Dy) 2 + v r+1/2,c v r+1/2,c 1 DxDy ◆ µ r 1/2,c 1/2 ✓ u r,c 1/2 u r 1,c 1/2 (Dz) 2 + v r 1/2,c v r 1/2,c 1 DxDy ◆◆ (2.42)
and,

D y r 1/2,c = 1 r r 1/2,c  2 ✓ µ r,c (v r+1/2,c v r 1/2,c ) µ r 1,c (v r 1/2,c v r 3/2,c ) (Dz) 2 ◆ + ✓ µ r 1/2,c+1/2 ✓ v r 1/2,c+1 v r 1/2,c (DR) 2 + u r,c+1/2 u r 1,c+1/2 DRDz ◆ R r 1/2,c+1/2 R r 1/2,c µ r 1/2,c 1/2 ✓ v r 1/2,c v r 1/2,c 1 (DR) 2 + u r,c 1/2 u r 1,c 1/2 DRDz ◆ R r 1/2,c 1/2 R r 1/2,c ◆ (2.43)
Note that equation 2.26 is an implicit equation for u ? as the right hand side has D ? , for which we need to solve a set of linear simultaneous equations using the multigrid method.

Surface Tension

The surface tension force is a force which acts on the interface between the two fluids. The surface force per unit area can be expressed as,

F F F lg = T k n n n (2.44)
where T is the ratio of surface tension coefficient to the density and k is the curvature of the interface.

k 1 = 8 < : H 00 (R) [1+(H 0 (R)) 2 ] 3/ 2
if normal of the interface in predominantly in z orientation,

H 00 (z) [1+(H 0 (z)) 2 ] 3/ 2
if normal of the interface in predominantly in R orientation,

The second principal curvature is given by (Guo et al., 2015),

k 2 = 8 > > > > < > > > > : H 0 (R) R[1+(H 0 (R)) 2 ] 1/ 2
if normal of the interface in predominantly in z orientation,

+ 1 R[1+(H 0 (z)) 2 ] 1/ 2 if normal of the interface in predominantly in +R orientation, 1 R[1+(H 0 (z)) 2 ] 1/ 2
if normal of the interface in predominantly in R orientation,

The sum of two principal curvatures is

k = k 1 + k 2 (2.46)
For two-dimensional Cartesian coordinate system, there is only one principal curvature, hence the expression reduces to,

k = 8 > < > : H 00 (x) [1+(H 0 (x)) 2 ] 3/ 2
if normal of the interface in predominantly in y orientation,

H 00 (y) [1+(H 0 (y)) 2 ] 3/ 2
if normal of the interface in predominantly in x orientation, where derivatives are approximated using the central difference discretisation,

H 0 (R or z) = H right H le f t 2D(R or z) , H 00 (R or z) = H right 2H center + H le f t (D(R or z)) 2 (2.47)
The simulations presented in this study are sufficiently refined not to encounter the inconsistent heights (Popinet, 2009). Nevertheless, curvature for the interfaces on circle which have different orientation with respect to adjacent cells is computed as average of curvatures in R and z (or x and y for two-dimensional Cartesian) orientations, i.e.

k = 1 2 (k R + k z ) (2.48)
The time step used is constrained by advection time scale and capillary time scale and is given by,

Dt = CFL ⇥ min Dx 2U max , r (r I + r O )(Dx) 3 4pT ! (2.49)
where CFL < 1.

Pressure Poisson Equation

A pressure Poisson equation is determined subject to the incompressibility constraint on the velocity field and this pressure field is used to obtain the final divergence free velocity field i.e.

(x) = H right H le f t 2Dx , H 00 (x) = H right 2H center +H le f t (Dx) 2 return k = H 00 (x) [1+(H 0 (x)) 2 ] 3/ 2
end procedure From equation 2.29 and Gauss divergence theorem,

1 V Z A ✓ 1 r rp n+1 ◆ • ê dA = 1 Dt r • ū⇤ (2.53)
It essentially states that, any solution of velocity field should comply with mass conservation, and the field thus obtained must be divergence free. Discretising the equation 2.53 over the finite volume (Fig. 2.13).

1 R r,c+ 1 2 (DR) 2 R r,c+1 (p n+1 r,c+1 p n+1 r,c ) r r,c+1 + r r,c + R r,c (p n+1 r,c p n+1 r,c 1 ) r r,c + r r,c 1 ! + 1 (Dz) 2 p n+1 r+1,c p n+1 r,c r r+1,c + r r,c + p n+1 r,c p n+1 r 1,c r r,c + r r 1,c ! = 1 2Dt " 1 R r,c R r,c+ 1 2 u ? r,c+ 1 2 R r,c 1 2 u ? r,c 1 2 DR ! + v ? r+ 1 2 ,c v ? r 1 2 ,c Dz !# (2.54)
The equation 2.54 represents a system of linear equations for unknown p n+1 's of left hand side. The boundary conditions for pressure are Neumann for no slip, free slip and inlet boundary conditions. For outflow, the pressure on outflow boundary is taken as zero.

Note:

The pressure p here is an approximation to the physical pressure generally used in fractional step methods. A detailed discussion can be found in Rempfer (2006) on the role of the variable p in the discretised Navier-Stokes equations for incompressible flows.

Successive Over Relaxation method (SOR)

Successive over relaxation is a technique for solving a linear system of equations, and is derived by extrapolating the Gauss-Seidel method [START_REF] Tryggvason | Direct numerical simulations of gas-liquid multiphase flows[END_REF].

The extrapolation takes the form of a weighted average between the previous iterate and the computed Gauss-Seidel iterate for equation 2.54 successively which is given by,

p a+1 r,c = " 1 R r,c+ 1 2 (DR) 2 ✓ R r,c+1 r r,c+1 + r r,c + R r,c r r,c + r r,c 1 ◆ + 1 (Dz) 2 ✓ 1 r r+1,c + r r,c + 1 r r,c + r r 1,c ◆ 1 ( 1 R r,c+ 1 2 (DR) 2 R r,c+1 p a+1 r,c+1 r r,c+1 + r r,c + R r,c p a+1 r,c 1 r r,c + r r,c 1 ! + 1 (Dz) 2 p a+1 r+1,c r r+1,c + r r,c + p a+1 r 1,c r r,c + r r 1,c ! 1 2Dt " 1 R r,c R r,c+ 1 2 u ? r,c+ 1 2 R r,c 1 2 u ? r,c 1 2 DR ! + v ? r+ 1 2 ,c v ? r 1 2 ,c Dz !#) + (1 b)p a+1 r,c
(2.55)

For two-dimensional Cartesian coordinate system equation 2.54 reduces to,

p a+1 r,c =  1 (Dx) 2 ✓ 1 r r,c+1 + r r,c + 1 r r,c + r r,c 1 ◆ + 1 (Dy) 2 ✓ 1 r r+1,c + r r,c + 1 r r,c + r r 1,c ◆ 1 ( 1 (Dx) 2 p a+1 r,c+1 r r,c+1 + r r,c + p a+1 r,c 1 r r,c + r r,c 1 ! + 1 (Dy) 2 p a+1 r+1,c r r+1,c + r r,c + p a+1 r 1,c r r,c + r r 1,c ! 1 2Dt ✓ u ⇤ r,c+1 u ⇤ r,c Dx + v ⇤ r+1,c v ⇤ r,c Dy ◆ + (1 b)p a+1 r,c (2.56)
where a is the iteration step and b is the relaxation parameter which is used to accelerate the convergence. b take values from 0-2, for over-relaxation it must be greater than 1 while for stability reasons it must be below 2. A choice of b = 1.2 1.5 is usually a good compromise between stability and convergence. The advantage of the SOR is its simplicity but it converges very slowly. The solution can be accelerated using the multigrid algorithm with V-cycle (Moin, 2010). We present here a short review of the multigrid algorithm along with test cases implemented for it.

Multigrid algorithm

The pressure ), multigrid methods has been used extensively to solve the PPE. We first discuss the characteristics of classical iterative methods. We then discuss ideas from linear algebra (available in literature) which are used to accelerate convergence. We subsequently discuss the algorithm implemented in our solver for solving the PPE. The multigrid method discussed here has been implemented into an in-house developed serial code.

General Iterative Methods

We consider here one dimensional Laplace equation as an example as explained in Murthy Jayathi (2002) for the Jacobi iteration method, adapted for the Gauss-Seidel method. Consider the equation

d 2 u dx 2 = 0 ( 2.57)
with boundary conditions, u(0) = 0 and u(L) = 0, which has a trivial solution u(x) = 0. We would now investigate a numerical method to obtain this solution.

The differential equation is discretised to obtain algebraic equations,

u i+1 2u i + u i 1 D 2 = 0 (2.58) or, u i+1 2u i + u i 1 = 0 (2.59) of form, Au = b (2.60)
Iterative methods typically do not solve u = A 1 b (due to high computational costs associated with matrix inversion), but instead formulate the problem as u = Pu + Q. Iterative algorithms differ in their forms of the matrices P and Q. For Gauss-Seidel method, we split the operator A, as shown in Fig. 2.17 in lower, diagonal and upper parts i.e. A = L + D + U using which we may rewrite equation 2.60 as From 2.62 we can see that P = (D + L) 1 ( U) and Q = (D + L) 1 . For problem 2.60, the error is defined as e j = u e u j , where e j is the the error, u j is the value of u at j th iteration and the exact solution for equation 2.57 is u e = 0. The error at the jth iteration is e j = u j . We can now discuss the behavior of Gauss-Seidel method by starting with arbitrary initial guesses. To find the characteristics for different error profiles, we solve the problem with initial guess u 0 ,

(D + L)u = Uu + b (2.61) u j+1 = (D + L) 1 ( U)u j + (D + L) 1 b (2.62)
u 0 = sin ✓ kpf L ◆ (2.63)
Fig. 2.18 shows the modes over the domain for k = 1, 2 and 8. We see that for lower values of k the error is closer to zero whereas for higher wavenumbers it is not. We can see for a given number of iterations the method converges faster for higher values of wavenumber in the error. To know the dependence of wavenumber on the convergence we have to analyse the characteristics of Gauss-Seidel method.

Convergence analysis of Gauss-Seidel method

As the iterative solution u j approaches the exact solution u e , by definition the error vector e j approaches zero. We look at how the error reduces with each iteration. We have the following,

u = Pu + Q exact solution u n = Pu n 1 + Q solution at n th iteration e n = u e u n
After n recursive substitutions, e n = P n e 0 (2.64)

If we now expand e 0 in the eigen-basis of the matrix P, we obtain

e 0 = SC k v k (2.65)
where v k are eigen-vectors of iteration matrix P and C k are are the components of e 0 in the eigen basis. Operating on equation 2.65 with the matrix P n using 2.64, we obtain

e n = Sl n k C k v k (2.66)
From 2.66, it can be seen that |l k | should be less than unity for the error to approach zero in successive iterations i.e. for convergence. The large of |l k | is known as the spectral radius of the iteration matrix P. It can also be concluded from here that the smaller the spectral radius faster the convergence. We have observed in the previous section (See figure 2.18) that there is faster convergence for high wavenumber errors and we also know that convergence is directly affected by the eigenvalues of the P matrix. We now look at the relation between the wavenumber and the eigenvalues of P. The eigenvalues of Gauss-Seidel iteration matrix is given by (Briggs et al., 2000), 2.67) where N is the number of nodes in the grid. From equation 2.67 we can now understand reasons for the results in previous section i.e. the magnitude of eigenvalue is lower for high wavenumber modes and hence we found faster convergence for high wavenumbers. In most practical applications we do not know the exact solution. A common guess is to start with zero vector. Such an arbitrary initial guess typically has projections on many wavenumbers. But notice from equation 2.67, the convergence depends on the ratio k N ( See figure 2.19 ). Thus for a given wavenumber say k = 1 for which convergence is low, it can be said that this low wavenumber component will converge faster on coarser grids. This observation lies at the heart of multigrid methods.

l k = cos 2 ✓ kp 2N ◆ k = 1,2,...N (

Implementation of multigrid algorithm

In the light of above results, the strategy is to solve the problem on successively coarser grids (which kills the low wavenumber component of the error vector) and then interpolate the result on a fine grid. This is better guess for the finer grid and we can repeat the process till the finest grid when we want the solution.

We solve the given problem Ax = b on the finest grid, and Ae = r on the coarser grids. The basic steps in the algorithm are:

1. Iterate for few steps the given problem on the finest grid. This step is relaxation.

2.

Transfer the residue from finest grid to coarser grid. This step is called a restriction.

3.

Interpolate the values of error and add it previous values to make a correction. This step is prolongation.

The algorithm for multigrid is followed as in Fig. 2.20. This is V cycle algorithm (Moin, 2010), where the restriction is done till the coarsest grid and then prolongation till the finest grid.

Sample Problem

In order to understand the implementation details of the algorithm, we have solved the following problems obtained from Moin (2010).

One dimensional Poisson's equation

d 2 u dx 2 = 1 2 [sinpx + sin16px] (2.68)
having boundary conditions u(0) = u(1) = 0. We can discretise equation 2.68 as, Let the RHS of equation 2.69 be q i at i th node. Then we can write for u i as,

u i+1 2u i + u i 1 h 2 = 1 2 [sinpx i + sin16px i ] (2.69)
u n+1 i = 1 2 (u n i+1 + u n+1 i 1 h 2 q i ) (2.70)
Equation 2.70 is iteration equation for Gauss-Seidel method but for multigrid we iterate only few times. Here we will iterate only once in each step in a V-cycle.

We have the finest grid having N = 64 with total nodes N + 1 i.e. 65. We now restrict the residue to the coarser grid which has N 2 + 1 nodes. We can use average restriction as,

r 2h i = 1 4 (r h 2i 1 + 2r h 2i + r h 2i+1 ) (2.71)
Then iterate for error on the coarser grid as, where k is an arbitrary constant. The equation is discretised using a 2nd-order finite difference approximation on a Cartesian grid having N nodes in both xand y-directions which correspond to a uniform grid spacing h. The value of f on the two-dimensional Cartesian mesh can then be approximated for each node i,j in the interior of the computational domain as

e n+1 i = 1
1 h 2 ⇣ f n i 1,j + f n i,j 1 4f n+1 i,j + f n i+1,j + f n i,j+1 ⌘ = f i,j , i, j = 2, ..., N 1 (2.78)
where the subscripts i and j represent the indices of the current node in the computational domain for the x and y directions, respectively.

Performance comparison with SOR in Parallel

The two dimensional Poisson equation is solved using SOR (MPI) and the results are compared with serial multigrid. The following tests performed on the a machine having the following specifications:

1. Processor: Intel(R) Xeon(R) CPU E5-1650 v3 @ 3.50GHz 2. CPU Cores: 6 (6 more with hyperthread-Virtual cores)

3. RAM: 64 GB 4. Cache: 15360 KB From Fig. 2.23 it can be observed that time of execution has decreased with increase in number of processors and speed up is more for larger size problem where speed-up is given by time taken by n processors to the time taken by one processor for the same problem t(1)/t(n). A decrease in speed up is seen when the number of processors are more than the physical cores and the threads took more time to process. Fig. 2.24 has emphasizes the fact that multigrid on a single processor is far more efficient than a multi-core SOR algorithm, upto six cores. This work can be further extended for parallelizing the multigrid algorithm and is underway.

Test cases

We now present results from test cases. Some of these test only the Navier-Stokes implementation or the curvature implementation as described earlier. Other test cases are designed to test our implementation of the entire solver. We compare with linearised analytical results wherever available (or new results developed by us) or against experimental data taken from literature. We ⌘ r I U 2 max R 0 /T ⇡ 10 3 where U max is the maximum velocity found in the domain at t = 20 s.

Two-dimensional standing wave

The interface of two immiscible incompressible viscous fluids of density and viscosity r u , µ u , r l , µ l and surface tension T on horizontally and vertically unbounded domain (See figure 2.29a). We neglect gravity and so the base state is trivial. At time t = 0 an interfacial perturbation of the form h(x, t) = a 0 cos(kx) is introduced as shown in figure 2.29a. Due to the restoring force of surface tension, a standing wave is set up (at sufficiently small value of e ⌘ a 0 k). The interface at later time t takes the form a(t) cos(kx). The linearised initial value problem for this was solved by Prosperetti (1981) who obtained an analytical expression for a(t). We compare results obtained from our in-house developed Navier-Stokes solver with this. Symmetry boundary conditions are used on on all four sides of the computational domain. The domain is [0,1] ⇥ [0,8]. The analytical expression for ã(s) (Laplace transform of a(t)) is given by Prosperetti (1981). We numerically invert expression 2.79 using the tool provided by McClure (2016) and compare with the time signal obtained from DNS conducted using our in-house solver.

ã(s) = 1 s a(0) + su(0) w 2 0 a(0) s 2 + L(s)s + w 2 0 ! (2.79)
where,

L(s) = 4k( r l r u s + k(µ u µ l )(r u (k l l ) r l (k l u )) + k 2 (µ l µ u ) 2 (k l l )(k l u )s 1 ) (r l + r u )(r l (k l u ) + r u (k l l )) , l l,u = r k 2 + s n l,u and w 2 0 = r l r u r l + r u gk + T r l + r u k 3 r l = 1, r u = 0.01, µ l = 0.01, µ u = 0.0001, T = 1, k = 2p, g = (0, 0)

Two-dimensional Rayleigh-Taylor instability

For the case when r l < r u , the wave becomes unstable, and evolution of growth amplitude is given by the same expression 2.79 used for planar wave,

r l = 0.01, r u = 1, µ l = 0.01, µ u = 0.0001, T = 1, k = 2p, g y = 70

Axisymmetric Rayleigh-Taylor instability

For the case when r l < r u , the wave becomes unstable. It is shown in chapter 3 that the analytical expression for a(t) is the same in cylindrical axisymmetric

r l = 0.01, r u = 1, µ l = 0.0001, µ u = 0.01, s = 10, k = 1, g z = 60 2

Conclusion

An in-house solver has been developed and described in this chapter. The code uses the finite volume algorithm coupled with a piecewise-linear and LVIRA based Volume-of-fluid algorithm together with the Direction averaged curvature method for surface tension. The interface is advected using geometrical flux calculations and is benchmarked with several test cases. For all cases, the interface remained sharp and the errors have been quantified using results obtained from literature. The computational effort required to solve the pressure Poisson equation is significantly reduced by implementing the multigrid algorithm. We have benchmarked the code using test cases involving inviscid and viscous capillary waves, for each of which solutions are known in the linearised approximation. The nonlinear term in the Navier-Stokes equation, is tested with data available in literature for the dam break and rising bubble, experiments. The solver thus developed has been used in the following chapters for studying linearised viscous capillary-gravity waves in various geometries.

Axisymmetric viscous interfacial oscillations

Some the results presented in this chapter have appeared in "Axisymmetric viscous interfacial oscillations -Theory and simulations", Palas Kumar Farsoiya, Y S Mayya and Ratul Dasgupta, J. Fluid. Mech and at the 70th Annual Meeting of the American Physical Society (APS), Division of Fluid Mechanics (DFD), Denver, USA. We acknowledge collaboration with Prof. Y. S. Mayya for part of the work. The presence of the continuous spectrum may be seen from the following physical argument. The normal mode analysis of the above problem leads to an algebraic dispersion relation which admits a finite number of roots only (see figure 3.2c later on depicting the two poles). As each root also corresponds to an eigenfunction (for the vertical part of the vorticity field), this implies that the discrete normal modes led to only a finite number of eigenfunctions, to represent all possible vortical initial conditions consistent with boundary conditions. For the cylindrical pool considered here, it will be seen that these eigenfunctions decay exponentially along the depth (see figure 3.2a). Consequently, any initial condition involving a finite amount of vorticity at large depth (i.e. at distances larger than k 1 where k is the wavenumber imposed at the interface) may not be adequately represented by this discrete set of eigenfunctions only. In other words, the discrete spectrum does not provide a complete set of eigenfunctions. This was recognised by Prosperetti (1976) who solved the viscous initial value problem for planar capillary-gravity waves with a zero perturbation vorticity initial condition. Interestingly, the zero vorticity initial condition is among the simplest of initial conditions which has projections on the discrete as well as the continuous spectrum eigenfunctions. Linear analysis shows that the interface initialised as a single Fourier mode evolves subsequently as a standing wave, but whose amplitude progresses in time as a linear sum of damped exponentials (discrete spectrum) and an integral over a continuous range of exponentials (continuous spectrum). A similar situation also occurs in axisymmetric radial geometry. Here the study by Miles (1968) was the first to formulate the Cauchy-Poisson problem for free surface waves initialised as a localised disturbance. We build on the studies by Miles (1968); [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF]Prosperetti ( , 1980b) ) formulating the IVP in radial geometry for a single Fourier-Bessel mode and taking into account inertia of the fluid above the interface as well. In addition to solving the linearised IVP and comparing predictions with Direct Numerical Simulations (DNS), we also explore the non-dimensional parametric space using DNS, reporting some novel observations in the nonlinear regime. All DNS results reported here have been obtained using the in-house developed code, described in the previous chapter.

Introduction

Viscous initial value problem

As shown in figure 3.1a (dotted line), the base state comprises of two quiescent fluids of density and viscosity r L , r U and µ L ; , µ U in radially and vertically unbounded geometry. We impose an axisymmetric interfacial perturbation of the form h(r, 0) = a(0)J 0 (kr) (See figure 3.1a). and its evolution is obtained as a function of time a(t). In the linearised approximation, we solve the initial value problem (IVP) to obtain an analytical expression for the time evolution of a(t).

We employ the viscous potential flow approach involving the decomposition of the velocity field into a curl-free and a (divergence free) vortical part based on similar approaches earlier [START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Menikoff | Initial value problem for rayleigh-taylor instability of viscous fluids[END_REF][START_REF] Miles | The cauchy-poisson problem for a viscous liquid[END_REF][START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF]. i.e. the perturbation velocity v and pressure p fields are decomposed as

v = v p + v v (3.1) p = p p + p v (3.2)
The approach relies on the following identity for a vector field v p (x

) (Kundu et al., 2012), r(r • v p ) r ⇥ (r ⇥ v p ) = r 2 v p . (3.3)
For a potential flow with v p = rf p , identity 3.3 reduces to r(r 2 f p ) = r 2 (rf p ). (3.4) This implies that a potential velocity field v p and a corresponding pressure field p p obtained from the Bernoulli's equation, satisfies the Navier-Stokes equation exactly [START_REF] Joseph | Viscous potential flow[END_REF], since according to identity 3.4 the viscous term in the Navier-Stokes equation vanishes identically for a potential flow field satisfying r 2 f p = 0. However this potential flow field v p and p p , does not satisfy continuity of tangential stresses and of tangential velocities (slip) at an interface separating two fluids of different viscosities and densities. An additional velocity and pressure field, v v and p v is needed and the composite fields viz. v p + v v and p p + p v are determined subject to the constraint that each field individually satisfies the linearised Navier-Stokes as well as continuity equations while the composite field satisfies the kinematic boundary condition, continuity of tangential and normal stresses and that of tangential velocities. Due to linearity, the composite solution will automatically satisfy the linearised Navier-Stokes equation as well. We determine v p and p p from solving the Laplace and linearised Bernoulli's equation respectively.

Potential Flow

The potential solution satisfies the axisymmetric Laplace equation. From variable separation, the radial part of the solution is given by Bessel function of the first kind (boundedness at r = 0, eliminates the Bessel function of the second kind). We set

f U p (r, z, t) = F(z)J 0 (kr) . a(t), f L p (r, z, t) = G(z)J 0 (kr) . a(t), (3.5)
with h(r, t) = a(t)J 0 (kr). We set the time dependence of f p equal to

∂ 2 f p ∂r 2 + 1 r ∂f p ∂r + ∂ 2 f p ∂z 2 = 0. (3.6)
Substituting equation 3.5 in this and using the equation for J 0 (r) (r ⌘ kr),

d 2 J 0 (r) dr 2 + 1 r dJ 0 dr + J 0 (r) = 0, (3.7)
we obtain

d 2 F dz 2 + k 2 F = 0, d 2 G dz 2 + k 2 F = 0. (3.8) Using f U (r, •, t) ! 0 and f L (r, •, t) ! 0,
and the linearised kinematic boundary condition

∂f U p ∂z z=0 = ∂f L p ∂z z=0 = ∂h ∂t (3.9)
we obtain

f U p (r, z, t) = k 1 exp[ kz]J 0 (kr)
.

a(t), f L p (r, z, t) = k 1 exp[kz]J 0 (kr) . a(t)(3.10)
The potential part of pressure p U p and p L p is given by the linearised Bernoulli's equation,

p L p (r, z, t) = r L ∂f L p ∂t r L gz p U p (r, z, t) = r U ∂f U p ∂t r U gz. (3.11)
If the flow is purely irrotational, then continuity of pressure at the linearised interface, leads to the following simple harmonic oscillator equation for a(t),

..

a(t) + ✓ r L r U r L + r U ◆ gk a(t) = 0 (3.12)
Note that the frequency of the oscillator in 3.12 is the dispersion relation for surface-gravity waves on an infinite pool. In subsequent analysis, we include viscous effects for which instead of imposing continuity of pressure, we impose a condition on viscous normal stresses (accounting for jump due to surface tension) at the interface. This will lead to a modified viscous equation for a(t), derived in the next section.

Viscous Flow

The viscous part of the flow satisfies the linearised Navier-Stokes equation in both the fluids

∂v v U ∂t = 1 r U rp U v + n U r 2 v v U , ∂v L v ∂t = 1 r L rp L v + n L r 2 v v L . (3.13)
Note that gravity has already been included as a body force in the potential part of the calculation and hence is excluded in the vortical calculation. We solve the viscous part of the flow in stream-function vorticity formulation. The curl of equations 3.13 gives us the vorticity equation

∂w U ∂t = n U r 2 w U , ∂w L ∂t = n U r 2 w L . (3.14)
The Stokes stream function y (Miles, 1968) is defined as

u U v ⌘ ∂y U v ∂z , v U v ⌘ 1 r ∂(ry U v ) ∂r u L v ⌘ ∂y L v ∂z , v L v ⌘ 1 r ∂(ry L v ) ∂r (3.15)
The azimuthal vorticity (w) q ⌘ w q , is given by

w U q ⌘ ∂u U v ∂z ∂v U v ∂r , w L q ⌘ ∂u L v ∂z ∂v L v ∂r .
(3.16)

Combining equation 3.15 and 3.16, the relation between vorticity and streamfunction is written as

∂ 2 y U v ∂r 2 y U v r 2 + 1 r ∂y U v ∂r + ∂ 2 y U v ∂z 2 = w U q , ∂ 2 y L v ∂r 2 y L v r 2 + 1 r ∂y L v ∂r + ∂ 2 y L v ∂z 2 = w L q . (3.17)
An equation for w q can be obtained from equation 3.14 using the expression for the Laplacian of a vector in cylindrical axisymmetric coordinates (Kundu and Cohen, 2002)),

∂w U q ∂t = n U ∂ 2 w U q ∂r 2 w U q r 2 + 1 r ∂w U q ∂r + ∂ 2 w U q ∂z 2 ! , ∂w L q ∂t = n L ∂ 2 w L q ∂r 2 w L q r 2 + 1 r ∂w L q ∂r + ∂ 2 w L q ∂z 2 ! (3.18)
The solution of equations 3.17 and 3.18 subject to boundary and initial conditions, determine the flow field in both fluids. We set the radial part of all quantities as Bessel function of the first kind J 1 (kr) viz.

w U q (r, z, t) = W U (z, t)J 1 (kr), y U v (r, z, t) = Y U (z, t)J 1 (kr), (3.19) w L q (r, z, t) = W L (z, t)J 1 (kr), y L v (r, z, t) = Y L (z, t)J 1 (kr). (3.20)
Substitution of equations 3.19 and 3.20 in equation 3.17, gives us an equation relating W to Y,

∂ 2 Y U ∂z 2 k 2 Y U = W U (z, t), ∂ 2 Y L ∂z 2 k 2 Y L = W L (z, t). (3.21)
Note that in deriving equation 3.21, we have used the equation governing J 1 (r) i.e.

d 2 J 1 (r) dr 2 + 1 r dJ 1 dr + ✓ 1 1 r2 ◆ J 1 (r) = 0. (3.22)
An equation involving W alone is obtained by substituting expressions for vorticity from equations 3.19 and 3.20 into equations 3.18 and then using equation 3.22. .23) equations 3.21 and 3.23 are the central equations. These turn out to be identical to those derived by Prosperetti (1981) for a planar geometry. We follow the approach of Prosperetti (1981) for solving these equations. These steps involve additional manipulations with Bessel functions not necessary in the planar case.

∂W U ∂t = n ✓ ∂ 2 W U ∂z 2 k 2 W U ◆ , ∂W L ∂t = n ✓ ∂ 2 W L ∂z 2 k 2 W L ◆ . ( 3

Viscous pressure

In order to obtain an equation for a(t), we will need an expression for the viscous part of pressure p v in both the fluids. This is obtained by integrating the vertical momentum equation. We demonstrate the algebra for the lower fluid. The (linearised) vertical momentum equation is (Kundu and Cohen, 2002)

1 r L ∂p L v ∂z = n L  1 r ∂ ∂r ✓ r ∂v L v ∂r ◆ + ∂ 2 v L v ∂z 2 ∂v L v ∂t (3.24)
Note that gravity is already included in the potential flow and hence is not present in the viscous part of the calculation. With

v L v = 1 r ∂(ry L v ) ∂r and equation 3.20, v L v = kY L (dJ 1 (r)/dr + J 1 /r).
Substituting this in equation 3.24, using the equation governing J 0 and the identity J 0

1 (x) + J 1 (x)/x = J 0 (x), prime indicating differentiation, we obtain 1 r L ∂p L v ∂z = kJ 0 (kr)  n L ✓ k 2 Y L v ∂ 2 Y L ∂z 2 ◆ + ∂Y L ∂t (3.25)
Equation 3.25 can be further simplified by using equation 3.21 on the right hand side to rewrite it as

1 r L ∂p L v ∂z = kJ 0 (kr)  n L W L + ∂Y L ∂t (3.26)
A similar equation can be written for the viscous pressure in the upper fluid.

Boundary/Initial conditions

In this section, we provide a derivation of the kinematic boundary condition, continuity of shear stress and tangential velocities at the interface and initial conditions.

Initial conditions

There is no motion in the fluid at t = 0 (we also assume . a(0) = 0 later). Hence,

W U (z, 0) = W L (z, 0) = Y U (z, 0) = Y L (z, 0) = 0. (3.27)

Decay at ±•

All quantities decay to zero at ±• and hence

Y U (•, t) = W U (•, 0) = 0, Y L ( •, t) = W L ( •, 0) = 0 (3.28)

Kinematic Boundary Condition

The (linearised) kinematic boundary condition is

v U (r, 0, t) = v L (r, 0, t) = ∂h ∂t . (3.29)
Writing the velocity as a sum of potential and viscous parts equation 3.29 becomes,

∂f U ∂z z=0 + v L v (r, 0, t) = ∂f L ∂z z=0 + v L v (r, 0, t) = ∂h ∂t . (3.30)
Subtracting the kinematic boundary condition for potential flow viz. equation 3.9 from 3.30 we find that

v U v (r, 0, t) = v L v (r, 0, t) = 0 (3.31)
Thus,

Y U (0, t) = Y L (0, t) = 0 (3.32)
In deriving equation 3.32 we have used the identity J 0

1 (x) + J 1 (x)/x = J 0 (x). Also note that equation 3.31 implies ∂v U v ∂r z=0 = ∂v L v ∂r z=0 . (3.33)

Continuity of shear stress

The shear stress t rz is continuous at the linearised interface implying, .34) where

µ U ✓ ∂v U ∂r + ∂u U ∂z ◆ z=0 = µ L ✓ ∂v L ∂r + ∂u L ∂z ◆ z=0 . ( 3
u U = u U p + u U v , v U = v U p + v U v and
µ U ∂v U p ∂r + ∂u U p ∂z + ∂u U v ∂z ! z=0 = µ L ∂v L p ∂r + ∂u L p ∂z + ∂u L v ∂z ! z=0 . (3.35) equation 3.
35 can be used to relate the vorticity on either side of the interface at z = 0. As the interface acts as a source of vorticity, the vorticity field is discontinuous at the linearised interface z = 0. The potential part of the flow satisfies .36) everywhere in the flow. This can be used in equation 3.35 which simplifies to,

∂u U p ∂z = ∂v U p ∂r , ∂u L p ∂z = ∂v L p ∂r , ( 3
µ U 2 ∂v U p ∂r + ∂u U v ∂z ! z=0 = µ L 2 ∂v L p ∂r + ∂u L v ∂z ! z=0 (3.37)
We now evaluate the following expression utilising equation 3.37.

⇣ µ U w U q µ L w L q ⌘ z=0 = µ U ✓ ∂u U v ∂z ∂v U v ∂r ◆ z=0 µ L ✓ ∂u L v ∂z ∂v L v ∂r ◆ z=0 (3.38)
Using equation 3.33, equation 3.38 can be simplified to obtain

⇣ µ U W U (0, t) µ L W L (0, t) ⌘ J 1 (kr) = ✓ µ U ∂u U v ∂z µ L ∂u L v ∂z ◆ z=0 (3.39)
Using equation 3.37 to simplify the right hand side of equation 3.39, we obtain

⇣ µ U W U (0, t) µ L W L (0, t) ⌘ J 1 (kr) = 2 µ L ∂v L p ∂r µ U ∂v U p ∂r ! z=0 = 2 µ L ∂ 2 f L p ∂r∂z µ U ∂ 2 f U p ∂r∂z ! z=0
(3.40) Equation 3.40 can be further simplified to,

µ U W U (0, t) µ L W L (0, t) = 2k(µ L µ U ) . a(t) (3.41)
where we have used the relation

dJ 0 (x) dx = J 1 (x)
to eliminate the radial dependence from both sides.

Continuity of tangential velocities

It is clear that the tangential velocities are not continuous at z = 0 for potential flow. Imposing the continuity of tangential velocities at the linearised interface implies,

∂f U p ∂r + ∂y U v ∂z ! z=0 = ∂f L p ∂r + ∂y L v ∂z ! z=0 (3.42)
Using expressions derived earlier, this can be written as .

a(t) + ∂Y U ∂z (0, t) = . a(t) + ∂Y L ∂z (0, t) (3.43)
where we have used

dJ 0 (x) dx = J 1 (x)
. Equation 3.43 can be rearranged to obtain

∂Y L ∂z (0, t) ∂Y U ∂z (0, t) = 2 .
a(t).

(3.44)

Laplace Transforms

Equations in the Laplace domain

Equations 3.21 and 3.23 are easily solved in the Laplace domain. We define the Laplace transform of W(z, t) and Y(z, t) as

L (W(z, t)) ⌘ W(z, s) = Z • 0 W(z, t) exp[ st]dt, Ỹ(z, s) ⌘ Z • 0 Y(z, t) exp[ st]dt.
Note that all quantities with a tilde on top are in the Laplace domain. The Laplace transform of equations 3.21 leads to,

∂ 2 ỸU ∂z 2 k 2 ỸU = WU (z, s), ∂ 2 ỸL ∂z 2 k 2 ỸL = WL (z, s). (3.45)
Laplace transforming equations 3.23 and using initial conditions in equation 3.27, we obtain

∂ 2 WU (z, s) ∂z 2 ⇣ k 2 + s n U ⌘ WU (z, s) = 0 ∂ 2 WL (z, s) ∂z 2 ⇣ k 2 + s n L ⌘ WL (z, s) = 0 (3.46)
Laplace transform of the viscous contribution to pressure in equation 3.26 leads to

1 r L ∂ pL v ∂z = kJ 0 (kr) h n L WL + s ỸL i (3.47)
where Y L (z, 0) = 0 is taken as initial condition. Similarly we can find an expression for pressure for the upper fluid,

1 r U ∂ pU v ∂z = kJ 0 (kr) h n U WU + s ỸU i (3.48)
Laplace transforming equations 3.28, 3.32, 3.41 and 3.44, we obtain

Ỹ(•, s) = WU (•, 0) = 0, ỸL ( •, t) = WL ( •, 0) = 0 (3.49) ỸU (0, s) = ỸL (0, s) = 0. (3.50) µ U WU (0, s) µ L WL (0, s) = 2k(µ L µ U )L ( . a(t)) , (3.51) ∂ ỸL ∂z (0, s) ∂ ỸU ∂z (0, s) = 2L ( . a(t)) , (3.52)

Solution in Laplace Domain

Our task is to first solve equations 3.45, 3.46 using equations 3.49 and 3.51. The solution of equation

3.46 is WU (z, s) = ÃU (s) exp [ zl U ] , WL (z, s) = ÃL (s) exp [zl L ] . (3.53)
Here A L (s) and A U (s) are constants of integration and following Prosperetti (1981), we define l U ⌘ q k 2 + s n U and so on. Equations 3.45 are solved using 3.49 and 3.50 to obtain, ỸU (z, s) = ÃU (s)

l 2 U k 2 (exp[ zl U ] exp[ kz]) (3.54) ỸL (z, s) = ÃL (s) l 2 L k 2 (exp[zl L ] exp[kz]) (3.55)
Equations 3.51 and 3.52 can now be written as

µ U ÃU (s) µ L ÃL (s) = 2k(µ L µ U )L[ . a(t)], (3.56) ÃL (s) k + l L + ÃU (s) k + l U = 2L[ . a(t)] (3.57)
Equations 3.56 and 3.57 can be solved for ÃL (s) and ÃU (s) to obtain,

ÃU (s) = 2L( . a(t))(k + l U )(µ U k + µ L l L ) µ L (k + l L ) + µ U (k + l U ) , (3.58) and ÃL (s) = 2L( . a(t))(k + l L )(µ L k + µ U l U ) µ L (k + l L ) + µ U (k + l U ) . (3.59)
The expression for the viscous part of pressure can also be simplified further. Equation 3.47 can be written as

1 r L ∂ pL v ∂z = kJ 0 (kr) h n L ÃL (s) exp [zl L ] + s ỸL i (3.60)
Equation 3.55 can be rewritten as

s ỸL (z, s) = ÃL (s)n L (exp[zl L ] exp[kz]) (3.61)
Using equation 3.61 in 3.60 we obtain

1 r L ∂ pL v ∂z = kJ 0 (kr) ÃL (s)n L exp[kz] (3.62)
which can be integrated from

z 0 = • to z 0 = z (z 0 being dummy variable) to obtain pL v (r, z, s) = µ L ÃL (s)J 0 (kr) exp[kz]. (3.63)
Using a similar procedure and integrating from z 0 = z to z 0 = • for the upper fluid,

pU v (r, z, s) = µ U ÃU (s)J 0 (kr) exp[ kz]. (3.64)
The potential part of the flow in Laplace domain is written below.

fU p (r, z, s) = k 1 exp[ kz]J 0 (kr)L ( . a(t)) , fL p (r, z, s) = k 1 exp[kz]J 0 (kr)L ( . a(t)) (3.65) pL p (r, z, s) = r L ⇣ s fL p f L p (r, z, 0) ⌘ r L gz (3.66) p U p (r, z, s) = r U ⇣ s fU p f U p (r, z, 0) ⌘ r U gz. (3.67)

Final assembly

We can now obtain an equation for ã(s) by taking into account the jump in normal stresses across the interface due to surface tension. This is (cf. Bush (2013)), (s zz is the zz component of stress tensor)

s U zz (r, 0, t) s L zz (r, 0, t) = T(r • n). (3.68)
where n is the local unit normal to an axisymmetric interface z = h(r, t) = a(t)J 0 (kr). The local radius of curvature is given by (Bush, 2013),

r • n = rh r 1 + h 2 r r 2 h rr r 2 (1 + h 2 r ) 3/2 .
(3.69) Equation 3.69 when linearised for small-amplitude oscillations becomes,

r • n ⇡ h rr + 1 r h r = a(t) d 2 J 0 dr 2 + 1 r dJ 0 dr = k 2 a(t)J 0 (kr). (3.70)
Equations 3.68 and 3.70 can be combined to obtain,

p L (r, 0, t) p U (r, 0, t) + 2µ U ∂v U ∂z (r,0,t) 2µ L ∂v L ∂z (r,0,t) = Tk 2 a(t)J 0 (kr) (3.71)
Note that we have used the Newtonian constitutive relation s zz = p + 2µ ∂v ∂z . Using the decomposition p = p p + p v and v = v p + v v for upper and lower fluids, equation 3.71 can be rewritten after Laplace transformation as,

pL p (r, 0, s) + pL v (r, 0, s) pU p (r, 0, s) pU v (r, 0, s) + 2µ U ∂ ṽU p ∂z (r,0,s) + 2µ U ∂ ṽU v ∂z (r,0,s) 2µ L ∂ ṽL p ∂z (r,0,s) 2µ L ∂ ṽL v ∂z (r,0,s) = Tk 2 L(a(t))J 0 (kr) (3.72)
From equations 3.66 and 3.67 we have

pL p (r, 0, s) = r L ⇣ s fL p f L l (r, 0, 0) ⌘ r L gJ 0 (kr)L(a(t)) = r L h sk 1 J 0 (kr)L( . a(t)) + k 1 J 0 (kr) . a(0) i r L gJ 0 (kr)L(a(t)) (3.73)
Note that the gravity term has been evaluated at z = h to stay consistent with linear approximation, Similarly,

pU p (r, 0, s) = r U ⇣ s fU p f U l (r, 0, 0) ⌘ r U gJ 0 (kr)L(a(t)) = r U h sk 1 J 0 (kr)L( . a(t)) k 1 J 0 (kr) . a(0) i r U gJ 0 (kr)L(a(t)) (3.74)
From equation 3.63 and 3.64, we obtain

pL v (r, 0, t) pU v (r, 0, t) = ⇣ µ L ÃL + µ U ÃU ⌘ J 0 (kr) (3.75) Also, 2µ U ∂ ṽU p ∂z z=0 2µ L ∂ ṽL p ∂z z=0 = 2k(µ U + µ L )J 0 (kr)L( . a(t)) (3.76)
It can be shown that

ṽL v = k ỸL (z, s)J 0 (kr) (3.77)
from which we obtain,

2µ U ∂ ṽU v ∂z z=0 2µ L ∂ ṽL v ∂z z=0 = 2kµ U ∂ ỸU ∂z z=0 J 0 (kr) + 2kµ L ∂ ỸL ∂z z=0 J 0 (kr) = 2kJ 0 (kr) ✓ µ U l U + k ÃU (s) + µ L l L + k ÃL (s) ◆ (3.78)
In the final assembly, we substitute expressions from equations 3.73, 3.74, 3.75, 3.76 and 3.78 into equation 3.72 and divide throughout by J 0 (kr) to obtain,

(r L + r U )k 1 s [s ã(s) a(0)] + ⇣ r L + r U ⌘ k 1 . a(0) ⇣ r L r U ⌘ g ã(s) ⇣ µ L ÃL + µ U ÃU ⌘ 2k ⇣ µ U + µ L ⌘ [s ã(s) a(0)] +2k ✓ µ U ÃU l U + k + µ L ÃL l L + k ◆ = Tk 2 ã(s) (3.79)
which can be rewritten as,

(r L + r U )s 2 ã (r L + r U )sa(0) + 2k 2 ⇣ µ U + µ L ⌘ s ã ⇣ r L + r U ⌘ . a(0) 2k 2 ⇣ µ U + µ L ⌘ a(0) + h (r L r U )gk + Tk 3 i ã + k ⇣ µ L ÃL + µ U ÃU ⌘ 2k 2 ✓ µ U ÃU l U + k + µ L ÃL l L + k ◆ = 0 (3.80)
We define x(s) and z(s) from equations 3.58 and 3.59 as

ÃU (s) = z(s)(s ã a(0)) ÃL (s) = x(s)(s ã a(0)) (3.81) Thus z(s) ⌘ 2(k + l U )(µ U k + µ L l L ) µ L (k + l L ) + µ U (k + l U ) , x(s) ⌘ 2(k + l L )(µ L k + µ U l U ) µ L (k + l L ) + µ U (k + l U ) (3.82)
For this study, we will consider . a(0) = 0. Hence using equations 3.81, equation 3.80 can be rewritten as

h⇣ r U + r L ⌘ s 2 + 2k 2 ⇣ µ U + µ L ⌘ s + ⇣ r L r U ⌘ gk + Tk 3 + k ⇣ µ L x + µ U z ⌘ s 2k 2 ✓ µ U z l U + k + µ L x l L + k ◆ s ã(s) = h (r L + r U )s + 2k 2 ⇣ µ U + µ L ⌘ +k ⇣ µ L x + µ U z ⌘ 2k 2 ✓ µ U z l U + k + µ L x l L + k ◆ a(0) (3.83)
which can be expressed as

ã(s) = h s + 1 r U +r L n 2k 2 (µ U + µ L ) + k(µ U z + µ L x) 2k 2 ⇣ µ U z l U +k + µ L x l L +k ⌘oi s 2 + 1 r U +r L n 2k 2 (µ U + µ L ) + k(µ U z + µ L x) 2k 2 ⇣ µ U z l U +k + µ L x l L +k ⌘o s + w 2 0 a (0)
(3.84) Equation 3.84 may be compactly written as (in the similar form as [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF] in Cartesian geometry)

ã(s) = 1 s 1 w 2 0 s 2 + L(s)s + w 2 0 ! a(0) (3.85)
in a form analogous to that obtained in Prosperetti (1981) in Cartesian geometry, where

w 2 0 = gk ✓ r L r U r L + r U ◆ + Tk 3 r L + r U (3.86)
and,

L(s) = 1 r U + r L ⇢ 2k 2 (µ U + µ L ) + k(µ U z + µ L x) 2k 2 ✓ µ U z l U + k + µ L x l L + k ◆ (3.87)

Analytical expressions for n

L = n U = n
Analytical expressions for the inversion of 3.85 into the time domain are possible Chandrasekhar, 1961b;[START_REF] Menikoff | Initial value problem for rayleigh-taylor instability of viscous fluids[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF] in the limit of equal kinematic viscosities of both fluids. In this limit l L = l U = l, and we obtain

(
z(s) = 2(r U k + r L l) r L + r U , x(s) = 2(r L k + r U l) r L + r U .
(3.88) Equation 3.84, can now be written as

ã(s) = ⇢ s + 2nk 2 r L +r U  r L + r U + 1 2k r L x + r U z ⇣ r L x+r U z l+k ⌘ a(0) s 2 + 2nk 2 r L +r U  (r L + r U ) + 1 2k (r L x + r U z) ⇣ r L x+r U z l+k ⌘ s + w 2 0 (3.89)
which may be rewritten as

ã(s) a(0) = s + L(s) s 2 + L(s)s + w 2 0 = L(s) Q(s) (3.90)
It is to be noted here that within the linearised framework, the expressions for ã(s) as presented in equation 3.85, turn out to be the same for rectangular and a cylindrical pool. This has been shown explicitly in Farsoiya et al. (2017). For inversion of the above expression in time domain, we use a Bromwich contour integral as shown in figure 3.2c. Care has to be taken while doing the contour integration, as L(s) in equation 3.90 has a square root branch point at s = nk 2 . On the complex s plane with p  arg(s)  p, we draw a "keyhole" contour in figure 3.2c with a branch cut running from s = nk 2 to • on the negative real axis. The integrals over the contours 2, 4 and 6 do not contribute in the limit R ! • and d ! 0. Using Cauchy residue theorem, it may be shown that

a(t) a(0) = 1 2pI Z c+I• c I• L(s) Q(s) exp[st]ds = Â m L(s m ) Q 0 (s m ) exp(s m t) n 2pI Z • 0 ⇢✓ L(s) Q(s) ◆ s=s + ✓ L(s) Q(s) ◆ s=s exp[ (x + k 2 )nt]dx (3.91)
where s i are the poles of the expression 3.90 and s ± = n(x exp[±Iq] k 2 ). For k > 0, any real value of c > 0 may be used and the prime in Q 0 denotes the derivative. Note that the ratio L(s)/Q(s) satisfies the relation

L(s) Q(s) = ✓ L(s) Q(s) ◆ .
Consequently, the expression inside the braces in the integral of equation 3.91 is a purely imaginary quantity implying 3.91 may be written as

a(t) a(0) = Â m L(s m ) Q 0 (s m ) exp(s m t) n exp[ nk 2 t] p Z • 0 = ⇢✓ L(s) Q(s) ◆ q=p exp[ nxt]dx (3.92)
Similarly, from equation 3.53 (See Appendix I.1 for the algebra),

W L (t, z) a(0)w 2 0 = Â m x(s m ) Q 0 (s m ) exp[l m z] exp(s m t) + 2n exp[ nk 2 t] p Z • 0 (b(x) cos[ p xz] + a(x) sin[ p xz]) exp[ nxt]dx (3.93)
where,

a(x) ⌘ r L kA r U p xB (r L + r U )(A 2 + B 2 ) , b(x) ⌘ r L kB + r U p xA (r L + r U )(A 2 + B 2 ) A ⌘ n 2 ⇣ x 2 k 4 ⌘ 2k 2 n 2 x k 2 r 2 L + r 2 U 4xr L r U (r L + r U ) 2 + w 2 0 , B ⌘ 4k 2 n 2 p x k 2 r 2 L + r 2 U + x k 2 r L r U k (r L + r U ) 2 (3.94)
and,

W U (t, z) a(0)w 2 0 = Â m z(s m ) Q 0 (s m ) exp[l m z] exp(s m t) + 2n exp[ nk 2 t] p Z • 0 (b(x) cos[ p xz] + a(x) sin[ p xz]) exp[ nxt]dx (3.95)
where,

a(x) ⌘ r U kA r L p xB (r L + r U )(A 2 + B 2 ) , b(x) ⌘ r U kB + r L p xA (r L + r U )(A 2 + B 2 ) A ⌘ n 2 ⇣ x 2 k 4 ⌘ 2k 2 n 2 x k 2 r 2 L + r 2 U 4xr L r U (r L + r U ) 2 + w 2 0 , B ⌘ 4k 2 n 2 p x k 2 r 2 L + r 2 U + x k 2 r L r U k (r L + r U ) 2 (3.96)
Note an important feature seen from equation 3.93 that at t = 0 it may be written as equation 3.97. The significance of this is as follows. For a given set of fluid parameters and wavenumber k, the vorticity eigenfunctions given by the real and imaginary parts of exp [l m z] (see fig.

W L (0, z) a(0)w 2 0 = Â m x(s m ) Q 0 (s m ) exp[l m z] + 2n p Z • 0 ✓ b(x) cos[ p xz] + a(x) sin[ p xz] ◆ dx (3.97)

Normal mode approximation

For normal mode analysis, we substitute

h(r, t) = a(0)J 0 (kr) 1 2 [exp[st] + c.c], f(r, z, t) = F(z)J 0 (kr) 1 2 [exp[st] + c.c] w q (r, z, t) = J 1 (kr) 1 2 [W(z) exp[st] + c.c], y v (r, z, t) = J 1 (kr) 1 2 [Y(z) exp[st] + c.c].
(3.98)

A dispersion relation may be found using standard techniques. However, as the IVP has been solved earlier we can obtain the dispersion relation directly from the denominator of equation 3.85. This is,

s 2 + L(s)s + w 2 0 = 0 (3.99)
where L is defined in equation

Direct Numerical Simulations

We solve the in-house axisymmetric incompressible Navier-Stokes equations as described in Chapter 2. Refer to figure 3.3a, the interface is initiated as a(0)J 0 (kr).

The domain extends upto the third extrema of J 0 along the radial direction although the results are independent of this. The vertical domain size was arrived at by checking for domain size independence in the vertical direction. Since our theoretical results are for axially and radially unbounded domain, we impose symmetry boundary conditions on all four sides of the computational domain in figure 3.3a. This implies zero normal velocities and Neumann conditions on tangential velocity, pressure and volume fractions on boundaries. In our input script file, we approximate J 0 (kr) by numerical integration using the following definition,

J 0 (kr) = 1 p Z p 0 cos (kr sin (l)) dl.
With a 0 ⌘ a(0), the non-dimensional numbers which govern the problem are non-dimensional amplitude e ⌘ a 0 k, density ratio r r ⌘ r L /r U , viscosity ratio

µ r ⌘ µ L /µ U , Laplace number La ⌘ Tr U /((µ U ) 2 k) and Galilei number Ga ⌘ g(r U ) 2 /((µ U ) 2 k 3 )
. The parametric space is large and the simulations carried out for the present study are listed in Table 3.1. The choice for the range of values for La and Ga arise from the fact that for air-water combination, these numbers are very large e.g. for air as the upper fluid and water as lower fluid in earth's gravitational field in CGS units, T = 72, r U = 1.5 ⇥ 10 3 , r L = 1, µ U = 1.8 ⇥ 10 4 , k = 1, g = 981 leads to Ga = 6.81 ⇥ 10 4 and La = 3.33 ⇥ 10 6 . While most of the simulations are done with a domain size of [0, 7.0155] ⇥ [0, 7.0155], for others (Table 3.1) larger domains were necessary (see section 3.5).

Results and Discussion

A comparison of the normal mode approximation (equation

⌘ gk ⇣ r L r U r L +r U ⌘ + Tk 3 r U +r L .
In the case of simulations, the location of the interface is tracked at the axis of symmetry i.e. r = 0 (although any other radial location also shows the same trend). Figure 3.4a and 3.4b correspond to underdamped oscillations while 3.4c and 3.4d are for overdamped oscillations. The excellent agreement between the solution to the IVP and the Navier-Stokes simulation is apparent. It is also seen that while normal mode is a good approximation to the IVP solution (equation 3.85) for low values of µ r , the difference between the two increases at higher values of viscosity ratio especially at intermediate times. (1982) for analogous observations for large amplitude droplet oscillations). Figure 3.5b shows the effect of increasing e on amplitude a(t) as a function of time. For e = 0.01 and 0.1, the IVP solution given by equation 3.85 is an excellent approximation. For e larger than 0.5, there are marked deviations due to non-linear effects. For increasing values of e, the interface at the axis of symmetry starts rising to a height greater than a(0). We call this phenomena 'jetting'. Note that for pure capillary waves in figure 3 The effect of increasing e is visualised as a collage in figure 3.6 and 3.7. Both these figures are to be read along columns. Each column shows the interface evolving in time for a particular value of e, all other non-dimensional parameters being the same. The background colour contours represent the vorticity field.

La Ga L 1 L 2 L 3 L 4 a 0 k r L r U µ L µ U Tr U (µ U ) 2 k g(r U ) 2 (µ U ) 2 k 3 7.
It is clearly seen from figure 3.6 that most of the vorticity is limited to a thin region around the interface, with the maximum value of vorticity being an increasing function of e. It is also seen that the peak value of vorticity (magnitude) seems to concentrate in the lighter fluid (see Govindarajan (2010) for related analysis of this). Further the magnitude of the vorticity fields is an increasing function of e as seen from the increasing scale of the color-bars. Impending droplet pinch-off is clearly visible in the last image of all three columns in figure 3.7. For e larger than 0.5, non-linear effects are clearly visible e.g. the instantaneous shape of the interface at the row corresponding to t = T 0 /2 (T 0 ⌘ 2p/w 0 ) in all the columns in figure 3.7 is very different from J 0 (kr), the assumption h(r, t) = a(t)J 0 (kr). However at the end of one oscillation (for e  0.5), the interface returns to a shape which is approximately J 0 (kr) again, 1 Cases 14 and 15 are also used for normal mode comparison with DNS. increased (e = 0.5) , systematic deviations are seen from linear predictions. While the shapes of the vorticity profiles at different radial locations appear broadly the same, a collapse of these profiles by scaling with J 1 (kr) no longer appears possible, hinting at the appearance of higher modes. It is also seen that there is a vorticity jump at the interface. The inset in the first column in figure 3.6 shows this jump as we cross the interface. This is consistent with the linearised predictions from equation 3.51 reflecting the fact that the interface acts as a source of vorticity (Brøns et al., 2014). The vorticity jump at the interface can however be minimised by choosing the viscosity ratio µ r = 1. To see this, we rewrite equation 3.51 as,

W U (0, t) = µ r W L (0, t) 2k(µ r 1)
.

a(t) (3.101)
It is clear that for µ r = 1, the vorticity jump between the upper and lower fluid at the interface goes to zero i.e W L (0, t) = W U (0, t). This prediction is borne out by simulations. In the insets of Figs. 3.9a, 3.9b and 3.9c as we make µ r = 1 (figure 3.9b), the difference between the colour of vorticity contours reduces as we cross the interface, although it does not become zero. Note that equation 3.101 predicts that for µ r 6 = 1, the vorticity jump across the interface depends both on viscosity ratio as well as the density ratio, though . a(t). It should be remembered that equation 3.101 is a linearised prediction applicable to the vorticity values immediately above and below the interface only. figure 3.10 shows the effect of variation of parameters on underdamped oscillations. figure 3.10a shows the influence of Galilei number. At larger values of Ga, one can collapse all the curves by scaling t with w 0 . However this collapse fails as one reduces the value of Ga. Note however that for surface tension driven, pure capillary oscillations a collapse is possible using the time scale proposed in Denner (2016). Increasing the viscosity ratio in figure 3.10b increases the damping but the time period remains relatively unaffected. Increasing the density ratio in figure 3.10c has an effect both on the time period (as one can estimate from the inviscid dispersion relation) as well as on the damping. Note that for r r = 1, the effect of gravity drops out and these become underdamped pure capillary oscillations. For r r = 1000 (µ r = 100), the inviscid-irrotational approximation of an undamped harmonic oscillator with frequency w 0 is an excellent approximation for the first few oscillations. The effect of change of Laplace number is shown in figure 3.10d where we see some damping but hardly any change in the time period. In order to understand the modal structure of the interface and study the emergence of new modes due to nonlinearity, we use Hankel transform of the interface defined by the following 3.1) for evaluating the Hankel transform. In all cases, we checked and found that the interface evolution was independent of the domain size to working precision. The interface evolution is shown in figure 3.11 while the corresponding Hankel transform of the interface is shown in figure 3.12. The columns in both figures represent e = 0.5, 1.5 and 3.5 from left to right with each figure showing the interface (figure 3.11) and Hankel transform (figure 3.12) at various instants of time. It is seen that the second mode k = 2 starts appearing with a magnitude comparable to k = 1 for e = 1.5. For e = 3.5, the second, third and fourth modes (k = 2, 3 and 4) also appear. It is interesting to note that for e greater than 0.5, at the end of one inviscid time-period, higher modes survive along with the initial mode k = 1. Note that in figure 3.11 last column, the interface becomes a multivalued function of radius near r = 0. The Hankel transform for this case was obtained by excluding this small region of multivaluedness near r = 0. The bubble formation can be seen when the cavity collapse process starts. A vortex ring travelling downwards can be seen in figure We study here the evolution of the interface for a localised initial perturbation at the interface. Due to the linearity assumption, we can represent the interface at t = 0 as a superposition of a continuum of Fourier-Bessel modes. The consequent evolution of each mode of wavenumber k, follows from 3.85 with an integral superposition at all later times. Note that since our theory is for radially unbounded geometry, all possible value of wavenumbers are allowed. These results has been presented in Farsoiya and Dasgupta (2017a,b). For an initial localised interfacial perturbation of the form h(r, 0), we define the zeroth order Hankel transform as h(k, 0). For a single Bessel mode, we have seen earlier the expression for ã(s) in equation 3.85 which represents the amplitude of the Bessel mode in the Laplace domain. When a continuum of modes is excited due to a localised interfacial perturbation, the expression for h(k, s) (the tilde indicates Laplace transformed variables) is

H [h(r, t)] = Z L 0 rh(r, t)J 0 (kr)dr, ( 3
h(k, s) = h s + 1 r U +r L n 2k 2 (µ U + µ L ) + k(µ U z + µ L x) 2k 2 ⇣ µ U z l U +k + µ L x l L +k ⌘oi s 2 + 1 r U +r L n 2k 2 (µ U + µ L ) + k(µ U z + µ L x) 2k 2 ⇣ µ U z l U +k + µ L x l L +k ⌘o s + w 2 0 h(k, 0) (3.103)
where

w 2 0 = gk ⇣ r L r U r L +r U ⌘ + Tk 3
r L +r U and z(s) and x(s) are given by

z(s) ⌘ 2(k + l U )(µ U k + µ L l L ) µ L (k + l L ) + µ U (k + l U ) , x(s) ⌘ 2(k + l L )(µ L k + µ U l U ) µ L (k + l L ) + µ U (k + l U ) . l L = r k 2 + s n L , l U = r k 2 + s n U
From linear superposition, one may thus obtain the interface shape at any later time using inverse Hankel and inverse Laplace transform as

h(r, t) ⌘ Z • k=0 kdk J 0 (kr) L 1 { h(k, s)} (3.104)
The Laplace and the inverse Hankel transforms are computed numerically in

MATLAB (2015).

Results

The in-house solver is used to verify the above expressions. An initial localised perturbation given by h(r, 0) = a 0 e r 2 (1 r 2 ) (taken to be of the same form as Miles (1968)), with approximately air-water properties, r U = 0.001g cm 3 , µ U = 10 4 P, r L = 1g cm 3 , g = 981 cm s 2 , T = 72 dyne cm 1 . As seen in figure 3.14, analytical predictions agree quite well with DNS. Note that the simulation is taken only upto time when outgoing waves are not reflected from the computational domain. When this happens, mismatch between theory (assumes radially unbounded domain) and DNS which has a finite computational domain, are to be expected and are indeed seen.

Conclusion

In this chapter, we have solved the linearised initial value problem of an axisymmetric perturbation initialised as a Bessel mode (or as a localised perturbation) on the interface of two radially unbounded, viscous immiscible fluids. We find that a small amplitude Bessel mode evolves temporally in the same way as a spatial Fourier mode in planar geometry. We compare our analytical results to the normal mode approximation and to DNS. For small amplitudes, the DNS results compare extremely well to the solution of the IVP. The normal mode approximation is found to be good at very early and late times with deviations at intermediate times. The vorticity jump across the interface is found to depend in general both on the density and viscosity ratios of the two fluids. We show that by choosing the viscosity ratio to be unity, the magnitude of this jump can be reduced. Upon increasing the initial amplitude of perturbation, nonlinear effects produce qualitatively new features not present in the linearised IVP solution. Notable among these is the formation of a jet at the axis of symmetry which can rise to a height much greater than the parent perturbation. Depending on the viscosity ratio, the jet may undergo end pinch-off (i.e. pinch-off from the tip of the jet) while advancing or receding, leading to daughter droplet formation. A relation between the maximum jet height and the initial perturbation amplitude is obtained. Parametric studies are reported for the effect of change of parameters on underdamped oscillations. The emergence of new modes is studied using the Hankel transform of the interface. The analysis used for single Bessel mode is used to solve axisymmetric Cauchy-Poisson problem (for arbitrary initial conditions) in closed form, taking into account viscosity and inertia in both the fluids. For general initial conditions which excite a large (an uncountable infinite) number of modes at t = 0, superposition of which produces a travelling wave. An excellent agreement is shown by analytical expressions with DNS.

It is worthwhile to discuss the importance of the continuous spectrum obtained in our study. While we have only studied initial conditions without any vorticity, practical applications e.g. in oceans can involve localised regions of vorticity (turbulence) within one wavelength depth. The discrete eigenfunctions found in this study decay exponentially fast and are hence not adequate for representing such vortical initial conditions. It is necessary to solve the IVP for such initial conditions where the continuous spectrum is expected to have a large contribution.

Azimuthal oscillations of cylindrical gas filament surrounded by viscous fluid

Some the results presented in this chapter have appeared in "Azimuthal capillary waves on a hollow filament-discrete and the continuous spectrum", Palas Kumar Farsoiya, Anubhab Roy and Ratul Dasgupta, J. Fluid. Mech, 2020, 883, A21

We acknowledge collaboration with Prof. Anubhab Roy, IIT Madras on this work.

Introduction

In the previous chapter, we have analysed the discrete and continuous spectrum for capillary-gravity waves on a radially unbounded cylindrical pool of infinitely deep liquid. In this chapter, we demonstrate that a similar continuous spectrum also exists for a cylindrical filament comprising of incompressible gas surrounded by viscous liquid, assumed to be radially unbounded

w 2 0 = T R 3 0 2 4 kR 0 k 2 R 2 0 + m 2 1 r I I m (kR 0 ) I 0 m (kR 0 ) r O K m (kR 0 ) K 0 m (kR 0 ) 3 5 . (4.1)
Here w is the (circular) frequency, k is the wavelength of the interfacial, axial Fourier mode (k 2 R + ), m is the index of the interfacial, azimuthal Fourier mode (m 2 Z + ), R 0 is the radius of the unperturbed filament, T is surface tension and I m (•) and K m (•) are the mth order, modified Bessel functions of the first and second kind respectively. Note from 4.1 that while the RP instability criterion (kR 0 < 1) is the same for a liquid and a hollow filament, the most unstable wavelength (and its temporal growth rate) turns out to be higher for the hollow filament compared to a liquid filament of the same density (Rayleigh, 1892a).

As is also seen from 4.1, the RP instability may be bypassed rather easily.

As the denominator of 4.1 is always positive, imposing a three dimensional (µ cos(mq) cos(kz), m > 0, k 0) or even a purely azimuthal interfacial perturbation (k = 0, m 2) on the filament, produces simple harmonic oscillations at the interface. The case of purely azimuthal perturbations (kR 0 ! 0, m > 0) is particularly interesting. For such perturbations, the dispersion relation 1998)). Akin to inviscid results, in the viscous case as well, the wavelength of the most unstable mode (and its growth-rate) is found to be higher for the hollow filament (Chandrasekhar, 1981) when compared to a liquid filament of the same density and viscosity.

w 2 0 = T R 3 0  m(m 2 1) r I + r O . ( 4
An important feature of the temporal spectrum in these viscous problems on unbounded domains, is that it may not be purely discrete. For a given Fourier mode imposed on the interface, the viscous dispersion relation may admit finite number of roots only, thus yielding a finite set of eigenfunctions. These eigenfunctions constituting the discrete spectrum, typically decay exponentially in space.

Consequently any perturbation (e.g. in vorticity) imposed in the far field cannot be expressed as a linear combination of the discrete eigenmodes alone. In other words, the discrete spectrum eigenfunctions do not form a complete set and there is a continuum of eigenfunctions (labelled as the continuous spectrum), whose contribution may be large depending on fluid parameters and initial conditions. This is a common feature in viscous, interfacial perturbation problems occurring on unbounded domains e.g. capillary standing waves occurring on a rectangular The in-house code using which DNS is conducted is described in chapter 2. We compare our analytical predictions to DNS in section 4.5. For small amplitude perturbations, excellent agreement is demonstrated between analytical solutions and DNS. The chapter concludes with a summary and scope of further work.

(

Governing equations and boundary conditions -Hollow filament

As seen in figure 4.1a and 4.1b, the base state consists of an infinitely long cylindrical gaseous filament of radius R 0 . The fluid outside the filament is taken to be radially unbounded and is quiescent in the base state with density and kinematic viscosity r and n respectively. We ignore the dynamic effect of the gas inside the filament assuming it to be a low density, low viscosity, incompressible quiescent gas, exerting negligible stress at the interface (hollow filament approximation). The interface is thus free of tangential stresses at all times and we solve only for the fluid outside the filament. We take into account surface tension while neglecting gravity in subsequent analysis. Subject to these approximations, we solve for the flow generated in the fluid outside the filament due to linearised perturbations taken to be of the standing wave form. The incompressible Navier-Stokes and continuity equations govern the total velocity u tot and total pressure fields p tot viz.

Du tot Dt = 1 r rp tot + nr 2 u tot , r • u tot = 0, (4.3)
All variables (with superscript tot standing for total) may be written into a sum of base and perturbation variables as,

u tot = 0 + u, p tot = T R 0 + p (4.4)
where T is the gas-liquid surface tension and R 0 is the unperturbed radius of the circular filament. As we do not solve for the gas inside the filament, the base state pressure inside the filament is taken to be zero (gas inside exerts negligible pressure) in our analysis implying that the base state pressure in the fluid outside the filament is negative. This is implied in the second of equations 4.4. Linearised equations governing the perturbations may be obtained using the standard procedure of substituting 4.4 into 4.3, subtracting the base state

∂u ∂t = 1 r rp + nr 2 u, r • u = 0, (4.5)
As the perturbation treated here is pure azimuthal, our analysis is twodimensional and we choose to solve equations 4.5 in the streamfunction-vorticity formulation in cylindrical plane polar coordinates. The curl of the momentum equation in 4.5 leads to the equation for vorticity (w ⌘ r ⇥ u). Expressing velocity in terms of a vector streamfunction y using u ⌘ r ⇥ y (the negative sign is for convenience) and using standard vector identities, leads to the following equations written in the (perturbation) streamfunction-vorticity formulation.

∂w ∂t = nr 2 w, r 2 y = w. (4.6)
For two-dimensional flow the only non-zero component of vorticity and streamfunction is the out-of-plane z component i.e. w = (0, 0, w(r, q, t)) and y = (0, 0, y(r, q, t)) (see figure 4.1b). Projecting equations 4.6 into a cylindrical system (r-q-z), we obtain the z component of these equations viz.

∂w ∂t = n ✓ ∂ 2 w ∂r 2 + 1 r ∂w ∂r + 1 r 2 ∂ 2 w ∂q 2 ◆ (4.7) ∂ 2 y ∂r 2 + 1 r ∂y ∂r + 1 r 2 ∂ 2 y ∂q 2 = w, (4.8)
We may obtain the radial and azimuthal components of the two-dimensional perturbation velocity field u = (u r , u q , 0) from the perturbation streamfunction y(r, q, t) using u r = 1 r ∂y ∂q , u q = ∂y ∂r .

These relations follow from the definition u ⌘ r ⇥ y used earlier in obtaining 4.6. In order to impose boundary conditions, we define a scalar field F and a unit normal and tangent n and t respectively (forming a right handed coordinate system), to the perturbed interface. where the stress strain-rate relation, s tot = p tot I + µ " ru tot + ru tot T # applies to a Newtonian fluid. Projecting equations 4.10-4.12 on a plane polar (r, q) coordinate system, using the decomposition 4.4 and expressing (perturbation) velocity components in terms of (perturbation) streamfunction while retaining upto linear terms in perturbed quantities we obtain,

F(r, q, t) ⌘ r R 0 h(q, t), n ⌘ rF |rF| . ( 4
∂h ∂t + ✓ 1 r ∂y ∂q ◆ r=R 0 = 0, (4.13) µ ✓ ∂ 2 y ∂r 2
1 r ∂y ∂r

1 r 2 ∂ 2 y ∂q 2 ◆ r=R 0 = 0, (4.14) p(R 0 , q, t) + 2µ ✓ 1 r ∂ 2 y ∂r∂q 1 r 2 ∂y ∂q ◆ r=R 0 = T R 2 0 ✓ h + ∂ 2 h ∂q 2 ◆ , (4.15)
and lim

r!• w(r, q, t) ! finite, lim r!• y(r, q, t) ! finite. (4.16)
In further analysis, we use governing equations 4.7 and 4.8 along with boundary conditions 4.13-4.16.

Linear stability analysis -Normal modes

Here we conduct a temporal stability analysis on the base state described earlier.

An azimuthal Fourier mode of the form cos (mq) is imposed at the interface (see figure 4.1a) and we solve for the eigenvalues and the radial part of the eigenfunctions, using the normal mode approach. This implies setting the temporal dependency for all field variables to be of the form exp (st). It will emerge from this analysis that for a given value of m and Laplace number La ⌘ TR 0 r µ 2 , the problem admits a discrete spectrum with a finite set of poles on the complex s plane as well as a continuous spectrum, represented by the negative s axis. The discrete part of the spectrum is analysed in the following sub-section.

Normal modes -Discrete spectrum

Due to variable separability, we seek normal mode solutions in the form of standing waves and set, h(q, t) = a 0 cos(mq) where c.c. stands for complex conjugate. Here W(r), Y(r) and P (r) are the eigenfunctions of the radial part of the Laplacian operator (written in cylindrical coordinates) while s is related to it's eigenvalue. We assume that a 0 and m are real (the latter restricted to only integer values for periodicity, m 2 Z + ) while s is allowed to be complex (temporal analysis). Due to s being complex, W(r), Y(r) and P (r) are complex functions of a real argument as will be seen in subsequent algebra. Substituting


W(r) = CK m (lr) + DI m (lr) (4.25)
where C, D are (complex) constants of integration and I m , K m are mth order, modified Bessel functions of the first and second kind respectively. We define l ⌘ q s n , noting that l is like a wavenumber having the dimensions of inverse length. In subsequent algebra we write our equations in terms of l, replacing all instances of s with nl2 . In order to prevent divergence of W(r) as r ! •, we set D = 0 in 4.25. We are thus implicitly assuming that <(l) > 0 2 , since for fixed m and z ! •, I m (z) is asymptotic to exp(z) p 2pz

(National Institute of Standards and Technology, 2018) which diverges as z ! • only if <(z) > 0. With W(r) = CK m (lr) from 4.25, the solution to equation 4.24 can be written as a linear combination of the two independent homogeneous solutions v 1 (r) = r m and v 2 (r) = r m and the particular integral (see equation 2.2.13 in Prosperetti (2011). This is,

Y(r) =  a + C Z r • q m+1 K m ( q) 2m dq r m +  b C Z r R 0 q m+1 K m ( q) 2m dq r m (4.26)
Here a and b are real constants of integration (since coefficients of the left hand side of equation 4.24 are real) whose value depends on the choice of the lower limits of integration and q ⌘ q q s n = ql. We set a = 0 in equation 4.26 to prevent divergence as r ! • and obtain the following expression for Y(r) after lengthy algebraic manipulations.

Y(r) = br m + ✓ C l 2 ◆ K m (r) C R 0 2ml ✓ r R 0 ◆ m K m+1 R0 (4.27)
with r ⌘ rl, R0 ⌘ R 0 l etc. In order to satisfy boundary conditions, we need an expression for perturbation pressure p(r, q, t). This is obtained from the (linearised) momentum equation for the radial component of velocity. Expressing the linearised radial momentum equation in terms of streamfunction y and using expression 4.27 and 4.19, we obtain

1 r ∂p ∂r = m cos(mq) " 1 2
( nl 2 br m 1 nlC 2m ✓ r R 0 ◆ m 1 K m+1 R0 ) exp(nl 2 t) + c.c. # (4.28)
Equation 4.28 can be integrated from r to • with the boundary condition p(•, q, t) = 0. We thus obtain the following expressions for the radial part of the field variables as defined in equations 4.18-4.20.

W(r) = CK m (r) (4.29) Y(r) = br m + ✓ C l 2 ◆ K m (r) C R 0 2ml ✓ r R 0 ◆ m K m+1 R0 (4.30) P (r) = r ( nl 2 br m + nlC R 0 2m ✓ r R 0 ◆ m K m+1 R0 ) (4.31)
Expressions 4.29-4.31 when used in boundary conditions 4.13-4.15 lead to the following homogeneous equations.

nl 2 a 0 + m R 0 Y(R 0 ) = 0 (4.32) ✓ d 2 Y dr 2 1 r dY dr + m 2 r 2 Y ◆ r=R 0 = 0 (4.33) P (R 0 ) + 2µm ✓ 1 r dY dr 1 r 2 Y ◆ r=R 0 Ta 0 R 2 0 ⇣ 1 m 2 ⌘ = 0 (4.34)
Substituting expressions for Y(r) and its derivatives as well as P (r) from 4.29-4.31 into equations 4.32-4.34 and after lengthy algebraic manipulations, we obtain three homogeneous equations in a 0 , C and b. These are

nl 2 R 0 a 0 K m 1 ( R0 ) 2 R0 C + mR m 2 0 b = 0 (4.35) K 0 m 1 ( R0 ) C + 2m(m + 1)R m 2 0 b = 0 (4.36) T rR 2 0 (1 m 2 )a 0 + n R0 " R2 0 2m K m+1 ( R0 ) (m 1)K m 1 ( R0 ) # C nl 2 R m 0 1 + 2m(m + 1) R2 0 ! b = 0 (4.37)
For obtaining non-trivial values of a 0 , C and b, we must have

nl 2 K m 1 ( R0 ) 2l m R 0 0 K 0 m 1 ( R0 ) 2m(m + 1) R 2 0 T rR 2 0 (1 m 2 ) n R0 G( R0 ) nl 2 1 + 2m(m + 1) R2 0 ! = 0, where G( R0 ) ⌘ R2 0 2m K m+1 ( R0 ) (m 1)K m 1 ( R0 ).
Solving the determinant and after some further algebraic manipulations, we obtain the viscous dispersion relation for azimuthal perturbations expressed in non-dimensional form.

R4 0 + ( 2m(m + 1) 1 G( R0 ) R0 K 0 m 1 ( R0 ) ! ) R2 0 + La m(m 2 1) 1 (m + 1) R0 K m 1 ( R0 ) K 0 m 1 ( R0 ) ! = 0 (4.38)
where the Laplace number is defined as La ⌘ TR 0 r µ 2 . Equation 4.38 may be interpreted as a (non-dimensional) dispersion relation which determines R0 for given values of La and m. A number of consistency checks have been performed on this equation. In Appendix I.3, we show that it reduces to the viscous dispersion relation (equation I.49) in Cartesian coordinates (Lamb, 1932) under the limit m !

•, R 0 ! • such that m/R 0 ! k,
where k is a constant and represents wavenumber in Cartesian coordinates. In Appendix I.4, we expand equation 4.38 in powers of n 1/2 for small n (large La) showing that it reduces to the inviscid limit 4.2. The first viscous correction occurs at O(n) and provides the damping rate for the oscillation. For arbitrary values of Laplace number, 4.38 may be solved numerically to determine R0 for given values of m and La. Note that the relation Km (z) = K m ( z) (Abramowitz and Stegun, 1965) implies that roots of equation 4.38 when complex always occur in conjugate pairs. For given values From this pair and using the relation

s = n ✓ R0 R 0 ◆ 2
, we obtain a corresponding conjugate pair of values for s. Table 4.1 provides a sample list of such pairs. In figures 4.3a and 4.3b, we show such complex conjugate pairs of values of ŝ ⌘ s/w 0 as a function of La for m = 2 and 3 respectively. It is seen that <( ŝ) is always negative implying that all perturbations are damped. For La! •, the curves representing the locus of all pairs of roots start from s = ±Iw 0 (given by 4.2) and with increasing Laplace number the curves turn around touching the x axis asymptotically as La! 0. This is verified from 4.38 where for La= 0, R0 = 0 is a root of the equation. As the curves in figures 4.3a and 4.3b do not touch the real axis for any non-zero value of La, this implies that the response of the filament always remains underdamped. Such lack of overdamped behaviour has also been reported earlier for capillary oscillations on a spherical gas bubble (Prosperetti, 1980a) and shows that a hollow filament has many similarities to a Serial no. bubble. This analytical prediction of the lack of an overdamped response at large viscosity (i.e. small Laplace number) is tested using DNS.

r µ T R 0 m La⌘ TR 0 r µ 2 R0 l = R0 R 0 s = ✓ R0 R 0 ◆2 n 1 1 p 2 0.

Normal modes -continuous spectrum

In the previous section, we have numerically found that for a given value of m and Laplace number La, the dispersion relation 4.38 allows only two values of l and correspondingly s (a complex conjugate pair). This complex conjugate pair leads to two radial vorticity eigenfunctions viz. the real and the imaginary parts of K m (lr). These constitute the discrete spectrum (DS hereafter) eigenmodes. It is clear that these two eigenfunctions which decay exponentially at large radii, cannot express far field vorticity perturbations. In order to ensure completeness, we need another set of eigenfunctions. We show here that (for the same m and La as above) there is an additional, uncountably infinite set of eigenfunctions labelled as the continuous spectrum (CS hereafter). Note that in the previous section, we had assumed that <(l) > 0. We now allow for < (l) = 0 in our analysis writing l = q s n ⌘ Ix with x 2 R + (see Lamb (1932) for a similar argument in Cartesian geometry). This implies s is a real (negative) number in contrast to earlier analysis where it was complex. We thus set h(q, t) = a 0 exp(st) cos(mq), (4.39) w(r, q, t) = exp(st) sin(mq)W(r), (4.40) y(r, q, t) = exp(st) sin(mq)Y(r), (4.41) p(r, q, t) = exp(st) cos(mq)P (r) (4.42) and replace s in subsequent algebra with nx 2 . In order to avoid profusion of symbols in the analysis, we retain the same symbols for all variables and constants of integration as in the previous section. In terms of x, equations 4.23 and 4.24 may be written as

d 2 W dr 2 + 1 r dW dr + ✓ x 2 m 2 r 2 ◆ W = 0, d 2 Y dr 2 + 1 r dY dr m 2 r 2 Y = W, (4.43)
respectively. The solution to the first equation in 4.43 is

W(r) = CJ m (xr) + DY m (xr) (4.44)
where C, D are real constants of integration (as coefficients of equation 4.39 are real) and J m , Y m are Bessel functions of the first and second kind.

Both J m and Y m decay to zero as r ! •. and thus unlike the discrete spectrum analysis, it is not necessary to set C or D to zero in 4.44. The solution to the second equation in 4.43 with W(r) from 4.44 can be obtained using steps which are identical to the previous section. The expression for Y(r) satisfying boundedness constraints at r ! • is

Y(r) = br m 1 x 2 ( CJ m (r) + DY m (r) ) + R 0 2mx ( CJ m+1 ( R0 ) + DY m+1 ( R0 ) ) ✓ R 0 r ◆ m , (4.45)
where r ⌘ xr, R0 = xR 0 . The expression for pressure is obtained by integrating the radial component of the momentum equation. After some algebra we obtain,

P (r) = r " sb + R 0 2xm ( CJ m+1 ( R0 ) + DY m+1 ( R0 ) )# ✓ R 0 r ◆ m (4.46)
Note that we now have four unknowns viz. C, D, a 0 and b with three boundary conditions (like earlier) to determine these. Substituting expressions 4.39-4.42 in 4.13-4.15 and using expressions for W(r), Y(r) and P (r) from 4.44, 4.45 and 4.46 respectively, we obtain the following linear equations. Due to having three equations for four unknowns (a 0 , b, C, D), we can only determine three ratios from these equations. We have chosen these to be a 0 /b, C/b and D/b. These are,

a 11 ✓ a 0 b ◆ + a 12 ✓ C b ◆ + a 13 ✓ D b ◆ = mR m 1 0 , (4.47) a 21 ✓ a 0 b ◆ + a 22 ✓ C b ◆ + a 23 ✓ D b ◆ = 2m(m + 1)R m 2 0 , (4.48) a 31 ✓ a 0 b ◆ + a 32 ✓ C b ◆ + a 33 ✓ D b ◆ = rnx 2 (4.49)
with the coefficients given by

a 11 ⌘ nx 2 , a 12 ⌘ 1 x ( 1 2 J m+1 ( R0 ) m R0 J m ( R0 ) ) , a 13 ⌘ 1 x ( 1 2 Y m+1 ( R0 ) m R0 Y m ( R0 )
)

a 21 ⌘ 0, a 22 ⌘ J 0 m 1 ( R0 ), a 23 ⌘ Y 0 m 1 ( R0 ) a 31 ⌘ T R 2 0 ⇣ 1 m 2 ⌘ , a 32 ⌘ ✓ 2µm R0 J m 1 ( R0 ) rR 0 2xm J m+1 ( R0 ) ◆ , a 33 ⌘ ✓ 2µm R0 Y m 1 ( R0 ) rR 0 2xm Y m+1 ( R0 ) ◆
In contrast to the analysis in section 4. is obtainable for every value of 0 < x < • by solving equations 4.47-4.49. Note that since s = nx 2 , this implies that • < s < 0 for the continuous spectrum.

Eigenfunctions and completeness

Having demonstrated the existence of the discrete and continuous spectrum in sections 4.3.1 and 4.3.2, we now examine the radial part of the vorticity eigenfunctions. For a given value of m and La, there is a complex conjugate pair of l, l obtained from the dispersion relation and correspondingly the DS eigenmodes for the radial part of vorticity (equation 4.25 with D = 0.) K m (lr), K m ( lr) (4.50)

These modes decay to zero exponentially fast as r ! • since K m (z) ⇠ q p 2z exp ( z) as z ! •. Figure 4.4a and 4.4b shows these eigenfunctions plotted as a function of radius for two different values of La. It is seen that the length scale of decay of the eigenmodes increases with increasing viscosity. In addition we also have the CS eigenmodes which decay algebraically to zero. 1980b). The completeness theorem (for the radial part of vorticity field) may be written as a linear superposition of the DS and the CS modes. For an interfacial Fourier mode of index m, the radial part of any arbitrary initial distribution of vorticity may be (formally) represented as

W(r, 0) = Â i C i K m (l (i) r) + Z • 0 ⇢ D(x)J m (xr) + E(x)Y m (xr) xdx (4.52)
where C i , D(x) and E(x) are to be obtained from the inner product of W(r, 0) with K m , J m and Y m using suitable orthogonality conditions. The integral in 4.52 is purposely written in the form of an inverse Fourier-Bessel transform. From normal mode analysis, the time evolution of W(r, t) may be written as

W(r, t) = Â i C i K m (l (i) r) exp(s (i) t) + Z • 0 ⇢ D(x)J m (xr) + E(x)Y m (xr) x exp( nx 2 t)dx (4.53)
where it is seen that 4.52 is recovered from 4.53 at t= 0. The integral in equation 4.53 can produce terms with non-exponential time dependence, a contribution which is missed if the continuous spectrum is not taken into account. It is interesting to note is that if the initial condition has a projection on a single continuous spectrum eigenmode only, i.e. D(x), E(x) µ d(x x 0 ) then W(r, t) becomes expressible as a sum over exponentials. A non-exponential contribution to W(r, t) from the continuous spectrum arises only if D(x) has support on a range of x. It will be seen in the next section, that perturbing the interface with a single Fourier mode and zero perturbation vorticity everywhere in space, represents such an initial condition. Similar to 4.52 which is a decomposition in space, it is also possible to write a corresponding completeness theorem for any function of time by expressing it as a linear combination of the discrete spectrum frequencies (allowed by the dispersion relation) and an integral superposition over the continuous spectrum frequencies e.g. the amplitude of the standing wave perturbation is expected to be

a(t) a 0 = Â i k i exp(s (i) t) + Z 0 • X(s) exp(st)ds (4.54)
where k i and X(s) are to be obtained from initial conditions. Note that the s i in equation 4.53 and 4.54 are obtained from the l (i) using the relation

s (i) = n ⇣ l (i) ⌘ 2
ensuring that only those l (i) satisfying < ⇣ l (i) ⌘ > 0 are used. Equations 4.53 and 4.54 are written based on the normal mode analysis presented here and in section 4.3.1. By solving the corresponding initial value problem (IVP), it will be shown in section 4.4 that one can analytically obtain expressions of exactly the same form as 4.53 and 4.54. Furthermore, the IVP solution will provide expressions for C i , D(x), E(x) in 4.53 as well as k i , S(s) in 4.54.

Initial value problem (IVP)

In this section, we solve the IVP with initial perturbation comprising a single interfacial azimuthal Fourier mode of index m with zero initial velocity (vorticity) everywhere in the outer fluid. We further assume that there is no impulse at the interface initially. For m 2 and La > 0, it is seen from equation 4.50 and 4.51 that the vorticity field associated with the DS and CS eigenfunctions decays to zero exponentially and algebraically respectively. Consequently the aforementioned initial condition we solve for here, requires a linear superposition of the DS as well as the CS eigenfunctions, to achieve a cancellation of vorticity everywhere in space. This implies that the initial condition excites the DS and the CS modes. We consequently anticipate the temporal evolution of the interface to be a summation of discrete exponential terms and an integral over a continuous range of exponential terms (see Prosperetti (1976) for similar conclusions for viscous capillary waves on a deep pool of liquid). For solving the IVP, we set h(q, t) = a(t) cos(mq), (4.55) w(r, q, t) = sin(mq)W(r, t), (4.56) y(r, q, t) = sin(mq)Y(r, t) (4.57) p(r, q, t) = cos(mq)P (r, t) (4.58)

Substituting expressions 4.55-4.58 in equation 4.7 and 4.8, we obtain

∂W ∂t = n  ∂ 2 W ∂r 2 + 1 r ∂W ∂r m 2 r 2 W (4.59) ∂ 2 Y ∂r 2 + 1 r ∂Y ∂r m 2 r 2 Y = W (4.60)
Define the Laplace transform and inverse transform pair f (t) and f (s) as,

f (s) ⌘ Z • 0 f (t) exp( st)dt (4.61)
With initial conditions described earlier in this section, we have a(0) = a 0 , . a(0) = 0, W(r, 0) = 0 and the Laplace transform of equations 4.59 and 4.60 is (Laplace transformed variables are indicated with a tilde on top)

∂ 2 W ∂r 2 + 1 r ∂ W ∂r ✓ s n + m 2 r 2 ◆ W = 0, ∂ 2 Ỹ ∂r 2 + 1 r ∂ Ỹ ∂r m 2 r 2 Ỹ = W(r, s) (4.62)
The algebra henceforth is entirely similar to that in section 4.

Ỹ(r, s) = b(s)r m + ✓ C(s) h 2 ◆ K m (hr) C(s)R 0 2mh ✓ r R 0 ◆ m K m+1 (hR 0 ) (4.64)
where h ⌘ p s n . C(s), D(s) are complex functions while b(s) is a real function of s. Equations 4.63 and 4.64 are obtained under the constraint <(h) > 0 in order to prevent divergence as r ! •. We discuss this point when discussing the Laplace inversion in Appendix I.7. Using steps identical to the discrete spectrum derivation, we can obtain an expression for the radial part of perturbation pressure in the Laplace domain. This is,

P (r, s) = r ( s b(s)r m + s C(s)R 0 2mh ✓ r R 0 ◆ m K m+1 (hR 0 ) ) (4.65)
Laplace transforming the boundary conditions 4.13-4.15 and using expressions for W(r, s), Ỹ(r, s) and P (r, s) from 4.63-4.65 into these, we obtain three inhomogeneous equations in the three unknowns viz. ãm (s), C(s) and b(s). These are

s R 0 ã(s) K m 1 (hR 0 ) 2hR 0 C(s) + mR m 2 0 b(s) = a 0 R 0 (4.66) K 0 m 1 (hR 0 ) C(s) + 2m(m + 1)R m 2 0 b(s) = 0 (4.67) T rR 2 0 (1 m 2 ) ã(s) + n hR 0 " h 2 R 2 0 2m K m+1 (hR 0 ) (m 1)K m 1 (hR 0 ) # C(s) sR m 0 1 + 2m(m + 1) h 2 R 2 0 ! b(s) = 0 (4.68)
The expression for ã(s) can be obtained from 4.66-4.68. After some algebra, it maybe shown that this is (written in the same form as earlier in chapter 3),

ã(s) = " s + M(s) s 2 + sM(s) + w 2 0 c(s) # a 0 , < ✓r s n ◆ > 0 (4.69) with G(s) ⌘ s n R 2 0 2m K m+1 ✓r s n R 0 ◆ (m 1)K m 1 ✓r s n R 0 ◆ , c(s) ⌘ 1 (m + 1) p s n R 0 K m 1 p s n R 0 K 0 m 1 p s n R 0 , M(s) ⌘ 2nm(m + 1) R 2 0 1 G(s) p s n R 0 K 0 m 1 p s n R 0 ! and w 2 0 ⌘ Tm(m 2 1) rR 3 0
Note that w 0 is the inviscid frequency discussed in the introduction in equation 4.2, under the approximation r O >> r I .

In order to obtain an expression for a(t), we need the Laplace inversion of 4.69.

For this it is necessary to know the poles and branch points (and the consequent branch cut structure) of ã(s). The possibility of ã(s) having branch points arises due to K m (z) having a branch point at z=0 (and at z = •) (Abramowitz and Stegun, 1965) as seen from the following series representation (Parnes, 1972),

K m (z) = 1 2 m 1 Â r=0 ( 1) r (m r 1)! r! ✓ 2 z ◆ m 2r + ( 1) m+1 1 m! ⇣ z 2 ⌘ m h ln ⇣ z 2 ⌘ 1 2 Y(m + 1) + 1 2 g + O ✓ ⇣ z 2 ⌘ m+2 ◆ (4.70)
where as can be seen in the series expansions provided in Appendix I.6. We depict the imaginary parts of ã(s) given by 4.69 for m = 2, 3 in figures 4.5a and 4.5b. It is clearly seen that the negative real axis is a line of discontinuity for

Y(m + 1) = 1 + 1 2 + 1 3 + . . .
= ⇣ ã(s) a 0 ⌘
for m = 2 and 3 respectively. The peaks in these figures are due to the two poles of ã(s) and the domain of the imaginary axis (=(s)) has been restricted to enable visualisation of the branch cut. We now proceed to invert 4.69 for which we re-write it as,

ã(s) = L(s) Q(s) a 0 , Q(s) ⌘ s 2 + sM(s) + w 2 c(s), L(s) ⌘ s + M(s) (4.73)
As is to be expected, the denominator Q(s) in 4.73 is the dispersion relation 4.38 written in terms of s, with s replacing s. Denoting the roots of Q(s) = 0 as s i (i.e. roots satisfying the constraint < p s n > 0, see 4.69) and using Cauchy residue theorem, the formal inversion of 4.73 may be expressed as (using the keyhole contour discussed in Appendix I.7)

a(t) a(0) = 1 2pI Z c+I• c I• L(s) Q(s) exp[st]ds = Â i L(s i ) Q 0 (s i ) exp(s i t) 1 2pI Z • 0 ⇢✓ L(x exp(Ip)) Q(x exp(Ip)) ◆ ✓ L(x exp( Ip)) Q(x exp( Ip)) ◆ exp[ xt]dx (4.74)
with c > 0 and the prime in Q 0 denoting a derivative. Note that the ratio L(s)/Q(s) satisfies the relation

L(s) Q(s) = ✓ L(s) Q(s) ◆ .
Consequently, the expression inside the braces in the integral of equation 4.74 is a purely imaginary quantity implying 4.74 may be written as

a(t) a(0) = Â i L(s i ) Q 0 (s i ) exp(s i t) 1 p Z • 0 = ⇢✓ L(x exp(Ip)) Q(x exp(Ip)) ◆ exp[ xt]dx (4.75)
Comparing 4.75 with expression 4.54 with the replacement s = x, it is seen that both are the same expressions. Equation 4.75 serves to determine k i and X(s) in 4.54. Expression 4.54 was obtained from normal mode analysis while 4.75 obtained through the IVP solution, validates our analysis. In an entirely analogous manner, the inversion for W(r, s) in 4.63 (see expression for C(s) in Appendix I.6) may be accomplished using Cauchy residue theorem as

Ŵ(r, t) ⌘ R 0 2(m + 1)a 0 w 2 0 W(r, t) = Â i K m ⇣ r q s i n ⌘ K 0 m 1 (R 0 q s i n )Q 0 (s i ) exp(s i t) 1 p Z • 0 = ⇢✓ R(r, x exp(Ip)) Q(x exp(Ip)) ◆ exp[ xt]dx (4.76) with R(r, s) ⌘ K m (r p s n ) K 0 m 1 (R 0 p s n )
. It is shown in Appendix I.8 that the integral on the right hand side of 4.76 may be re-written in terms of Bessel functions J and Y allowing us to rewrite 4.76 as

Ŵ(r, t) = Â i K m ⇣ r q s i n ⌘ K 0 m 1 (R 0 q s i n )Q 0 (s i ) exp(s i t) 1 p Z • 0 A(x)J m r p x n B(x)Y m r p x n A(x) 2 + B(x) 2 exp[ xt]dx (4.77) with A(x) ⌘ ⇣ x 2 w 2 0 ⌘ Y m (a) + 2x 2 ⇢ Y 00 m (a) m(m + 1) a 2 Y m 2 (a) 2w 2 0 a Y m 1 (a) B(x) ⌘ ⇣ x 2 w 2 0 ⌘ J m (a) + 2x 2 ⇢ J 00 m (a) m(m + 1) a 2 J m 2 (a) 2w 2 0 a J m 1 (a)
where a ⌘ R 0 p x n . Comparison of expressions 4.77 with 4.53 with the replacement p

x/n = x shows that they are of the same form with 4.77 determining C i , D(x) and E(x) in 4.53. Expressions 4. 75 and4.77 are the central results of our analysis and they verify equations 4.54 and 4.53 respectively obtained earlier from normal mode analysis. It is interesting to note that although the Bessel function Y m diverges at the lower limit of integration x = 0 in 4.77, the integrand is not singular at the lower limit. In order to validate analytical predictions 4.75 and 4.77, we compare with DNS in the next section.

Results and discussion

The DNS results in this section are produced using the in-house solver described and benchmarked in chapter 2. The DNS geometry is shown in figure 4.6a. The simulation domain size L was chosen to be sufficiently large relative to the undisturbed filament radius R 0 , to effectively mimic a radially unbounded domain. For the present purpose, we found that choosing L/R = 10 was sufficient to eliminate finite domain effects on the time scale of our simulations. DNS is governed by five non-dimensional parameters viz. m, Laplace number La, density ratio r = r I /r Unlike the DS approximation, the NMA has the advantage that by construction it satisfies initial conditions for a(t). It is seen in figures 4.7a-4.8b that while the NMA is a better approximation at early times, the DS approximation improves over it at later stages showing less disagreement with DNS, even at lower values of La. We also observe that for La= 2, the interface continues to show underdamped behaviour (see figure 4.7a and 4.8a), consistent with our predictions from the dispersion relation 4.38 in section 4.3.1.

It is seen from figure 4.9a that nonlinear effects are quite pronounced at e = 0.5. Similarly, systematic effect of inertia of the gas inside is seen in figure 4.9b as r r is increased, with large deviations from the present theoretical predictions seen at r r = 1.

For comparing the vorticity fields, we note that our in-house solver is written in two-dimensional Cartesian coordinates (x-y) while the analytical expressions for Ŵ(r, t) (equation 4.77) have been obtained using plane-polar (r-q) coordinates (see figure 4.1a). In figures 4.10a and 4.10b we show the contour of the vorticity field w(r, q, t) obtained from DNS at t = p/w 0 . In both the figures, the vorticity field inside the filament (i.e. in the region 0  r  R 0 + h) has been removed to aid comparison with the analytical prediction from IVP in 4.11a and 4.11b.

The four lobes seen in these figures corresponds to an azimuthal dependence of sin(2q). Note that the sign is opposite as W(r, t) is negative at the instant of time (see figures four different instants of time. The parameters correspond to case 1 in table 4.2.

A good match is seen in every case. The slight disagreement seen between DNS and analytical prediction 4.77 near the interface r = R 0 + h ⇡ 1, is to be expected as the interface in DNS experiences a small amount of tangential stress exerted by the fluid inside the filament. This is not taken into account in the theory where the interface is treated as being free of tangential stress at all times (see boundary condition 4.11). It is also seen from figures 4.12b-4.12e that the discrete spectrum approximation to the vorticity field (DS) given by the summation terms alone in 4.77, follows the same qualitative trend as DNS but has a large quantitative mismatch.

Conclusions

In this chapter, we have analysed the temporal spectrum of linearised, azimuthal Fourier modes on a hollow cylindrical filament, surrounded by radially unbounded, quiescent, viscous fluid. The base state is found to be stable to all azimuthal Fourier modes of a standing wave form. Using normal mode analysis, we obtain the dispersion relation governing the discrete spectrum for the problem.

= p/4 at different times -Case 1 in table 4.2.
The dispersion relation is a transcendental equation whose numerical solution predicts that all azimuthal Fourier modes (investigated here) respond in an underdamped manner, independent of the value of fluid viscosity. Using normal mode analysis further, we show that in addition to the discrete spectrum, the problem also admits a continuous spectrum. The discrete spectrum eigenfunctions decay exponentially at large radii and thus cannot represent far field perturbations. A continuum of uncountably infinite eigenfunctions are also found which decay to zero algebraically at large radii and are essential for representing far field perturbations. The completeness theorem for the radial part of the vorticity field is expressed as a weighted sum over the discrete eigenmodes and an integral over the continuous ones. A similar completeness statement is also obtained for any function of time, including as a particular case, the amplitude of the standing wave. We further solve the linearised, initial value problem (IVP) for a perturbation in the form of a purely azimuthal Fourier mode with zero perturbation vorticity. Analytical expressions governing the time dependence of the amplitude of the standing wave and its associated vorticity field, are obtained in the Laplace s domain. Analysis of these expressions shows that in addition to poles these also have branch points on the complex s plane. Laplace inversion of these expressions using Cauchy residue theorem shows that the discrete spectrum comes from the residue at the zeroes (which manifest themselves as poles) of the dispersion relation while the continuous spectrum comes from the difference in contribution from either side of the branch cut. The zero vorticity initial condition investigated in this study has projections on the discrete as well as the continuous spectrum. Consequently, both set of eigenmodes are excited at time t=0. The subsequent motion of the interface resembles a standing wave whose amplitude is a sum of discrete exponentials (due to the discrete spectrum) and an integral over a continuous range of exponentials (due to the continuous spectrum). Explicit expressions are obtained for the amplitude and the vorticity field associated with a purely azimuthal Fourier mode. Theoretical predictions are compared to DNS conducted using an in-house developed code and excellent agreement is observed at small value of the nonlinearity parameter e. In contrast, the discrete spectrum approximation is found to be a poor approximation at early times, especially at O(1) Laplace numbers.

Faraday waves on cylindrical fluid filament (in collaboration with Sagar Patankar)

The results presented in this chapter have appeared in "Faraday waves on cylindrical fluid filament-generalised equations

Introduction

In 1954) to obtain the amplitude of forcing h at which the (flat) interface becomes linearly unstable to a standing wave of wavenumber k, for a given forcing frequency W. The stability chart on the amplitude of forcing and wavenumber (h-k) space display tongue shaped curves. For sufficiently small forcing amplitude, choosing parameters inside the tongues, implies an unstable oscillatory response with exponential growth. The frequency response of the surface is nW/2 with n = 1, 2, 3 . . . (Benjamin and Ursell, 1954), labelled as subharmonic and harmonic alternatively. Outside these tongues for finite forcing amplitude h, a standing wave perturbation produces a bounded, oscillatory (but not necessarily periodic) response (Benjamin and Ursell, 1954).

129

(in collaboration with Sagar Patankar) An important extension of Faraday waves to spherical base state geometry has been made in Adou and Tuckerman (2016) who show that in this geometry, the inviscid problem is also governed by the Mathieu equation whose coefficients depend only on the lower index l of the spherical harmonic Y m l (q, f). Extension of these inviscid results to the viscous case has been made in both Cartesian and spherical geometry (Adou and Tuckerman, 2016; Kumar, 1996; Kumar and Tuckerman, 1994) where using Floquet theory, it has been demonstrated that it is possible for the first unstable response to be harmonic rather than subharmonic. In addition to witnessing substantial theoretical progress, the subject of interfacial parametric instabilities has also seen applications in vibration induced atomisation (James et al., 2003) leading to novel technologies for nebulisation [START_REF] Tsai | Faraday instability-based micro droplet ejection for inhalation drug delivery[END_REF]. In this chapter, we study Faraday waves on a cylindrical filament surrounded by an inviscid immiscible liquid. The filament is assumed to be axially unbounded and is subjected to a radial, time periodic body force. We study the stability of this time-periodic base state showing that this is also governed by the Mathieu equation.

Mathieu Equation for Faraday waves on cylindrical fluid filament

As shown in figure 5.1, the cylindrical filament of radius R 0 is subject to radial time periodic body force in both the fluids. For the unperturbed filament, the tangential coordinates are z and q and the normal coordinate is r. A time periodic body force G(r, t) acts radially on the filament of radius R 0 and thus the radius of curvature in the base state R = R 0 . The interface separates two inviscid, immiscible, quiescent fluids of density r I (inner fluid) and r O with surface tension coefficient T. We study the (linear) stability of this time-periodic base state, by imposing an interfacial perturbation in the form of a standing wave. For studying parametric forced standing waves on the cylindrical filament (see figure 5.1 for geometry of forcing), we set G(r, t) = h cos(Wt) r R 0 as a linearly varying function of radial position r and acting radially. We follow Adou and Tuckerman (2016) in regularising the singularity (multivaluedness) on the axis of the filament r = 0 using a linear spatial variation of the body force i.e. G(r, t) µ r/R 0 .

Choosing a non-linear variation of the form G(r, t) µ (r/R 0 ) z with z = 2 . . . does not alter our results qualitatively, since in a linearised description the coefficient of equation 5.1 is evaluated at r = R 0 . The Mathieu equation in cylindrical geometry is (see derivation in Appendix I.2)

d 2 a dt 2 + 2 4 w 2 m (k) + r I r O hk cos(Wt) r I I m (kR 0 ) I 0 m (kR 0 ) r O K m (kR 0 ) K 0 m (kR 0 ) 3 5 a(t) = 0 (5.1)
where w 2 m (k) is given by

w 2 m (k) ⌘ T R 3 0 2 4 kR 0 k 2 R 2 0 + m 2 1 r I I m (kR 0 ) I 0 m (kR 0 ) r O K m (kR 0 ) K 0 m (kR 0 ) 3 5
(5.2) Equation 5.1 can be rewritten compactly as (in the same notation as Bender and Orszag (2013)),

d 2 a d t2 + (A + 2B cos( t)) a( t) = 0, (5.3)
where

t ⌘ Wt, A = ⇣ w m (k) W ⌘ 2 and B ⌘ h 2W 2 " (r I r O )k r I I m (kR 0 ) I 0 m (kR 0 ) r O K m (kR 0 ) K 0 m (kR 0 ) # . Equation 5.3
is the well-known Mathieu equation whose solution can either be obtained through numerical integration or through perturbative techniques. It is known that depending on parameters (i.e. values of B and A), there can be bounded or exponentially growing solutions to equation 5.

Results and discussion

We compare analytical predictions made in earlier sections, against numerical results obtained from direct numerical simulations (DNS) of the Euler equations. a(0) = 1 and . a(0) = 0 in all cases. A good match with linearised predictions is seen in both cases. In the unstable case in figure 5.5b, it is seen that the DNS amplitude follows an initial exponential growth but eventually saturates after few cycles, leading to disagreement with the solution to the Mathieu equation 5.3, which has no mechanism for saturation.

Case Comments e ax ⌘ a(0)k e az ⌘ a(0)m R 0 r r ⌘ r I r O G ⌘ r I hR 2 0 T b ⌘ r I R 3 0 W 2 T (A, B)

RP modes in the presence of forcing

For a given value of forcing, a range of modes which are RP unstable (i.e. m = 0, kR 0 < 1) in the absence of forcing, can be stabilised through forcing in a small range of forcing amplitudes. We test this prediction here. Figure 5.6a is a stability chart for the special region kR 0 < 1 (here 0 < k < 6.0 with R 0 = 1/6). The white region in this chart corresponds to unstable behaviour while the hashed region is predicted from linear theory to be stable. Case 7 in figure 5.6a has the same perturbation wavenumber as Case 1 but is theoretically predicted to be stabilised through forcing. An important point to note here is that the stabilisation of axisymmetric RP modes in figure 5.6a, is predicted to be purely through (axisymmetric) forcing i.e. the forcing function G(r, t) in our simulations is independent of q. This is in contrast to the unforced case (h = 0), where the RP instability can be suppressed only through threedimensional initial perturbations (m 2). The predicted stabilisation of Case 7 is indeed observed in DNS at early times in figure 5.7a where both axisymmetric and 3D DNS, show stabilisation and agreement with the solution to the Mathieu equation. This agreement is however short-lived, as nonlinear effects start to (in collaboration with Sagar Patankar) appear around t ⇡ 2.5 producing instability. Interestingly, this instability can be understood from our linear theory. It is seen that at t = 2.638 (See Fig. 5.9a) the interface does not resemble a single Fourier mode, although we had started with a single wavenumber k = 4 at t = 0. The lower panel shows the Fourier transform f (k) of the interface in the upper panel, revealing the formation of two additional Fourier modes viz. k = 8 and 12 due to nonlinearity at t = 2.638. This is in contrast to the Fourier transform of the interface at t = 0 where there is a single peak corresponding to a Fourier mode k = 4. The inset of figure 5.9a shows the corresponding stability chart, with the additional wave numbers produced due to nonlinearity, shown as points. It is clearly seen that at this forcing amplitude, the k = 8 wavenumber lies inside a tongue and is unstable. This is the origin of the destabilisation seen in figure 5.7a. Note that while the time of onset of nonlinearity is roughly the same in both axisymmetric and 3D DNS in figure 5.7a, the subsequent temporal routes are very different as seen in the difference between time signals.

Figure 5.7b corresponds to Case 8 in the stability chart 5.6a and is predicted to be linearly unstable. Once again, a good match at early times is observed although around t ⇡ 3, nonlinear effects lead to disagreement with linear predictions. As a further validation of our linearised predictions, Case 9 and 10 in the stability chart 5.6a (for k = 2) are predicted to be linearly unstable and stable respectively, in the presence of forcing. This prediction is seen to be validated in figures 5.8a and 5.8b (see inset of both figures for linearised predictions). Once again at early times, both axisymmetric and 3D DNS show good agreement with the Mathieu solution but show disagreement at later times due to the production of higher modes and non-axisymmetric effects.

It is worthwhile to ask, if the predicted linear stabilisation of RP modes via forcing may be observed in experiments? We hypothesize that the effect of viscous damping will have a strong suppressing effect on the emergence of the higher wavenumbers and subsequent destabilisation seen in our inviscid DNS. Consequently, for a viscous fluid we anticipate that the stabilisation via forcing, predicted in this study may indeed be observed for much longer times than is seen in our inviscid results. Note that the stability boundaries of figure 5.6a will also be modified due to viscosity.

Parametric studies

For visualisation of the instantaneous surface of the cylindrical filament, we have provided a collage of three-dimensional images in figures 5.10, 5.11,5.12 and 5.13 which correspond to cases 6, 12, 15 and 16 in table 5.2 respectively. Figure 5.10 shows the temporal evolution of unstable Faraday waves at the cylindrical of the interface on the x-y plane does not resemble a cosine mode. Interestingly, the slice on the z-y plane shows that the interface remains circular to a good approximation. Figure 5.14c shows the instantaneous streamlines for the snapshot corresponding to non-dimensional time t = 14p. For ease of visualisation, the thickness of the streamlines have been increased giving them the appearance of tubes. It is seen that except at the necks of the filament, velocity vectors in the azimuthal (q) direction are small. Figure 5.11 shows a collage of instantaneous snapshots for Case 12 in table 5.2. This is also a linearly unstable case albeit for higher values of forcing parameter h compared to case 6. We see the clear emergence of sheet like structures after three forcing time periods ( t = 6p).

The disintegration seen at t = 6.8p is an outcome of the sheet closing upon the filament while shrinking, see movie in supplementary material in Patankar et al.

(2018) for details. The corresponding instantaneous slices on the x-y and z-y planes are seen in figure 5.15. The tiny dots seen at t = 6.6p inside and outside the filament in figures 5.15, correspond to droplets and bubbles formed due to pinch-off and entrapment of surrounding fluid inside the filament respectively, see movie. Due to a larger forcing amplitude, nonlinear effects are visible earlier in figure 5.11 when compared to 5.10. Note that both the simulations presented in figures 5.10 and 5.11 are initiated with axisymmetric initial perturbations (i.e. m = 0, k 6 = 0). The effect of three-dimensional initial perturbations (m 6 = 0, k 6 = 0) are seen in the collage provided in figures 5.12 and 5.13 where in addition to axial perturbation k, azimuthal perturbation of wavenumber m=2 and m=3 are also provided at time t=0 respectively. It is seen in both these cases that the instantaneous (nonlinear) patterns which appear, retain memory of the symmetry of the initial perturbation. Figures 5.12 and 5.13 panel c) shows two and three lobes respectively, retaining the memory of the initial perturbation. In both cases, the instability eventually leads to formation of sheets and their fragmentation, albeit in very different ways.

The corresponding instantaneous interface profiles on the x-y and z-y planes are shown in figures 5.16a, 5.16b and 5.17a,5.17b respectively. Figure 5.17c shows the instantaneous streamlines corresponding to t = 3p in figure 5.17a. Due to the initial three-dimensional nature of the perturbation, three-dimensional velocity fields are clearly seen. A systematic parametric study of the effect of change of non-dimensional parameters G, b, e az , e ax and r r (see table 5.2) are reported in figures 5.18a, 5.18b, 5.19a, 5.19b and 5.20 respectively. Figure 5.18a

shows the effect of change of non-dimensional forcing amplitude G ⌘ 11, 12, 13, 14 in table 5.2). Note that all the four simulations are linearly unstable and hence the amplitude grows with time. It is seen that increasing G causes nonlinear effects to appear relatively early. In approximately two (non-dimensional) forcing time periods i.e t = 4p, significant deviation from the Mathieu solution is observed for case 13, while case 14 shows disintegration after t ⇡ 10. Within the same time window, case 11 which has the lowest value of G continues to show good agreement with the Mathieu solution. 5.2. This is qualitatively similar to the RP instability, where introducing a non-zero azimuthal perturbation suppresses the instability. Both the unstable cases deviate from the Mathieu solution at approximately the same value of t ⇡ 8. Case 17 which is stabilised due to the higher value of m in the non-axisymmetric initial perturbation, displays a good agreement with DNS upto later times. Figure 5.19b shows that effect of change of e ax ⌘ a(0)k (Cases 5, 24, 25, 26 and 27 in table 5.2). These five simulations lie in the linearly stable regime and it is seen that the one with the smallest value of e ax agrees with the linearised solution to Mathieu equation for the longest time. It is clearly seen that for e ax > 0.6437, the agreement with linear predictions is rather poor. Figure 5.20 shows the effect of change of density ratio where we have varied density ratio by one order of magnitude (while restricting ourselves to r r > 1 for this study). All the three cases are in the linearly stable regime. It shows that change of density ratio has a small effect on the time of onset of nonlinearity. Linearised predictions from the Mathieu equation are seen to agree well over a change of an order of magnitude in the density ratio.

r I hR 2 0 T (cases

Conclusions

We have conducted a linear stability analysis of an interface separating two inviscid, immiscible, quiescent fluids subjected to a time periodic body force. It is found that the amplitude a(t) of standing wave interfacial perturbations, is governed by a Mathieu equation. This equation is similar to the known results on Faraday waves in two geometries viz. rectangular Cartesian (Benjamin and Ursell, 1954) and spherical (Adou and Tuckerman, 2016). The principal novel findings of the study can be summarised as follows. We also highlight potential applications at the end.

• A cylindrical capillary filament subject to an axisymmetric radial, timeperiodic, body force is shown to be able to sustain Faraday waves and parametric instability.

• The amplitude of standing waves on the filament is governed by the Mathieu equation in the linearised limit. For small amplitude perturbations and for small forcing amplitude, DNS results show very good agreement with the Mathieu solution.

• The stability chart of the Mathieu equation for the filament, contains a special region (m = 0, kR 0 < 1) which is Rayleigh-Plateau (RP) unstable in the absence of forcing. It is predicted that by choosing the forcing amplitude to lie in a certain range, RP modes can be stabilised.

• Stabilisation of RP modes via forcing is observed in DNS at early times. It is also seen that nonlinearity causes the emergence of higher wavenumbers, some of which are unstable due to forcing. These higher modes lead to eventual destabilisation of the filament. We show that this destabilisation can be understood from linear theory. Axisymmetric and three-dimensional DNS are found to adopt different routes to instability.

• Similar to the suppression of RP instability using azimuthal perturbations, we show that unstable three-dimensional Faraday waves can be stabilised by imposing perturbations of increasing azimuthal wavenumber.

• Parametric studies in a five dimensional space of numbers, clearly show the limits within which linearised theory yield accurate predictions. In particular, large values of e ax and G are found to be associated with early onset of nonlinear effects.

• Linearised theory is found to be unable to predict many interesting features which appear in DNS. Notably, in the unstable regime, depending on the symmetry of the initial perturbation, sheets and jets emerge from the cylindrical filament suffering further instabilities and leading to fragmentation, pinch-off etc.

• Faraday waves on a cylindrical filament using a radial forcing can have potentially many applications e.g. to stabilise the RP instability in falling jets or as a mechanism to atomize fuel jets into droplets inside internal combustion engines. The second is an interfacial solver which implements the one-fluid algorithm (labelled as OFS: one fluid solver) for simulating two-phase flows developed by the author and is described earlier in chapter 2. The principal difference between the two codes is that the FSS solves only for the heavier phase modelling surface tension as a boundary condition. In contrast the OFS necessarily solves for both phases implementing surface tension as a body force. For high density ratio, it is expected that both codes provide very similar results. The benefit of the FSS approach is that it can save substantial computing time by not solving the lighter phase. One of the purposes of this study is to establish the lower limit of density ratio where substantial differences between OFS and FSS maybe seen and consequently the saving in computational time due to the FSS approach, is offset by the accompanying loss in accuracy of results. We discuss four test cases involving inviscid, capillary waves with analytical solutions in the linearised limit. Access to analytical results also aid in benchmarking both codes.

Free perturbations

We consider linearised analytical solutions to inviscid, irrotational, free, capillary oscillations occurring on quiescent base states of two different geometries. The analytical solution for the planar case is available in the literature Debnath (1994) and the analytical solution for the circular drop has not been reported earlier and is presented in Singh et al. (2019b).

Cauchy-Poisson problem -Cartesian

As shown in figure 6.1a, we consider two horizontally unbounded fluids (i.e.

•  x  •), inviscid, immiscible fluids of density r I and r O (r O <r I ) and surface tension coefficient T, at rest forming a pool of depth H 1 and H 2 . We restrict the present analysis to capillary waves, neglecting gravity. Superscripts I and O are used to distinguish between inner (lower) and outer (upper) fluid respectively. As shown in figure 6.1a, the difference between the interface in the base and the perturbed state is given by h(x, t). A localised interfacial perturbation is introduced initially of the form h(x, 0) assumed to be of the form of a modified Gaussian [START_REF] Miles | The cauchy-poisson problem for a viscous liquid[END_REF]. Due to the restoring force of surface tension, the interface returns to its undisturbed equilibrium position, overshooting it due

⇣ (x L/2) 2 ⌘ ⇣ 1 2 (x L/2) 2 ⌘ b) h(q, 0) = 0.1 exp ⇥ q 2 ⇤ ⇥ 1 2q 2 ⇤
to inertia producing oscillations as a result. Interfacial travelling waves are created which propagate horizontally to the left and right. In

h(x, t) = 1 p 2p R • • dk exp[ikx] h0 (k) cos [w 0 (k)t] , (6.1) f I (x, y, t) = w 0 p 2p R • k= • dk exp[ikx]G I (k, y) sin [w 0 (k)t] , G I (k, y) ⌘ exp[|k|y] + exp[ |k|y 2|k|H 1 ] 1 exp( 2|k|H 1 ) h0 (k) |k| , (6.2) f O (x, y, t) = w 0 p 2p R • k= • dk exp[ikx]G O (k, y) sin [w 0 (k)t] , G O (k, y) ⌘ exp[ 2|k|H 2 + |k|y] + exp[ |k|y] 1 exp( 2|k|H 2 ) h0 (k) |k| . (6.3)
where w 2 0 (k) ⌘

T|k| 3 r I coth(|k|H 1 )+r O coth(|k|H 2 )
is the dispersion relation for free oscillations of a single Fourier mode of wavenumber k. For general h0 (k), the Fourier integrals in expressions 6. 1-6.3 are not analytically tractable and need to be numerically computed. For the special initial condition when the interface is perturbed by a single Fourier mode i.e. h(x, 0) = a 0 cos(k 0 x) and its Fourier transform, h0 (k) = p p/2 [d(k k 0 ) + d(k + k 0 )], expressions 6.1-6.3 simplify to h(x, t) = a 0 cos(k 0 x) cos [w 0 (k 0 )t] (6.4)

f I (x, y, t) = a 0 w 0 k 0 cos(k 0 x)L I (k, y) sin [w 0 (k 0 )t] ,
(6.5)

f O (x, y, t) = a 0 w 0 k 0 cos(k 0 x)L O (k, y) sin [w 0 (k 0 )t] .
(6.6)

L I (k, y) ⌘ exp[k 0 y] + exp[ k 0 y 2k 0 H 1 ] 1 exp[ 2k 0 H 1 ] L O (k, y) ⌘ exp [k 0 y 2k 0 H 2 ] + exp [ k 0 y] 1 exp [ 2k 0 H 2 ]
with w 0 (k 0 ) in these equations is obtained by replacing k with k 0 in the dispersion relation. We test equation 6.1-6.3 and 6.4-6.6 using data obtained from the OFS and the FSS in section 6.1. Tests for two initial conditions are reported viz. a localised perturbation of the form h(x, 0) = a 0 exp 1968) and a single Fourier mode i.e. h(x, 0) = a 0 cos(k 0 x).

⇣ (x x 0 ) 2 ⌘ ⇣ 1 2 (x x 0 ) 2 ⌘ Miles (

Cauchy-Poisson problem -circular

The Cauchy-Poisson problem discussed in the earlier sub-section, can be generalised to a two-dimensional circular base state as shown in figure 6 1 T/ r I + r O R 3 0 taking into account the (linearised) inertia of the outer fluid. We extend the analysis of Bohr (1909) to provide an analytical solution for a localised interfacial perturbation at time t = 0. This solution is the circular equivalent of the Cauchy-Poisson problem discussed in the earlier section. This analysis is new and has not been presented in the literature before. As seen in figure 6.1b, the base state consists of immiscible, quiescent fluids inside and outside a circular drop of radius R 0 . The fluid inside has a density r I and the one outside has r O , assumed to be unbounded. The surface tension coefficient is T and in the base state, the pressure inside and outside are uniform differing by an amount P I b P O b = T/R 0 (subscript b for base state variables). At time t = 0, an interfacial perturbation of the form h(q) is imposed on the interface and leads to travelling waves on the interface, due to surface tension. With the Fourier transform of the initial perturbation being defined as, f (m) ⌘ F (h(q, 0)), it is shown in Singh et al. (2019b) that the velocity potential and the shape of the interface at later time t is given by,

h(q, t) = 1 p 2p R • m= • dm f (m) exp[imq] cos [w 0 (m)t] , (6.7) f I (r, q, t) = w 0 R 0 p 2p R • m= • dm exp[imq] ✓ r R 0 ◆ |m| f (m) |m| sin [w 0 (m)t] , (6.8) f O (r, q, t) = w 0 R 0 p 2p R • m= • dm exp[imq] ✓ r R 0 ◆ |m| f (m) |m| sin [w 0 (m)t] ,
(6.9)

where w 2 0 (m) ⌘ T(|m| 3 |m|) (Bohr, 1909). For general f (m), the integrals in 6.7-6.9 need to be solved numerically. However, for the special initial condition when the interface is perturbed initially, by a single Fourier mode i.e. h(q, 0) = a 0 cos(m 0 q) and f (m

R 3 0 (r O +r I )
) = p p/2 [d(m m 0 ) + d(m + m 0 )],
the integral in these expressions simplify leading to, h(q, t) = a 0 cos(m 0 q) cos [w 0 (m 0 )t] , (6.10) .12) with w 0 obtained by replacing m with m 0 in the dispersion relation provided below equation 6.9. We test equations 6.7-6.10 against data obtained from the OFS and FSS for the localised initial condition h(q, 0) = a 0 exp ⇥ q 2 ⇤ ⇥ 1 2q 2 ⇤ as well as a single Fourier mode h(q, 0) = a 0 cos(m 0 q) in section 6.1. Note that the modification to the Gaussian in h(q, 0) is chosen such that upto linear order i.e. O(h), the area of the perturbed droplet is the same as the circular droplet i.e.

f I (r, q, t) = a 0 w 0 R 0 m 0 cos(m 0 q) ⇣ r R 0 ⌘ m 0 sin [w 0 (m 0 )t] , (6.11) f O (r, q, t) = a 0 w 0 R 0 m 0 cos(m 0 q) ⇣ r R 0 ⌘ m 0 sin [w 0 (m 0 )t] . ( 6
R 2p 0 [R 0 + h(q, 0)] 2 dq = pR 2 0 + O(h 2 ).

Parametrically forced perturbations

In the last section, we considered free oscillations occurring under the restoring force of surface tension. In this section, we obtain linearised analytical solutions for inviscid, irrotational, parametrically forced, capillary waves occurring on quies- cent, time-periodic base states. The two geometries considered are the same as before with an additional time-periodic body force now acting in both the fluids.

Faraday waves -Cartesian

Faraday waves Faraday (1831b) arise when a container partially filled with fluid is vertically vibrated at a frequency W and amplitude h Benjamin and Ursell (1954). In a frame of reference attached to the vibrating container, the effective acceleration appears as a time-periodic body force, acting normal to the unperturbed interface. For any driving frequency, above a certain critical forcing amplitude, standing waves appear at the interface which may destabilise leading to atomisation. As shown in figure 6.2a, the base state comprises a flat horizontally unbounded interface (dashed lines) subject to a time-periodic body force of the form h cos(Wt) (shown using arrows) which produces a time-periodic pressure field with zero velocity everywhere. An interfacial Fourier mode in the form of a standing wave of wavenumber k is imposed and we seek the amplitude and Tuckerman (1994). In the next section, we obtain the stability chart of equation 6.13 using numerical Floquet analysis and also solve it numerically. We compare its predictions with the FSS and OFS results, for the case r I > r O .

a(t)

Faraday waves -circular

Faraday waves in a Cartesian geometry described in the previous sub-section, can be generalised to a circular drop analogous to the way in which the planar Cauchy-Poisson problem was generalised to the circular case. Shown in figure 6.2b is the cross-section of a cylindrical filament (r I > r O ) being subjected to a time periodic body force in the radial direction i.e. normal to the unperturbed filament. The time-periodic body force is chosen analogous to the Cartesian case i.e. acting in a direction normal to the unperturbed interface. In order to prevent multi-valuedness at the origin (i.e. as r ! 0), we choose the body force to be spatially dependent and of the form h cos(Wt) (r/R 0 ) (see Adou and Tuckerman (2016) for Faraday waves on a spherical droplet where a similar forcing is used). The base state comprises of a circular droplet (in two dimensions) subject to a radial, time-periodic body force, applied in both the fluids. Due to surface tension, the pressure jump condition needs to be satisfied at the interface at all times, given by p .

I b (R 0 , t) p O b (R 0 , t) = T/
a(t) + m R 0 " T(m 2 1) (r I + r O )R 2 0 + ✓ r I r O r I + r O ◆ h cos(Wt) # a(t) = 0. (6.16)
For zero forcing, equation 6.16 reduces to a simple-harmonic oscillator equation.

A stability chart of equation 6.16 can be obtained using similar techniques as equation 6.13. Knowing a(t) from equation 6.16, the velocity potential in both fluids are obtained from

f I (r, q, t) = m 1 R 0 cos(mq) ✓ r R 0 ◆ m . a(t), (6.17) f O (r, q, t) = m 1 R 0 cos(mq) ✓ r R 0 ◆ m . a(t). (6.18)
The stability chart of equation 6.16 predicts that oscillatory Fourier modes (m = 2, 3, 4 . . .) which are stable in the absence of forcing (Rayleigh, 1879a) can be destabilised through sufficiently high forcing h. We compare DNS results obtained from the FSS and the OFS against these analytical predictions in the next section.

Results

Having presented four sets of analytical results in the previous section, we compare predictions from these to DNS results obtained from the OFS and the FSS. All simulations conducted as a part of this study, are summarised in table 6 et al., 2011) are also seen to be small, as is validated from the Weber number

We ⌘ r I U 2 max R 0 /T ⇡ 10 6 where U max is the maximum velocity found in the domain at t = 20 s. Note that t = 20 s represents a long time, as for the given choice of parameters, the time period of the azimuthal mode (m = 2) of oscillation is ⇡ 2.56 s. Table 6.2 shows the We for FSS, OFS and Basilisk.

Free perturbations

In this section, we compare predictions made earlier in equations 6.1-6.4 (Cartesian geometry) with DNS results.

Cartesian

In figures 6.4a-6.4d we compare analytical predictions obtained from equation 6.4 to data obtained from FSS and OFS. A cosine wave of amplitude a 0 and wavenumber k is initiated (cases 2-5 in table 6.1) as shown in the simulation geometry in figure 6.2a (without the forcing shown by arrows) and the amplitude of the wave is tracked as a function of time at x = 0 (any other location will also produce similar answers). The nonlinearity parameter e ⌘ a 0 k is a measure of the magnitude of nonlinear effects and is chosen to be small (approximately 0.063) in order to facilitate comparison with the linearised prediction in equation 6.4. As seen in figures 6.4a-6.4c, for density ratios (r r ⌘ r I /r O ) ranging from 10 3 -50, the predictions from the FSS agree reasonably well with the analytical solution, although numerical damping manifests itself, from the fourth oscillation onwards. At a density ratio of 10 (figure 6.4d), the FSS results show substantial deviation after the first oscillation. For further validation of these results, we have superimposed results obtained from Basilisk Popinet (2014), in this figure. It is seen that starting at a density ratio of 10 and below, inertial effects of the upper fluid are strongly felt on the oscillations and cannot be neglected at a first approximation. A qualitatively similar trend is seen in the horizontal velocity profiles in figures 6.5a and 6.5b which are plotted at different times within one inviscid, oscillation time-period T 0 ⌘ 2p/w 0 for two density ratios. For r r = 10 3 , the horizontal velocity profiles obtained from the FSS, the OFS and equation 6.5 are indistinguishable from each other as seen in figure 6.5a. This is not so at r r = 10 in figure 6.5b where the accuracy of the FSS is seen to be much lower than that of the OFS. In order to compare the effect of free perturbations further, we consider a localised interfacial perturbation of the form shown in figure 6.1a. As discussed earlier in section 6.0.1, this is referred to as the Cauchy-Poisson problem Debnath (1994). The initial perturbation is chosen to be a modified Gaussian of the (in collaboration with Manpreet Singh, M.Tech Thesis, IIT Bombay, 2018) [START_REF] Miles | The cauchy-poisson problem for a viscous liquid[END_REF]). In order to compare with DNS results, we evaluate the numerical Fourier transform of h(x, 0), further numerically evaluating the integral in expression 6.1. This is plotted as a function of time in figures 6.6a-6.6l where a good agreement with DNS results obtained from the FSS is apparent. Note that the linearised theory assumes a horizontally unbounded domain, while the simulation uses a finite domain (L ⇥ L). Consequently, the lack of agreement towards the end of the simulation in figure 6.6k is expected due to reflection at the domain boundaries. Such finite domain effects are not accounted for in the analytical formulae 6.1.

form h(x, 0) = a 0 exp h (x x 0 ) 2 i h 1 2 (x x 0 ) 2 i (such that it is an volume conserving perturbation i.e. R • • dx h(x, 0) = 0, see

Circular

We next simulate a single Fourier mode on the circular drop shown in figure 6.2b (without the forcing shown using arrows). Figure 6.7a, shows the time signal obtained from imposing a single azimuthal Fourier mode on the circular drop. The amplitude of the interface is tracked at q = p/2 radians, where q is measured anticlockwise with respect to the positive direction of the x axis. In this geometry, the nonlinearity parameter e az ⌘ a 0 m/R 0 , is also chosen to be small (= 0.1) to compare with the linearised prediction in equation 6.10. Here too we observe that from r r = 10 3 50, the agreement with analytical solution as well as OFS is good. However starting at density ratio r r = 10, the agreement between the FSS starts showing substantial deviations due to inertial effects of the outer fluid (r O ) being no longer small. As a validation of our code in this geometry, we have also added the time signal obtained from Basilisk Popinet (2014). Figure 6.8a and 6.8b show the radial velocity profiles within the droplet plotted as a function of radial distance. It is seen that at r r = 10, the FSS shows substantial disagreement with theory and OFS. Note that the agreement between theory and OFS for this case, is less than that seen in the earlier geometry in figure 6.5b, due to the slightly larger value of the nonlinearity parameter here. For further comparison with DNS results, we generalise the perturbation from a single azimuthal Fourier mode, to a localised perturbation of the form h(q, 0) = a 0 exp ⇥ q 2 ⇤ ⇥ 1 2q 2 ⇤ , cf. figure 6.1b. The response of the interface to this perturbation is shown in figures 6.9a-6.9i where we can see a good agreement with analytical predictions from equation 6.7. The Fourier integral in equation 6.7 is solved numerically in MATLAB (2015). Note that unlike k in the Fourier integrals in 6.1, m is restricted to be an integer. Consequently, dm in 6.7 and 6.8 is also taken to be unity while evaluating the integrals numerically.

In the FSS, the effect of inertia for the outer fluid is completely neglected, while that from the inner fluid is taken into account. In order to estimate the contribution from the nonlinear inertia of the outer fluid, we obtain the time signal for free oscillations of a single Fourier mode in the Cartesian and circular geometry, with higher values of e or e az . Figure 6.10a and 6.10b show that at high density ratio and increased e (or e az ), the frequency of oscillation shows a nonlinear (finite amplitude correction to the dispersion relation) correction manifested as a drift in the time signal, when compared to the linearised predictions from equation 6.4 or 6.10. Interestingly, the OFS shows higher numerical dissipation as compared to the FSS. As seen in figure 6.10b, this is also true for data obtained from Basilisk Popinet (2014), although the latter has lower numerical dissipation than our in-house developed OFS.

Parametrically forced perturbations

In this section, we discuss the case of forced perturbations. Analogous to the earlier section, we present results in two geometries viz. Cartesian in figure 6.2a and circular in figure 6.2b. As indicated by the arrows, a time-periodic body force acts normal to the unperturbed interface in both the fluids.

Cartesian

For the Cartesian case, shown in figure 6.2a, the base state consists of a flat interface, zero velocity in both the fluids and a time-periodic pressure field. In order to study the linear stability of this base state, an interfacial perturbation in the form of a standing Fourier mode of amplitude a(t) is imposed at t = 0. Equation 6.13 governs the amplitude a(t) of this standing wave of wavenumber k, for a given forcing frequency W. Using numerical Floquet analysis, we obtain the stability chart for equation 6.13 (see Singh et al. (2019b) for details, Floquet analysis was done by the M.Tech student and detail steps for obtaining these charts are available at Singh et al. (2019a)). Shown in figure 6.11a, are the stability tongues of equation 6.13. On the forcing-amplitude h versus wavenumber space k, inside the tongue-like regions labelled as U (unstable), equation 6.13 shows exponential growth and instability. In the regions labelled S (stable), the response is bounded and oscillatory. This is validated in figure 6.11b and 6.11c which are chosen to be linearly stable and unstable respectively (cf. figure 6.15b). In the stable case (case 14, figure 6.11b), the FSS shows very good agreement with numerical solution of the Mathieu equation 6.13. In the unstable case (case 15, figure 6.11c) some deviation from linear theory is visible at later times, primarily due to nonlinear effects associated with growing amplitude due (in collaboration with Manpreet Singh, M.Tech Thesis, IIT Bombay, 2018) to instability, see inset of figure 6.11c. The effect of reduction of density ratio can be seen in figure 6.11d (see stability chart for this, figure 6.15a) where it is observed that the FSS display more deviations when compared to the OFS.

Circular

Faraday waves seen in Cartesian geometry can also arise on a circular drop as seen in figure 6.2b, when the drop is subjected to a radial time-periodic body force, normal to the unperturbed interface in both fluids. The base state comprises of quiescent fluids, with a circular interface and a time-periodic body force acting radially in both fluids. To study its linear stability, we impose an azimuthal Fourier mode of amplitude a(t) on the interface at t = 0. Note that in our DNS in both FSS and OFS, only a quarter droplet is simulated as shown in figure 6.2c. The boundary conditions for the simulations are provided in table 6.1. The equation governing a(t) is the Mathieu equation 6.16 (see chapter 5 for a study of three-dimensional perturbations on a cylindrical filament). An important distinction between equation 6.13 and 6.16 is that, for the former, due to the unboundedness assumption in the horizontal x direction, any value of k (0 < k < •) in equation 6.13 is permitted. However for the azimuthal Fourier mode, due to periodicity requirements, m in equation 6.16 must be a positive integer. Consequently, in the stability chart for equation 6.16 obtained using numerical Floquet analysis, only those values of forcing strength are meaningful, which correspond to integer values of m (shown using symbols and dotted lines in figure 6.12a). Figure 6.12b and 6.12c show two cases, which correspond to stable and unstable Faraday waves on the circular droplet (cases 18 and 19 in stability chart in figure 6.15d). In the stable case, very good agreement with the FSS is seen. Note that the OFS shows some numerical dissipation for this case. However for the unstable case in figure 6.12b, substantial disagreement is visible with the FSS after a few oscillations, primarily due to nonlinear effects. It is important to note that such disagreement is to be expected in the unstable case, as our theoretical prediction through the Mathieu equation 6.16 is valid in the linear regime only. The effect of density ratio on Faraday waves on a droplet is seen in figure 6.13a and 6.13b (for figure 6.13b, the stability chart is presented in the figure 6.15c). Like the Cartesian case, it is seen that at high values of density ratio viz. r r = 10 3 and small e az , the agreement between analytical predictions from equation 6.16 and FSS/OFS is good, although nonlinear effects are seen eventually. This is not so at lower values of r r (figure 6.13b) where the FSS and the OFS show substantial disagreement with each other. Note that for r r = 10,

= CFL ⇥ min ✓ Dx 2U max , q (r I +r O )(Dx) 3 4pT ◆ 0 1 2

Conclusions

In this chapter, we have presented four analytical test cases for free and forced, capillary, interfacial oscillations. These test cases cover stable as well as unstable oscillations. Explicit analytical formula for the time evolution of the interface as well as velocity potential in both fluids, in the linearised limit, have been obtained in all cases. These test cases are valuable additions to the repository of analytically known results for capillary oscillations. They are expected to be of utility, especially for benchmarking numerical models for implementing surface-tension (see a recent extensive review by Popinet (2018), section 6 where the importance of analytical solutions as test cases for validating numerical surface tension models, is discussed). Using these test cases as reference, we have compared results obtained from two in-house developed solvers viz. the FSS (free-surface solver) and the OFS (one-fluid solver), in two base-state geometries -Cartesian and circular (two-dimensional). In all studies presented here, we find that the FSS produces good agreement with analytical predictions, when the density ratio is greater than 50 and the nonlinearity parameter e, e az < 0.1. At larger values of e (or e az ) and high density ratio (⇡ 10 3 ), the FSS is found to have lower numerical dissipation than the OFS, leading to better agreement with linearised theory. For the problems studied here, we conclude that when r r > 50, it is computationally efficient to use an FSS rather than an OFS. At density ratios of 10 and lower, inertial effects of the second fluid are substantial and large errors result, when they are neglected.

Case FSS OFS F=OFS/FSS Machine -Intel ® Xeon(R) CPU E5-1650 v3 @ 3.50GHz

These empirical estimates provide useful guidelines while deciding, whether to use an OFS or an FSS based approach when solving a two phase flow problem.

The FSS is expected to be significantly faster compared to the OFS, especially when using a very refined grid. The net savings in computational time when comparing the FSS and OFS, is determined as a result of competition between having to solve in fewer cells (in case of FSS) but with the additional overhead of an extra pressure-poisson equation (compared to the OFS). In typical simulations presented in this study, we found that the FSS is always faster compared to OFS by a factor which varied from 1.5 30 , depending on the case being considered. Table 6.5 shows some typical comparisons for the computational time between the OFS and the FSS.

the filament. We show that the amplitude of parametrically forced, capillary, standing waves is governed by the Mathieu equation. The results are then verified by three dimensional direct numerical simulations using the open source solver Gerris. [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF].

• In Chapter 6, we propose four test cases for studying the effect of density ratio on inviscid, capillary waves. Analytical results are compared against results obtained from two in-house solvers for solving the Euler equation with surface tension and good agreement is reported.

As we reach the end of this thesis, it is useful to reflect on future research directions. The two problems studied in Chapters 3 and 4 are intimately related, suggesting that the continuous spectrum is present whenever the domain is infinite. This needs a rigorous analysis currently underway. More importantly, we need to reformulate both studies in Chapters 3 and 4 with non-zero initial distributions of vorticity (analogous to, but more general than the study by

Prosperetti and Cortelezzi (1982)).

For vortical perturbations at large distances, one intuitively expects that such initial conditions will have large projections on continuous spectrum eigenmodes when compared to the discrete spectrum eigenmodes. This is a particularly important step towards formulating weakly nonlinear theories for understanding interactions between vorticity field and interfacial waves. Such an analysis is important for understanding phenomena like jetting observed from large amplitude axisymmetric capillary-gravity waves, as reported in Chapter 3.

The Faraday wave problem analysed in chapter 5 assumes potential flow. The observed stabilisation of hitherto unstable Rayleigh-Plateau modes is rather short-lived and introduction of viscosity is expected to prolong the stabilisation.

A viscous theoretical formulation of this problem is currently underway. It is known from literature that in the presence of viscosity, the amplitude is not governed by a damped Mathieu equation (Adou and Tuckerman, 2016; Kumar and Tuckerman, 1994). Insights into why this is so is necessary and is proposed as future work.

The in-house code which solves the Navier-Stokes equation with an interface between two fluids has been a substantial effort to implement. In particular, including surface tension forces and calculating curvature to obtain satisfactory results from a large number of benchmark tests, while ensuring that the code was reasonably fast. We have the following algebra

r L x + r U z = 2(r 2 L k + Ir U r L p x) r L + r U + 2(r 2 U k + Ir L r U p x) r L + r U (I.8) = 2k(r 2 L + r 2 U ) r L + r U + 4Ir L r U p x r L + r U (I.9) 1 2k 1 l + k = l k 2k(l + k) = (l k) 2 2k(l 2 k 2 ) = (I p x k) 2 2k(x + k 2 ) (I.10) x + 2kI p x k 2 2k(x + k 2 ) = x k 2 2k(x + k 2 ) + I p x (x + k 2 ) (I.11) (r L x + r U z) ✓ 1 2k 1 l + k ◆ = r 2 L + r 2 U r L + r U ! ✓ x k 2 x + k 2 ◆ 4r L r U r L + r U x (x + k 2 ) +2 p x k 2 (r 2 L + r 2 U ) + (x k 2 )r L r U k(x + k 2 )(r L + r U ) ! I = (r 2 L + r 2 U )(x k 2 ) 4xr L r U (r L + r U )(x + k 2 ) + 2 p x k 2 (r 2 L + r 2 U ) + (x k 2 )r L r U k(x + k 2 )(r L + r U ) ! I (I.12) L = 2nk 2 r L + r U  ⇣ r L + r U ⌘ + 1 2k ⇣ r L x + r U z ⌘ ✓ r L x + r U z l + k ◆ = 2nk 2 + 2nk 2 (r 2 L + r 2 U )(x k 2 ) 4xr L r U (r L + r U ) 2 (x + k 2 ) + 4nk 2 p x k 2 (r 2 L + r 2 U ) + (x k 2 )r L r U k(x + k 2 )(r L + r U ) 2 ! I (I.13) s 2 + sL(s) + w 2 0 = n 2 (x + k 2 ) 2 2n 2 k 2 (x + k 2 ) 2n 2 k 2 (r 2 L + r 2 U )(x k 2 ) 4xr L r U (r L + r U ) 2 + w 2 0 4n 2 k 2 p x k 2 (r 2 L + r 2 U ) + (x k 2 )r L r U k(r L + r U ) 2 ! I = A I B (I.14) 1 s 2 + sL + w 2 = A + IB A 2 + B 2 1 s 2 + sL + w 2 x(s) = 2(r L k + Ir U p x) (r L + r U )(A 2 + B 2 ) (A + BI) = 2 (r L + r U )(A 2 + B 2 ) (r L kA + r L kBI + Ir U p xA r U p xB) = 2 (r L + r U )(A 2 + B 2 ) (r L kA

I.2 Faraday waves on a cylindrical filament, Chapter 5

The derivation is provided here for Faraday waves on a cylindrical fluid filament. fI (r, z, q, t) = 0 + f I (r, z, q, t), pI (r, z, q, t) = P I b (r, t) + p I (r, z, q, t), fO (r, z, q, t) = 0 + f O (r, z, q, t), pO (r, z, q, t) = P O b (r, t) + p O (r, z, q, t),

The equation of filament in the base state is r = R 0 and thus x 0 ? = R 0 . The base state momentum equations in the (normal) radial direction are (base state velocity is zero in both fluids) In the presence of an interfacial perturbation, the difference of total pressure (base + perturbation) on either sides of the interface equals surface tension coefficient T times the surface divergence of the unit normal at the perturbed interface [START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF], P I b (R 0 + h, t) + p I (q, z, R 0 , t) P O b (R 0 + h, t) p O (q, z, R 0 , t) = T (r • q) r=R 0 +h , (I.28)

(
where q is the unit normal to the perturbed cylinder and the perturbation pressure p is evaluated at the unperturbed cylinder radius r = R 0 in order to retain terms only upto O(a(0)) only. In order to compute curvature at the perturbed interface we define (Bush, 2013)

Q(q, z, r, t) ⌘ r R 0 h(q, z, t) = r R 0 a(t) cos(mq) cos(kz) (I.29) At a linear approximation, the divergence of q equals the Laplacian of Q (see equations 5. O(a(0)). An expression for the difference of base pressures on the left hand side of equation I.28 can be obtained using equation I.18. After some simple algebraic manipulations, this is

P I b (R 0 + h, t) P O b (R 0 + h, t) = T R 0 ⇣ r I r O ⌘ Z R 0 R 0 +h
G(r 0 , t)dr 0 (I.31)

where we have used the base state pressure-jump condition to replace P I b (R 0 , t) P O b (R 0 , t) with T/R 0 (see below equation I.18). Defining G(r, t) ⌘ ∂Y(r,t) ∂r and using Taylor series expansion, we find

Z R 0 R 0 +h G(r 0 , t)dr 0 = Y(R 0 , t) Y(R 0 + h, t) = ✓ ∂Y ∂r ◆ R 0 h(q, z, t) + O(h 2 )
⇡ G(R 0 , t)h(q, z, , t).

⇠ ⇣ k2 + 1 ⌘ 1/2 k ⌥ 1 m k ⇣ k2 + 1 ⌘ 1/2 ⌥ 1 m ln  k + ⇣ k2 + 1 ⌘ 1/2 ± 1 m k ⇣ k2 + 1 ⌘ 1/2 n k + ⇣ k2 + 1 ⌘ 1/2 o ⇠ ⇣ k2 + 1 ⌘ 1/2 k ln  k + ⇣ k2 + 1 ⌘ 1/2 ⌥ ⇣ k2 + 1 ⌘ 1/

I.4 Inviscid limit of dispersion relation for hollow gas filament, chapter 4

We show here that the dispersion relation 4.38 has the correct inviscid limit and for small value of viscosity. Re-writing 4.38 in dimensional variables and

R 4 0 s 2 n 2 + ( 2m(m + 1) 0 @ 1 R2 0 2m K m+1 ( R0 ) (m 1)K m 1 ( R0 ) R0 K 0 m 1 ( R0 ) 1 A ) R 2 0 s n + TR 0 r µ 2 m ⇣ m 2 1 ⌘ 1 (m + 1) R0 K m 1 ( R0 ) K 0 m 1 ( R0 ) ! = 0 (I.50)
which may be written as

I.5 Expression of vorticity for hollow gas filament, Chapter 4

The expression for C(s) in equation 4.63 is obtained using Cramers rule on equations 4.66-4.68. This is The series representation of K m (z) for m = 0, 1, 2 and 3 are obtained below from 4.70. For m > 0, it is seen that z = 0 is a pole and a logarithmic branch point. The integrals 2 and 6 maybe shown to be zero as R ! •. We evaluate the integral over the little circle 4. This may be written using s = d exp(Iq) as, Due to a jump in the integrand across the branch cut, the integrals 3 and 5 do not cancel each other. With s = x exp (Ip) and s = x exp ( Ip) respectively, I 3 and I 5 become

C(s) = s a 0 m R 0 0 0 2m(m + 1) R 2 0 T rR 2 0 (1 m 2 ) 0 s 1 + 2m(m + 1) h 2 R 2 0 ! s K m 1 (hR 0 ) 2h m R 0 0 K 0 m 1 (hR 0 ) 2m(m + 1) R 2
K 0 (z) = h g + ln ⇣ z 2 ⌘i + O ✓ ⇣ z 2 ⌘ 2 ◆ K 1 (z) = 1 z + 1 4 h 2 ln ⇣ z 2 ⌘ 1 + 2g i z + O ✓ ⇣ z 2 ⌘ 3 ◆ (I.
I 3 = 1 2pI Z • 0 L(x exp (Ip)) Q(x exp (Ip)) exp [ xt] dx, I 5 = 1 2pI Z • 0 L(x exp ( Ip)) Q(x exp ( Ip))
exp

Proof:

We have the identities pIJ m (z) = exp ( mpI/2) K m (z exp ( pI/2)) exp (mpI/2) K m (z exp (pI/2)) (I.66) pY m (z) = exp ( mpI/2) K m (z exp ( pI/2)) + exp (mpI/2) K m (z exp (pI/2))

2 + k 2 a 2 1 ⌘

 21 of Congress Prints and Photographs Division Washington, D.C. 20540 USA. Reproduction numbers -LC-USZ62-120355, LC-USZ62-120893. (Reproduced under fair use policy for educational use) of St. Andrews archives, Lord Rayleigh -Biographical. NobelPrize.org. ©Nobel Media AB 2018. accessed: 27 May 2019 (Reproduced with special permission for educational use) where p is angular frequency of oscillation as in equation 26 in "On the Capillary Phenomena of Jets" (Rayleigh, 1879a).

( a )

 a Lord Kelvin (1824[START_REF] Pederson | On the surface-tension of liquids investigated by the method of jet vibration[END_REF] Source: pBarron (2007)(Reproduced under fair use policy for educational use) (b) Perturbed spherical interface

c( a )

 a Diagram of the accelerated mass of liquid. The arrow shows the direction of the acceleration f cos wt, and the curve C is the boundary of the free surface.

  Stability chart for the solutions of Mathieu's equation d 2 a dT 2 + (p 2q cos 2T)a = 0.

Figure 1 . 1 :

 11 Figure 1.1: Reproduced from Benjamin and Ursell (1954) under fair use policy for educational use.

Figure 1 . 2 :

 12 Figure 1.2: Geometrical configurations of waves studied in the present work

3 .

 3 Lagrangian methods in which the grid follows the fluid. e.g. Droplet breakup by Oran and Boris (2005)

Figure 2 . 1 :

 21 Figure 2.1: Reconstruction using volume fractions.

Figure 2 . 2 :

 22 Figure 2.2: Normal form of a line in the cell

Figure 2 . 3 :

 23 Figure 2.3: Extrapolation of the line drawn at the central cell into a 3 ⇥ 3 stencil.

( a )

 a Initial Condition for Volume Fraction field (b) After advecting Volume Fraction for 20 time steps (c) After Volume Fraction density for 20 time steps(d) After Volume Fraction density for 50 time steps

Figure 2 . 4 :Figure 2 . 5 :Figure 2 . 6 :

 242526 Figure 2.4: Advection of Volume Fraction using finite difference method

  possibilities arise depending on whether the Heaviside functions is zero or unity. This gives rise to different shapes under the line. The four broad possibilities are 1. 0  a  0.5 and 0  m  0.5 2. 0.5 < a  1.0 and 0  m  0.5 3. 0  a  0.5 and 0.5 < m  1.0 4. 0.5 < a  1.0 and 0.5 < m  1.0 Each of the above four cases give rise to two sub-cases which are dealt with below.

Figure 2 . 7 :

 27 Figure 2.7: Advection test for velocity field u=1,v=0

Figure 2 . 8 :Figure 2 . 9 :Figure 2 . 10 : 2 . 3 . 5 4. 5 .

 28292102355 Figure 2.8: Advection test for velocity field u=2,v=1

  Initial condition for shear test (b) After advecting 250 steps (c) After advecting 500 steps (d) After advecting 1000 steps

Figure 2 . 11 :Figure 2 . 12 :

 211212 Figure 2.11: Advection test result for shear velocity field

Figure 2 . 14 :

 214 Figure 2.14: Staggered grid for velocity and pressure

Figure 2 . 15 :

 215 Figure 2.15: Three sub stencils for reconstruction of u c+1/ 2

Figure 2 . 17 :

 217 Figure 2.17: The matrix A

Figure 2 . 18 :

 218 Figure 2.18: Profile after 10 Gauss-Seidel iterations, Blue -Initial guess, Orange -After 10 iterations

Figure 2 . 19 :

 219 Figure 2.19: Eigenvalues l with respect to k N

Figure 2 .

 2 Figure 2.20: Multigrid Algorithm

  Error after reaching max |r| < 10 3

Figure 2 . 22 :Figure 2 . 23 :

 222223 Figure 2.22: Solver results of for N = 512

Figure 2 . 24 :

 224 Figure 2.24: Comparison of SOR and Multigrid Time of execution

Figure 2 . 25 :

 225 Figure 2.25: Lid driven cavity

Figure

  Figure 2.26: Square cavity with regularised boundary conditions

Figure 2 . 27 :Figure 2 . 28 :

 227228 Figure 2.27: Regularised lid driven cavity

Figure 2 . 29 :

 229 Figure 2.29: Planar wave viscous oscillations

Figure 2 . 30 :

 230 Figure 2.30: Planar viscous Rayleigh-Taylor instability

1 . 2 . 3 .Figure 2 . 32 :

 123232 Figure 2.31: Axisymmetric viscous Rayleigh-Taylor instability

Figure 2 . 33 :

 233 Figure 2.33: Comparison with Bhaga and Weber (1981), Red:in-house code

  Waves at the interface of two immiscible viscous fluids have been of interest to the fluid dynamics community for technological applications like drug delivery (Tsai et al., 2014) as well as from a fundamental perspective in understanding nonlinear systems and phenomena like atomisation (Li and Umemura, 2014; Puthenveettil and Hopfinger, 2009). Dating from the seminal study by Michael Faraday (Faraday, 1831a), forced standing interfacial oscillations have received continued interest and attention in fluid dynamics literature having become a canonical system for studying nonlinear physics and pattern formation in parametrically forced systems (Perinet et al., 2009; Rajchenbach and Clamond, 2015). A study of such oscillations, especially the ones occurring in cylindrical geometries have significant practical utilities. They appear during sloshing in partially liquid-filled cylindrical containers subject to vertical vibration (Henderson and Miles, 1990). Atomisation as a result of instabilities of standing waves appearing on a vibrated sessile droplet, have been proposed as a technique for spray cooling of microprocessors (James et al., 2003). Axisymmetric free interfacial oscillations similarly have many applications. At large scales where surface tension and viscous effects are negligible, they are relevant for understanding seiches in harbours (Mack, 1962). The first solution to an IVP for surface waves, was formulated by Cauchy and Poisson (Cauchy, 1827; Poisson, 1818) in response to a prize problem an-61 nounced by the French Academy of Sciences (Craik, 2004). They independently solved the IVP for inviscid, irrotational, axisymmetric as well as two-dimensional planar geometry. Their analysis comprised of expressing the shape of the interface as an integral superposition of Fourier/Bessel modes (Craik, 2004). Many real-world problems concerning surface waves involve viscous dissipation and boundary layer formation at the interface. The viscous normal mode analysis by Lamb (1932) for planar waves obtained the dispersion relation (also see Harrison (1908)) showing the influence of viscosity. It was known to Lamb that in addition to the discrete spectrum (governed by the viscous dispersion relation), the problem also admits a continuous spectrum (see discussion around equation 26 in Lamb (1932), chapter XI, article 349).

Figure 3 . 1 :

 31 Figure 3.1: Unbounded domain for analytical solution

  anticipating the linearised kinematic boundary condition at the interface equation 3.9. The axisymmetric Laplace equation is obtained from the continuity equation r • v p = 0 written in cylindrical axisymmetric coordinates with u p ⌘ ∂f p ∂r and v p ⌘ ∂f p ∂z . This is (Kundu and Cohen, 2002),

  likewise for the lower fluid. Combining now equation 3.33 and equation 3.34 we find,

  3.2a) in 3.93, 3.95 are constrained by the dispersion relation (which dictates the allowable value of l m ) and hence are labelled discrete. In contrast, it may be shown using a careful normal mode analysis that the other set of vorticity eigenfunctions cos p xz , sin p xz in 3.93, 3.95 are not constrained by the dispersion relation (i.e. any value of x is allowed with x 2 [0, •]) and are hence called continuous spectrum eigenfunctions. The discrete spectrum eigenfunctions decay to zero exponentially with depth (See figure3.2a) while the continuous spectrum ones, remain finite as z ! ±• (See figure 3.2b). As expected, the initial perturbation vorticity (zero in this case) is a linear sum over the discrete spectrum eigenfunctions exp[l m ]z and the continuous spectrum eigenfunctions cos p xz , sin p xz as shown in equation 3.97. One can also obtain such an equation for the upper fluid.

Figure 3 . 2 :

 32 Figure 3.2: c) Complex s plane with poles and a square root branch point for r L = 1.0, r U = 0.01, n = 0.5, k = 1, g = 1, T = 0.1

Figure 3 . 3 :

 33 Figure 3.3: Bounded domain with symmetry boundary conditions for DNS

  Figure 3.5a validates the analytical prediction made in Sec 3.2. For sufficiently small amplitudes, our axisymmetric DNS results agree very well with the planar DNS results obtained from Denner et al. (2017) (Row # 42 and # 113 in data provided in Denner et al. (2017). The axisymmetric wavenumber k of the Bessel mode was chosen to be the same as the planar wavenumber in Denner et al. (2017)). This comparison is done for a pure capillary wave with zero gravity. At larger amplitudes, the axisymmetric simulation shows deviations from the planar one. In agreement with the earlier observations of Denner et al. (2017), increasing the amplitude, decreases the frequency (Prosperetti, 1980a) of oscillation in the axisymmetric case as well (see Becker et al. (1991) and Trinh and Wang

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Comparisons between normal mode approximation, initial value problem (IVP) and Navier-Stokes simulation for Case 14, 15, 29 and 30 [i.e. a), b), c) and d)] respectively in Table 3.1.

  .102) where L is the radial domain length used in simulations. Hankel transform H [h] is done for the interface h(r, t) using an in-house subroutine written in MATLAB (2015). The code was validated by reproducing the results given in (Mathematica, 2017) for L = 1000 (L ! • recovers the Hankel transform). The numerical noise in the evaluation of equation 3.102 reduces as L ! •. As the domain size [0, 7.0155] ⇥ [0, 7.0155] produced a very noisy transform, we performed additional larger domain simulations ([0, 35.3323] ⇥ [0, 35.3323], cases 24, 25 and 26 in Table

0 Figure 3 . 6 :

 036 Figure 3.6: Interface with vorticity contours for cases 1, 2 and 3 in Table 3.1. For visualisation, the interface is mirrored about the axis of symmetry. From left to right e = 0.01, 0.1 and 0.5 respectively for the columns. Note the onset of 'jetting' in the third column.

0 Figure 3 . 7 :

 037 Figure 3.7: Interface with vorticity contours for cases 4, 5 and 6 in Table 3.1. From left to right e = 1.0, 1.2 and 1.5 respectively for the columns. Due to strong jetting in case of e = 1.5, the domain size had to be enlarged. Domain size independence of results has been verified for all cases.

0 Figure 3 . 8 :Figure 3 . 9 : 1 Figure 3 . 10 : 0 Figure 3 . 11 :

 0383913100311 Figure3.8: Vorticity ŵ ⌘ w q /J 1 (kr) vs z for cases 1, 2 and 3 in Table3.1. From left to right, the columns are for e = 0.01, 0.1 and 0.5 respectively. The IVP solution from equation 3.93 and 3.95 and computed in Mathematica.

0 Figure 3 . 12 :

 0312 Figure 3.12: Hankel Transform for Cases 24, 25 and 26 in Table3.1 corresponding to the interface profiles in figure3.11. Note that the vertical scale of the figures are not the same.

Figure 3 . 13 :

 313 Figure 3.13: Vorticity for case 26 at t = T 0

Figure 3 . 14 :

 314 Figure 3.14: Initial disturbance h(r, 0) = a 0 e r 2 (1 r 2 )

4 . 1

 41 reduces to[START_REF] Bohr | Determination of the surface-tension of water by the method of jet-vibration[END_REF][START_REF] Fyfe | Surface tension and viscosity with lagrangian hydrodynamics on a triangular mesh[END_REF][START_REF] Rayleigh | On the instability of jets[END_REF]

Figure 4 . 1 :

 41 Figure 4.1: a) A perturbed gas filament b) Three dimensional visualisation of the filament.

  (a) m = 2, La=7.2 ⇥ 10 5 (b) m = 3, La=7.2 ⇥ 10 5 (a) m = 2, La=10 2 (b) m = 3, La=10 2 of m and La, a finite number of roots are found by numerically solving 4.38 in Matlab (MATLAB, 2018) using vpasolve. Among these roots, a single conjugate pair is found to satisfy the constraint <( R0 ) > 0 (see discussion above equation 4.26 for <(l) > 0 which translates non-dimensionally to <( R0 ) > 0). This is shown on the complex R0 plane in figures 4.2a-4.2b for various m and La. For all values of La and m examined here, we numerically find a single conjugate pair of roots satisfying the constraint <( R0 ) > 0 (roots in green in figure 4.2a-4.2b).

Figure 4 . 2 :

 42 Figure 4.2: Roots of 4.38 on the complex R0 plane for various m and La. Roots in red do not satisfy the constraint < R0 > 0.

3 Figure 4 . 3 :

 343 Figure 4.3: Values of ŝ ⌘ s w 0 on the complex plane for various values of La and a) m=2 b) m=3.

3 . 1 , 2 Figure 4 . 4 :

 31244 Figure 4.4: a), b) DS eigenmodes c) CS eigenmodes

+ 1 mg

 1 and g ⇡ 0.57721 is the Euler-Mascheroni constant and Y(1) = g. Expression 4.70 evaluated for specific cases of m = 1, 2, 3 and 4 are provided in the Appendix I.6. It is clear from these expressions that for m > 0, K m (z) diverges as z m at z=0 and in addition z=0 is also a logarithmic branch point. In order to show that ã(s) in 4.69 has a branch point at s = 0, we use the Series function in Mathematica (Wolfram Research, Inc., 2017) and expand ã(s) for m = 2, 3 etc. To prevent lengthy expressions we set µ = p 0.2, R 0 = 1.0, T = 1, r = 1 with m = 2, 3 in 4.69 (these parameters correspond to Cases 1 and 3 in table 4.2) for obtaining these expansions, although the conclusions are expected to remain true for other values as well. The expansions are (subscripts on ã(s) indicating the value of m) ã2 (s) a 0 = 0.894427 + ( 0.75 ln(s) 0.896309)s + O

Figure 4 . 5 :

 45 Figure 4.5: Discontinuity of = ( ã(s)) in 4.69 on the negative real axis, for m = 2 and 3. The peaks are due to the poles of ã2 (s) and ã3 (s). Parameters for ã(s) have been chosen corresponding to Case 1 and 3 in table 4.2.

Figure 4 .

 4 Figure 4.6: a) Simulation geometry b) Benchmarking of the DNS code for equation 4.2.

Figure 4 . 7 : 0

 470 Figure 4.7: For NMA a) C = 0.5 + 0.57I and b) C = 0.5 + 0.17I

Figure 4 . 8 :

 48 Figure 4.8: For NMA c) C = 0.5 + 0.66I and d) C = 0.5 + 0.17I

5 ϵ = 0

 50 oooo o o o o o o o o o oooo o o o o o o o o oo o ϵ = 0.1 * ϵ = 0.3 + ϵ = 0.* * ***** * * * * * * * * ** *** * * * *

Figure 4 . 9 :

 49 Figure 4.9: Effect of change of a) e b) r r

Figure 4 .

 4 Figure 4.10: DNS vorticity contours outside the filament at half inviscid time period.

Figure 4 . 11 :

 411 Figure 4.11: Vorticity contours outside the filament at half inviscid time period from the analytical prediction in 4.56 and 4.77.

*

  * * * * ** *** ******* *** *** *** *** ***** *****************+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + * * * ** ******** ** ** *** **** ************ + + + + + + + + + + ++ ++++++ ++ ++ ++ +++ +++++++++ DNS o Eq. 4.77 * DS + * * * * * ** *** ****** *** *** **** ******* + + + + + + + + + + + ++ +++++++++ ++ ++ ++ +++ +++++ DNS o Eq. 4.77 * DS + o o o o o o o o o o o o o o o o oo oo ooo oooooooo ooo

Figure 4 . 12 :

 412 Figure 4.12: Radial vorticity versus radius at q = p/4 at different times -Case 1 in table 4.2.

  and simulations", Sagar Patankar, PalasKumar Farsoiya and Ratul Dasgupta, J. Fluid. Mech, 857, 80-110, 2018. The numerical Floquet analysis for obtaining the stability charts was performed by Sagar Patankar and details of this is available in Patankar et al. (2018). We only discuss the stability charts in this chapter and refer the reader to Patankar et al. (2018) for details.

Figure 5 . 1 :

 51 Figure 5.1: Cylindrical filament with radial forcing in both fluids. The surface indicated with lines is the unperturbed filament.

3 .

 3 Equation 5.3 thus predicts that instability analogous to Faraday waves on Cartesian (Benjamin and Ursell, 1954) and spherical geometry (Adou and Tuckerman, 2016), are also possible on a cylindrical filament. On the B-A plane, the stability chart of equation 5.3 i.e. the boundaries separating the stable bounded solutions from the exponentially growing ones, can be determined either using perturbative methods for B << 1 or through numerical Floquet analysis provided in Patankar et al. (2018).

Figure 5 . 2 :

 52 Figure 5.2: The neutral stability curves of equation 5.3. The tongues separate stable from unstable regions, S represents stable, SH subharmonic and H harmonic. The dots (see inset on the right) are obtained from Bender and Orszag (2013). Note the small hashed region which is stable. We acknowledge Sagar Patankar for obtaining this chart and refer the reader to Patankar et al. (2018) for details.

Figure 5 . 3 :

 53 Figure 5.3: Free perturbations: Case 1, 2, 3 and 4 in Patankar et al. (2018). a) Comparison of DNS results with predictions from expression 5.1. b) Generation of droplets through RP instability for Case 1 in a).

Figure 5 . 4 :

 54 Figure 5.4: a) DNS geometry b) A sample adaptive grid. c) An axially perturbed column at t = 0 i.e. h(x, 0) = a(0) cos(kx). d) An azimuthally perturbed column at t = 0 i.e. h(y, z, 0) = a(0) cos(mq) with q ⌘ tan 1 (z/y).

Table 5 . 1 :

 51 Boundary conditions for DNS in Fig.5.4a.The simulations have been performed using the open source code Gerris[START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] (our in-house developed code is a two-dimensional solver in comparison to Gerris which can solve both two-dimensional and three-dimensional flows) which solves the incompressible Euler equations with surface tension and an additional time periodic body force, chosen to be in the radial direction in cylindrical coordinates. As a preliminary validation of the code, we have simulated axisymmetric (m=0) free oscillations (kR 0 > 1) as well as the RP instability (kR 0 < 1) by setting h = 0 in equation 5.1 (also see earlier validations in (Popinet, 2013)). In figures 5.5a and 5.5b, we compare DNS data with numerical solutions to equations 5.3 using the MATLAB (2015) solver ode45 and the perturbative solutions given by expressions (Patankar et al., 2018).

Figure 5 .

 5 5a corresponds to the linearly stable case 5 in table 6.1 while figure5.5b corresponds to the linearly unstable case 6. In both cases, the simulation parameters were chosen to ensure that B << 1 thus enabling comparison with the perturbative solutions (Patankar et al., 2018). Note that while solving the Mathieu equation 5.3, we have used the initial conditions (in collaboration withSagar Patankar)

Figure 5 . 5 :

 55 Figure 5.5: a) Comparison of DNS results (case 5 in table 5.2) with numerical solution to equation 5.3 and perturbative solution Patankar et al. (2018) b) Comparison of DNS results (case 6 in table 5.2) with numerical solution to equation 5.3 and perturbative solution Patankar et al. (2018).

Figure 5 .Figure 5 .

 55 Figure 5.6: a) Stability chart with case numbers from Table 5.2 b) DNS for case 1

Figure 5 .Figure 5 . 9 :Figure 5 . 10 :Figure 5 .

 5595105 Figure 5.8: a) Case 9: Linearly unstable b) Case 10: Linearly stable

Figure 5 Figure 5 .

 55 Figure 5.12: Case 15

Figure 5 .Figure 5 . 14 :

 5514 Figure 5.14: For case 6 in table 5.2 a) Instantaneous interface profiles on the x-y plane b) z-y plane. c) Instantaneous streamlines for the profile corresponding to t = 14p.

Figure 5 . 15 :Figure 5 . 16 :

 515516 Figure 5.15: For case 12 in table 5.2, instantaneous interface profiles on the a) x-y plane b) z-y plane

Figure 5 . 17 :Figure 5 .

 5175 Figure 5.17: For case 16 in table 5.2, instantaneous interface profiles on the a) x-y plane b) z-y plane c) Instantaneous streamlines for the profile corresponding to t = 3p

Figure 5 .Figure 5 . 20 :

 5520 Figure 5.19: a) Effect of change of e az . b) Effect of change of e ax .

  Figure 6.1: Free perturbation geometry a) h(x, 0) = a 0 exp

  the linearised limit, this problem was analytically first solved by Cauchy and Poisson (Cauchy, 1827; Poisson, 1818) in the inviscid, irrotational approximation (hence the name Cauchy-Poisson problem (Debnath, 1994)) by expressing the solution as a superposition over all the Fourier modes excited at t = 0. The solution provided here generalises the one given in Debnath (1994) including (linearised) inertial effects of both fluids and surface-tension. With h0 (k) ⌘ F(h(x, 0)) where F represents the Fourier transform, the interface shape and the velocity potential in both the fluids at any later time t is given by (Singh et al., 2019b),

0 = m 0 m 2 0 1 T 0 = m 0 m 2 0

 02100 .1b. The circular droplet of radius R 0 (dotted lines) shown in the figure represents the cross-section of a cylindrical filament in three dimensions. Ignoring the inertia of the outer fluid, Rayleigh (1879a) had shown that pure azimuthal Fourier modes of wavenumber m 0 on the circular interface are stable and produce oscillations of frequency dictated by w 2 /r I R 3 0 . The effect of the ambient fluid was included by Bohr (1909) who obtained the modified dispersion relation w 2

(Figure 6 . 2 :

 62 Figure 6.2: Forced perturbation simulation geometry a) Faraday waves on Cartesian geometry (not to scale) b) Faraday waves on a circle. Arrows indicate forcing. c) For the circular case, only one-fourth of the drop is simulated, see table 6.1.

  R 0 (subscript b for base state variables, see Singh et al. (2019b)). For the cylindrical filament, the existence of three-dimensional Faraday waves on the filament, has been demonstrated in chapter 5. The case of purely azimuthal perturbations was not discussed earlier and is treated here. It has been shown in Singh et al. (2019b) that for a standing azimuthal Fourier mode imposed on the circular interface, of the form a(t) cos(mq) (m being constrained to be a positive integer due to the periodicity requirement), the equation governing a(t) is the Mathieu equation (the derivation was obtained by the M.Tech student and is available at Singh et al. (2019b)), .

(Figure 6 . 3 :

 63 Figure 6.3: Case 1 At t = 10T 0 , See Table 6.2 for We values.

Figure 6 . 4 :Figure 6 . 5 :Figure 6 . 7 :

 646567 Figure 6.4: ! FSS, ! Analytical solution from equation 6.4 , x ! OFS, 4 ! Basilisk-centered Popinet (2014)

Figure 6 . 13 :Figure 6

 6136 Figure 6.12: a) Stability chart parameters: r r = 1000, R 0 = 1.0, T = 1.0 ! FSS, ! Mathieu solution: equation 6.16, x ! OFS, 4 ! Basilisk Popinet (2014)

Figure I. 1 :

 1 Figure I.1: Contour for Bromwich integral

Figure I. 2 :

 2 Figure I.2: Cylindrical filament with radial forcing in both fluids. The surface indicated with lines is the unperturbed filament.

0 G

 0 these equations along r we obtain (one of the integration limits is R 0 )P I b (r, t) = P I b (R 0 , t) r I Z R 0 r G(r 0 , t)dr 0 , P O b (r, t) = P O b (R 0 , t) + r O Z r R (r 0 , t)dr 0 . (I.18)Due to surface tension coefficient T, the pressure jump boundary condition across the interface in the base state isP I b (R 0 , t) P O b (R 0 , t) = T/R 0 .The perturbation velocity potential satisfies the Laplace equationr 2 f I = r 2 f O = 0 (I.19)With h(q, z, t) defined as the difference between the perturbed and the unperturbed surfaces (see fig. I.2), The linearised kinematic boundary condition is the unit vector in the radial direction of the unperturbed filament.Using variable separability, the interfacial perturbation and the perturbation velocity potentials are chosen to have the form of standing wave h(z, q, t) = a(t) cos(mq) cos(kz) (I.21)f I (r, z, q, t) = L I (r) cos(mq) cos(kz) r, z, q, t) = L O (r) cos(mq) cos(kz) into the Laplace equation, we obtain equations for L. i.e. verified that the solutions to equations I.24 are (see expression 6.2.28 in Prosperetti (2011))L I (r) = C I I m (kr) + D I K m (kr), L O (r) = C O I m (kr) + D O K m (kr) (I.25)where I m (kr), K m (kr) are the modified Bessel functions of mth order, of the first and second kind respectively. We would like the solution to stay bounded at r = • and r = 0. This implies that (Weisstein, 2017a,b) L I (r) = C I I m (kr) and L O (r) = D O K m (kr). In order to determine C I and D O , we use this in equation I.20 alongwith I.23 to obtain (rf • êr in cylindrical coordinates is ∂f ∂r )

2 || 2 r 2 + k 2 ⌘

 2222 40 and 5.41 in Bush (2013)). Using the decomposition in cylindrical coordinates r 2 = r equation I.29, upto linear order in a(0) we have,(r • q) r=R 0 +h ⇡ ⇣ r 2 Q ⌘ r=R 0 +h = c(R 0 + h) a(t)r 2 || cos(mq) cos(kz) z, t) a(t)l(R 0 )F(z, q).where c(r) ⌘ r 2 ? (r), l(r) ⌘ ⇣ m and c(R 0 ) is an O(1) term arising from base state curvature with 1 R 0 = c(R 0 ) through h. In writing the second line after equation I.30, we have Taylor expanded c(r) about R 0 retaining terms upto

  28 with I.30, I.31, I.32 (after cancelling out the base state contribution to pressure using 1R 0 = c(R 0 )), we obtain, ⇣ r I r O ⌘ G(R 0 , t)h(q, z, t) cos(mq) cos(kz) (I.33) The (linearised) Bernoulli's equation for the perturbation pressure at the linearised interface (Farsoiya et al., for perturbation pressures p from equation I.34 in equation I.33, using I.23 and setting the periodic forcing in both fluids to be of the form G(r, t) = h cos(Wt) r R 0 (Adou and Tuckerman, 2016), we obtain a Mathieu equation of the form .. a(t) + f (t)a(t) = 0 is the Mathieu equation governing Faraday waves on a cylindrical filament.

  order behaviour at large z, K m (z)

  of K m (z), Chapter 4

Figure I. 3 :

 3 Figure I.3: Bromwich contour for obtaining the Laplace inversion of 4.74 and 4.76. The two poles which lie on the principle sheet are complex conjugates of each other viz. s 1 = s2 .

L 2 For m 2

 22 (d exp(Iq)) Q(d exp(Iq)) exp [d exp[Iq]t] d dq (I.62)Interchanging the limit and the integral in I.62, we obtain an expression for lim d!0L(d exp[Iq]) Q(d exp[Iq]) . For |s| << 1, as d ! 0, M (d exp[Iq])) ⇡ 2, c (d exp[Iq]) Consequently I 4 vanishes in the limit d ! 0 in I.62.

,

  exp [st] ds = 2 Â i=1 L(s i ) Q 0 (s i ) exp [s i t] (I 3 + I 5 ) (I.64) which is the form written in 4.74. Note that as all roots of Q(s) lie on the left half plane, any c > 0 may be chosen in I.64. The Laplace inversion of vorticity given by 4.76 may be obtained using the same contour as figure I.3 and an analogous procedure. Q ⌘ s 2 + sM(s) + w 2 c(s) and the expressions for M(s), c(s) are provided below equation 4.69. The expressions for A and B are provided below 4.77.

0 m 0 m 0 m 1 0 m 1 G

 000101 J or Y. Subtracting the two we obtain,K m (z exp (pI/2)) = p 2 exp ( mpI/2) [Y m (z) + IJ m (z)] (I.70) Using the identify K (z) = K m+1 (z) + m z K m (z), we have K (z exp (pI/2)) = K m+1 (z exp (pI/2)) + m z exp (pI/2) K m (z exp (pI/2)) 70 and I.71 in I.73, we obtain c(x exp (Ip)) = Y (a) + IJ (a) m+1 a ⇢ Y m 1 (a) + IJ m 1 (a) the second step we have used I.69. We also have using I.70 continuous spectrum eigenfunctions for hollow gas filament, chapter 4 191We also obtain using I.74 and I.76, Q (x exp(Ip)) ⌘ h s 2 + sM(s)

Table 1 . 1 :

 11 Some landmark studies of inviscid interfacial waves

	Author(s)	Remarks
	Lamb (1932)	Free surface Cartesian interface with normal
		mode analysis
	Chandrasekhar (1961a) Two fluid normal mode analysis in Cartesian and
		Spherical geometry.
	Miles (1968)	Free surface axisymmetric Cauchy-Poisson prob-
		lem. Mentioned the continuous spectrum but
		did not discuss its significance.
	Present

Thesis and Far- soiya et al. (2017)

Table 1 . 2 :

 12 Studies of viscous IVP and normal mode analysis in interfacial waves

et al. (2011)),

  The incompressible Navier-Stokes equations with surface tension included as a body force are (see equation 2.2 in Tryggvason

	Test SLIC(E)	Hirt-	FCT-VOF(E) Youngs(E)	LVIRA(E)	LVIRA
			Nichols(E)			D f (%)
	1 2 3 4	1.30 ⇥ 10 2 4.55 ⇥ 10 2 1.28 ⇥ 10 2 9.18 ⇥ 10 2 1.9 ⇥ 10 1 3.99 ⇥ 10 2 4.59 ⇥ 10 2 6.66 ⇥ 10 2 3.14 ⇥ 10 2 8.38 ⇥ 10 2 9.62 ⇥ 10 2 3.29 ⇥ 10 2	3.08 ⇥ 10 3 1.5 ⇥ 10 3 2.98 ⇥ 10 2 1.05 ⇥ 10 2 8.60 ⇥ 10 3 6.90 ⇥ 10 3 1.09 ⇥ 10 2 9.7 ⇥ 10 3	0.14 1.33 0.93 0.49
		r • u = 0				(2.22)

Table 2 .

 2

1: Errors for various tests (Column 2-4 from Rudman (1997)) Figure 2.13: Integral over a control volume

  After adding it with its transpose, we are left with the symmetric part of the velocity gradient combining with axisymmetry (∂f/∂q = 0, w = 0(zero swirl)),

	we obtain,
	ru + ru T =	2 6 6 6 6 6 6 6 6 6 4	2 ∂v ∂R	∂u ∂R 0 + ∂u ∂z	0 2 u R 0	∂u ∂z 2 + 0 ∂z ∂R ∂v ∂v	3 7 7 7 7 7 7 7 7 7 5
			2	∂u	1	∂u	w	∂u	3
			6 6 6 6 6 6 6 6 6 4	∂R ∂w ∂R ∂v	R 1 R	∂q ∂w ∂q 1 ∂v +	R u R	∂z ∂w ∂z ∂v	7 7 7 7 7 7 7 7 7 5
				∂R		R	∂q	∂z

Brackbill et al., 1992) enables

	The Continuum Surface Force(CSF) model (us to
	add surface tension as a body force in Navier-Stokes equations,
	F F F lg d s = T k n n nd s	(2.45)
	where dirac delta function d s keeps that force only at the interface (

Tryggvason et al., 2011).

Kalland (2008); Lörstad and Fuchs (2004) is

  To calculate the surface tension force we need to compute the local curvature of the interface k. As volume fraction f is not a smooth function, a height function is calculated from f to find curvature values of the interface. The Direction Averaged Curvature (DAC) method as described in used for estimating k. As shown in Fig.2.16 (for the case when the orientation of normal of the interface is mainly towards positive +z orientation(|N z | > |N R | and N z > 0 in axisymmetric coordinates), the other three cases can be dealt similarly). N R and N z are the horizontal and vertical components of the normal to the interface. The first principal curvature is given by,

	0.00 0.00 0.00	0.00 0.00 0.00	0.00 0.00 0.00
	0.00 0.00 0.00	0.00 0.00 0.00	0.00 0.00 0.00
	0.00	0.05 0.40	0.00	0.05 0.40	0.00	0.05 0.40	1.40
										0.93
	0.12	0.88	1.00	0.12	0.88	1.00	0.12	0.88	1.00
	0.60	1.00	1.00	0.60	1.00	1.00	0.60	1.00	1.00	-0.34
							-1.00
	0.94	1.00 1.00	0.94	1.00 1.00	0.94	1.00 1.00
							-1.00
	1.00 1.00 1.00	1.00 1.00 1.00	1.00 1.00 1.00
	0.12 0.88 1.00	0.12 0.93 1.40	-0.34 0.93 1.40
	Figure 2.16: Variable Stencil for calculation of height functions for +z orientation
	Direction averaged curvature method

< 1 else if (orientation=DOWN/LEFT/RIGHT) then . respective heights for the orientations similar to UP end if H 0

	Algorithm 1 Curvature
	procedure k(r, c,orientation)	. cell index as argument
	H le f t , H center , H right = 0 . Normal pointing upwards as shown in figure 2.16
	if (orientation=UP) then jt = r
	. sum c upward to find the height of darker fluid from the cell bottom
	do
	H le f t/center/right = H le f t/center/right + f jt,c 1/c/c+1 jt = jt + 1
	while f r,c 1/c/c+1 > 0 jt = r 1
	. sum (1 f ) to find the height of darker fluid from the cell bottom
	do
	H le f t/center/right = H le f t/center/right + ( f r,c 1/c/c+1 1) jt = jt 1
	while f r,c 1/c/c+1,j
	r • ūn+1 = 0				(2.50)
	Using equation 2.27 and 2.50, we obtain
	⇣ r • ū⇤ ū⇤ Dt P n+1 r • ūn+1 = r • r • P n+1 = 1 Dt	⌘	= 0	(2.51) (2.52)

1995),Briggs et al. (2000)

	Poisson equation (PPE) solved as a part of solving the Navier-
	Stokes equation is the slowest part of the solution. Thus it is necessary to solve
	this equation efficiently. Fast Poisson solvers use advanced techniques, such as
	discrete Fast Fourier transform or cyclic reduction, multigrid methods or iterative
	Krylov-subspace methods. Developed in late eighties (Wesseling (

The results for this problem are compared with Moin (2010) (See Fig. 2.21). 2.5.2 Two dimensional Poisson equation

	10 0
									Moin(2010)
									Present Study
	10 -2							Moin (2010)
								Present Study
	10 -4
	10 -6
	10 -8
	10 -10
	10 -12
	10 -14	0	2	4	6	8	1 0	1 2	1 4	1 6
						V cycles V-Cycles
			Figure 2.21: Comparison with Moin(2010)
	backward for prolongation step as,
			e h i = e 2h i	i = 0,...,N		(2.73)
			e h 2i+1 =	e 2h i + e 2h i+1 2	i=0,...,N		(2.74)
	While prolongation iteration step is taken at every grid level and at finest grid,
	we add the interpolated error to the u, again start with the V-cycle with the new
	u as a guess. We are interested in solving the Poisson equation, Equation 2.75, in a square
	domain of size unity centered on unity with Neumann boundary conditions on
	all sides (Popinet, 2003).
				r 2 f = f (x, y)			(2.75)
	Source term f (x, y) is given by
			2 f (x, y) = p 2 (k 2 + l 2 ) sin(pkx) sin(ply) (e n i+1 + e n+1 h 2 r i ) i 1		(2.72) (2.76)
	Restriction and iteration continue till the coarsest grid. At coarsest grid the
	iteration the error equation is solved exactly. Then the error is interpolated with k = l = 3. Exact solution of the Poisson equation with this source term is
			f(x, y) = sin(pkx) sin(ply) + k		(2.77)

Vortex Property Ghia et al. (1982) Present

				Study
	Primary Center	y min	0.117929	0.1165
		w	2.04968	2.0087
		(x, y)	(0.5313, 0.5625)	(0.5313, 0.5625)
	Bottom Left	y max w	2.31129 ⇥ 10 4 0.36175	2.345 ⇥ 10 4 0.3626
		(x, y)	(0.0859, 0.0781)	(0.0859, 0.0781)
	Bottom Right	y max w	1.75102 ⇥ 10 3 1.15465	1.7417 ⇥ 10 3 1.1138
		(x, y)	(0.8594, 0.1094)	(0.8594, 0.1094)
	Bottom Right (Secondary)	y min	9.31929 ⇥ 10 8	6.452 ⇥ 10 7
		w (x, y)	8.52782 ⇥ 10 3 (0.9922, 0.0078)	0.0096 (0.9922, 0.0078)

Table 2 .2: Stream function and vorticity 2.6.1 Lid driven cavity

 2

	Lid driven cavity is considered a classic test case for testing implementation of
	one phase Navier-Stokes equations in CFD. Here we have a square cavity of side
	length unity and no slip conditions on three boundary except the top which
	moves with constant velocity. The normal derivatives of pressure is kept zero
	the incompressible NS equations are integrated until a steady-state is reached
	(at Re = uL/n = 1000). We benchmark velocity profiles obtained from the solver
	against the data

of Ghia et al. (1982) (See figure 2.25 and Table 2.2). 2.6.2 Regularised lid driven cavity

	This test case was suggested to the author by Prof. John Hinch during a GIAN
	course Hinch (2018) conducted at IIT Bombay. The principle reason for this is
	the stress singularity at the top corners of the ordinary lid-driven cavity due to
	discontinuity in tangential velocities. In order to resolve this, the velocity of the
	lid was modified from a constant value of unity to a spatially dependent function
	sin 2 (px) which goes to zero at both top corners (See figure 2.26). See Table 2.3
	for the value of tangential force on the lid at y = 1 and figures 2.27 for contours
	of pressure and vorticity.
	2.6

.3 Static circle test The

  curvature calculation method section 2.3.3 in the in-house solver is verified from the static circle test. This test (cf.

Tryggvason et al. (2011), sec

. 7.3.1 for

Table 2 .3: Force on the lid

 2

	1.0
	0.8
	0.6
	0.4
	0.2
	0.0 0.0	0 .2	0 .4	0 .6	0 .8	1 .0
			x

  ⇥ s 1 e s 2 t s 2 e s 1 t ⇤ for s 1 6 = s 2 ,(3.100) where s 1 and s 2 are the two roots of equation 3.99. Equation 3.99 is solved in Mathematica (Wolfram Research Inc., 2010) to determine the roots. For overdamped oscillations s 1 and s 2 are real while in the underdamped case they occur as a complex conjugate pair. In order to test the validity of the linearised IVP solution in equation 3.85 and the normal mode approximation in equation 3.100, we turn to numerical simulations.

	a(t) =	a(0) s 1 s 2

3.87

. Such a dispersion relation in planar geometry taking into account (linearised) inertia and viscosity in both fluids is provided in Chandrasekhar (1961b). As remarked earlier, since ã(s) turns out to have the same expression in Cartesian and cylindrical geometry (Farsoiya et al., 2017), it is expected that the dispersion relations are also the same. It may be checked numerically that equation 3.99 has two roots only (See figure

3

.2c where they appear as poles in s complex plane). We obtain what we call a normal mode approximation, by assuming that the amplitude of the standing wave a(t) may be expressed as C 1 exp (s 1 t) + C 2 exp (s 2 t) and determining C 1 , C 2 from initial conditions a(0) and . a(0) = 0. This approximation has the form,

  3.100), the IVP solution in equation 3.85 inverted numerically (McClure, 2016) and results from DNS is presented in figure 3.4. Note that t ⌘ tw 0 with w 2 0

et al. (2017).

  .5a, one does not observe jetting at the non-dimensional parametric values reported in Denner

	While observations of jetting
	(Worthington jets) have been reported earlier in the context of standing waves
	hitting a wall (Longuet-Higgins, 2001), standing waves generated through forced
	oscillations (

Zeff et al., 2000), droplet impact (Bartolo et al., 2006) and

  water entry of axisymmetric objects (

Gekle and Gordillo, 2010), observation

  of jetting in purely free oscillations is reported here for the first time. In figure3.5b, as e is increased further upto a value of 1.2, a droplet pinches off from the jet. This is referred to as jet break-up and droplet pinch-off (see[START_REF] Leal | Advanced transport phenomena: fluid mechanics and convective transport processes[END_REF], page 83 for discussion of a capillary mechanism of end pinch-off in cylindrical threads).

	Remarks	Case	e	r r	µ r	La	Ga	Domain	Grid
	e r r µ r La Ga Hankel Trans-form Normal mode 1 Comparison et (2017) al. Denner with	1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 29 30 31 32	0.01 0.1 0.5 1.0 1.2 1.5 0.01 do do do do do do do do do do do do do do do do 0.5 1.5 3.5 0.01 0.01 0.062 1.25	100 do do do do do 1 10 1000 100 do do do do do do do do do do do do do do do do 10 2 do 1 do	100 do do do do do do do do 0.1 1 10 10 3 10 4 2 ⇥ 10 4 5 ⇥ 10 4 10 5 10 2 do do do do do do do do 125 83 1 do	10 6 do do do do do do do do do do do do do do do do 5.0 ⇥ 10 6 1.0 ⇥ 10 7 1.5 ⇥ 10 7 10 6 do do do do do 6.25 2.78 126 do	10 5 do do do do do do do do do do do do do do do do do do do 5.0 ⇥ 10 5 1.0 ⇥ 10 6 1.5 ⇥ 10 6 10 5 do do 0.625 0.278 0 do	L 2 1 do do do do L 2 2 L 2 1 do do do do do do do do do do do do do do do do L 2 3 do do L 2 1 do L 2 4 do	512 2 do do do do 1024 2 512 2 do do do do do do do do do do do do do do do do 1024 2 do do 512 2 512 2 256 2 do
	e	r r	µ r

Table 3 .1: DNS Parameters

 3

Tripathi et al. (2014) and Dixit and

Castro-Hernández et al., 2011; Van Hoeve et al., 2011).

  The breakup of these cylindrical liquid or gaseous jets is often attributed to the Rayleigh-Plateau instability (

. Surface oscillations, instability, deformation and breakup of such cylindrical gaseous filaments (density r I ) surrounded by an immiscible fluid (density r O ) are of interest to numerous scientific applications e.g. microbubbles are obtained from breakup of a gaseous jet injected into a flowing liquid inside a T-junction (Garstecki et al., 2006) or a coaxial device (

Plateau, 1873; Rayleigh, 1892a; Van Hoeve et al., 2010) or

  to the pressure drop across the injected gas bubble (

Castro-Hernández et al., 2011; Garstecki et al., 2006).

  In order to simplify analysis, stability studies of cylindrical jets or filaments 1 are frequently conducted under one of two approximations for the density ratio. The first approximation r O /r I ! 0 de- scribes a liquid jet issuing into an ambient gas (like air) of significantly lower density. The converse limit r O /r I ! • corresponds to the so-called 'hollow jet' approximation (

Chandrasekhar, 1981; Eggers and Villermaux, 2008; Rayleigh, 1892a) and

  applies to gaseous jets issuing into a (much) denser fluid, a recurrent situation in many microfluidic and chemical engineering applications. The Rayleigh-Plateau instability (RP hereafter) is driven by surface tension and affects liquid as well as gaseous filaments. Inviscid, irrotational, linearised analysis (see (Eggers and Villermaux, 2008), and Lin (2003) for exhaustive reviews) predicts that axisymmetric interfacial Fourier modes with wavelength greater than the circumference of the unperturbed filament are unstable and can lead to breakup. This may be seen from the dispersion relation 4.1 governing three dimensional, irrotational perturbations on a cylindrical filament separating two inviscid, immiscible fluids (

Meister and Scheele, 1967; Patankar et al., 2018; Rayleigh, 1878, 1879b, 1892a),

et al. (1995); Moon et al. (2016); Netzel et al. (1964); Ronay (1978)

  , zero frequency modes respectively and are not of interest). The study of such azimuthal capillary perturbations (the topic of the present study) have led to interesting applications e.g. Rayleigh (1889) proposed the use of 4.2 to estimate changes in dynamic surface tension of an elliptic liquid jet composed of water and oleate mixed in various proportions (see Pederson (1907)). Neils Bohr (Bohr, 1909) analytically determined the effect of viscosity (and individually that of weak nonlinearity or density of the ambient air) to 4.2. A number of experimental and theoretical investigations have subsequently examined the accuracy and limitations of the modifications to equation 4.2 proposed by Bohr (see Bechtel and page 72 of the review by Eggers and Villermaux (2008)). In many applications, inviscid-irrotational analysis may not be sufficient and it is necessary to incorporate viscous effects into the relations 4.1 and 4.2. The dispersion relation governing small amplitude, axisymmetric perturbations on

.2) It is seen from 4.2 that azimuthal modes with m = 2, 3 . . . lead to oscillatory motion at the interface (in two dimensions m = 0 and 1 represent compressible and translationala viscous liquid cylinder was first obtained by Rayleigh (1892b) (also see Bauer (1984); Liang et al. (2011); Stone and Brenner (1996); Tomotika (1935) for further results). For axisymmetric perturbations on a hollow filament surrounded by a viscous liquid, the corresponding dispersion relation is presented in Chandrasekhar (1981) (also see Parthasarathy and Chiang (

Cortelezzi and Prosperetti, 1981; Lamb, 1932; Prosperetti and Cortelezzi, 1982; Prosperetti

	The possibility of the continuous spectrum for axisymmetric perturbations on a
	viscous liquid filament was first studied by Berger (1988) by solving the initial
	value problem with mixed initial conditions involving a surface deformation and
	an impulse (see section 2.2 in García and González (2008) for a survey of refer-
	ences on the inviscid, initial-value problem). The study by García and González
	(2008) has however negated some of Berger (1988)'s conclusions showing that
	the corresponding viscous dispersion relation admits infinite roots leading to a
	countably infinite set of eigenmodes. The results of García and González (2008)
	are consistent with similar conclusions drawn by Prosperetti (1980a) earlier, for
	capillary perturbations on a liquid droplet. While the normal mode approach
	has been compared to the IVP solution for a liquid filament (

, 1981) or cylindrical quiescent, viscous pool of liquid (Farsoiya et al., 2017) of infinite depth, capillary standing waves on a bubble surrounded by radially unbounded viscous liquid (Prosperetti, 1980a,b)) etc. In addition to the problems mentioned above, the continuous spectrum has also been studied while investigating the stability of thin film flows on a horizontally unbounded domain with bottom topography (Kalliadasis and Homsy, 2001).

Berger, 1988; García and

  We validate the normal mode predictions by solving the linearised initial value problem (IVP) in section 4.55.

	peretti, 1980a, 1976, 1980b, 1981), especially at intermediate times. We validate
	our normal mode analysis by further solving the initial value problem. All analyt-
	ical predictions are tested against Direct Numerical Simulations (DNS) conducted
	using an in-house developed code.
	This chapter is organised as follows: The linearised equations and boundary
	conditions are presented in section 4.2. We perform normal mode analysis in
	section 4.3.1 obtaining the viscous dispersion relation governing the discrete

González, 2008), this has not been done for a hollow filament surrounded by a viscous fluid. In particular, the existence of the continuous spectrum in addition to the discrete spectrum has not been discussed in the literature before. Knowledge of both the spectrum bears importance towards understanding and modelling of transient free, oscillations of a gaseous filament or a jet. In this study, we analyse the discrete and the continuous spectrum for purely azimuthal perturbations on a hollow cylindrical filament. Our analysis complements and extends earlier results on the continuous spectrum by Lamb (1932),

[START_REF] Prosperetti | Viscous effects on small-amplitude surface waves[END_REF][START_REF] Prosperetti | Motion of two superposed viscous fluids[END_REF]

and Farsoiya et al. (2017) in other geometries, to the present case of a hollow cylindrical filament. We demonstrate analytically the existence of the discrete and the continuous spectrum. It is shown that an initial condition in the form of an interfacial Fourier mode and zero perturbation vorticity, has projections on the discrete as well as the continuous spectrum eigenfunctions. Consequently, it produces standing waves whose temporal evolution is not well approximated by a sum of exponentials (Farsoiya et al., 2017; Prosspectrum. Using normal mode analysis we further demonstrate in section 4.3.2, the existence of an additional, uncountable infinite set of radial eigenmodes constituting the continuous spectrum. The discrete and the continuous modes together lead to a completeness relation for radial vorticity in section 4.

3.3.

  .9)As shown in figure4.1a, h(q, t) represents the interfacial perturbation. The boundary conditions to be enforced are the kinematic boundary condition, zero shear stress and jump in normal stress due to surface tension, all at the interface. Additionally because our domain is radially unbounded, there are boundedness conditions for all variables, as r ! •. The mathematical expressions for these are respectively[START_REF] Leal | Advanced transport phenomena: fluid mechanics and convective transport processes[END_REF]

	and	✓ ⇣ t • s tot • n Dh Dt u tot • n ⌘ r=R 0 +h ◆ r=R 0 +h = 0, ⇣ ⌘ n • s tot • n Tr • n r=R 0 +h = 0,	= 0,	(4.10) (4.11) (4.12)

Table 4 .

 4

1: Sample roots of 4.38 satisfying <( R0 ) > 0

  Parameters for ã(s) have been chosen corresponding to Case 1 and 3 in table 4.2.

	ã3 (s) a 0	= 0.596285 0.05s + (0.124226 ln(s) 0.122629)s 2 + O	⇣	s 3	⌘	. . .	(4.72)
	. . .	. . .	. . .	. . .	. . .	. . .			. . .	. . .

Note the present of terms of the form s m 1 ln(s) in 4.71 and 4.72 which have branch points at s = 0 and s = •. These branch points are inherited from K m ( p s n R 0 )

Table 4 . 2

 42 All DNS parameters are summarised in table 4.2. All viscous simulations have been performed choosing density ratio r = 10 4 and viscosity ratio in the range µ ⇡ 10 3 and 10 4 . For minimising nonlinear effects, we chose e = 0.1. Note that in table 4.2, the values for viscosity of the outer fluid (µ O ) were chosen to ensure integral values of La. All simulations are initiated with an interfacial perturbation of the form r = R 0 + a 0 cos(mq) with zero velocity everywhere in the domain. The resultant motion of the interface is tracked in time and the vorticity field is computed from the velocity field obtained in DNS. As a first step,

	Case	Grid			a 0		r O		µ O	r I	µ I	T R 0 m La⌘	TR 0 r µ 2	e ⌘ a 0 m R 0
	0 1 2 3 4 5 6 7 8 9 10	256 ⇥ 256 0.05 1 256 ⇥ 256 0.05 1 256 ⇥ 256 0.05 1 256 ⇥ 256 1/30 1 256 ⇥ 256 1/30 1 256 ⇥ 256 0.05 1 256 ⇥ 256 0.15 1 256 ⇥ 256 0.25 1 256 ⇥ 256 0.05 1 256 ⇥ 256 0.05 1 256 ⇥ 256 0.05 1	-0.2 10 4 10 4 0.5 1 2 10 4 -0.5 1 2 2/20 10 4 10 4 0.5 1 2 p p p 0.2 10 4 10 4 0.5 1 3 p 2/20 10 4 10 4 0.5 1 3 0.01 10 3 10 4 1 1 2 0.01 10 3 10 4 1 1 2 0.01 10 3 10 4 1 1 2 0.01 0.01 10 4 1 1 2 0.01 0.1 10 4 1 1 2 0.01 1 10 4 1 1 2	• 2.5 100 2.5 100 10 4 10 4 10 4 10 4 10 4 10 4	0.1 0.1 0.1 0.1 0.1 0.1 0.3 0.5 0.1 0.1 0.1
		u = 0,	∂v ∂x	= 0,	∂p ∂x	= 0,	∂c ∂x	= 0,	for left and right boundary	(4.78)
	v = 0,	∂u ∂y	= 0,	∂p ∂y	= 0,	∂c ∂y	= 0,	for top and bottom boundary (4.79)
	on the four sides of the domain.

O

, viscosity ratio µ ⌘ µ I /µ O and e ⌘ a 0 m/R 0 , a nonlinearity parameter. In order to compare : DNS Parameters results obtained from DNS with the linearised theory presented earlier, we require the nonlinearity parameter e << 1. In addition we need to ensure that the density and viscosity ratio are small i.e. r, µ << 1, in order to keep dynamic effects of the fluid inside the filament to a minimum. The boundary conditions applied are, we have validated the inviscid dispersion relation 4.2 by imposing an azimuthal Fourier mode m = 2 in our solver and integrating the inviscid equations of motion (parameters corresponding to case 0 in table 4.2). As shown in figure

4

.6b, the time signal obtained from DNS is fitted with a(t)/a 0 = cos(w 0 t) and

  Comparing 4.7a and 4.7b, it is seen that the solution to the IVP given by 4.75 agrees very well with DNS. In contrast, the discrete part of expression 4.75 (labelled as discrete spectrum (DS)) shows a poor agreement at t = 0 especially at large values of viscosity (smaller La). The same behaviour is also observed in figures 4.8a and 4.8b for a Fourier mode with m=3. We also compare with an approximation

	* * * * * * * * * * * * *** * * * * * * * * * * **** * * * * * * * * * * ** * * * * * * * * * * ***** * * * * * * * * ** ** * * * * * * * * * ** *** * * * * * * * * ** *** * * * * * * * * * + oo o + + + + + + + + + + o ++ + + + + + + + + + + ++ + + + + + ++ + + + o + +++ + o + + + o + o + ++ + o ++ o + + o o ooo o

shows good agreement. For figures

4.7a, 4.7b, 4.8a and 4.8b

, all time signals are obtained by tracking the interface at q = 0 (see figure

4

.6a), although choosing any other q would lead to identical signals. We compare data obtained from DNS with predictions from expressions 4.75 and 4.77. Note that for evaluating the discrete part of both expressions (i.e. the summed over terms on the right hand side), we have employed the Mathematica function FindRoot (Wolfram Research, Inc., 2017). The integral in both expressions is solved numerically using the Mathematica function NIntegrate (Wolfram Research, Inc., 2017). labelled as the normal mode approximation (NMA) in figures 4.7a-4.8b. This approximation utilises the observation that the discrete spectrum has two roots s and s, only. Consequently it is possible to obtain an approximate solution for the amplitude a(t) = C exp(st) + C exp( st) with C being a complex constant whose

  4.12b-4.12e). Figure4.10a corresponds to a higher viscosity and thus has a thicker vorticity layer compared to 4.10b. For comparison, we have also

plotted the analytical prediction from 4.56 and 4.77 in 4.11a and 4.11b. A good agreement is observed in the analytically predicted vorticity field vis-a-vis DNS data. For further quantitative comparison, figures 4.12b-4.12e show the vorticity profiles out side the filament (r R 0 + h) at q = p/4 as a function of radius, at

Faraday et al., 1837). The

  this chapter, we study inviscid parametrically forced oscillations or Faraday waves. The study of such oscillations, started with observations by Michael Faraday (

first linearised theoretical analysis of this by Benjamin and Ursell (1954) assumed inviscid, irrotational flow demonstrating that the amplitude of a spatial (Fourier or Bessel) mode is governed by the Mathieu equation. The known properties of Mathieu functions led Benjamin and Ursell (

Table 5 .2: DNS parameters. All stability charts are available in Patankar et al. (2018).

 5

	Linear stability

Test cases for interfacial solvers (in collaboration with Manpreet Singh, M.Tech Thesis, IIT Bombay, 2018)

  The results presented in this chapter are published in "Test cases for two interfacial solvers by Manpreet Singh, PalasKumar Farsoiya and Ratul Dasgupta in Int. J. Multiphase Flow,115, 75-92, 2019. The analytical derivations and Floquet charts have been obtained by Manpreet and are available in his M.Tech. thesis as well as in the supplementary material in Singh

et al. (2019b).

  Implementation of surface tension force always been chink in the armour of two phase interfacial solvers. There are two principal methods (

Brackbill et al., 1992; Popinet and Zaleski, 2002; Tong and Wang, 2007) which

  implement surface tension either as a body force or as a boundary condition respectively. However different flow regimes involving surface tension forces may not always be well described by the available algorithms. Thus there is need for analytical test cases to identify the range of validity of different algorithms (

Popinet, 2018).

  In this chapter, we propose four test cases involving inviscid, capillary, interfacial oscillations on two base state geometries viz. Cartesian and circular. Such interfacial capillary oscillations, both in the free and forced regimes, have many engineering applications and are also of interest in stability studies (

Adou and Tuckerman, 2016; Denner, 2016; Farsoiya et al., 2017; Kumar and Tuckerman, 1994).

  Analytical results for the Cartesian case is available in the literature Explicit formula for the velocity potential and shape of the interface as a function of time, are obtained in each case and compared to inviscid Direct Numerical Simulations (DNS) conducted using two in-house developed codes. One of these is labelled as the free surface solver (FSS hereafter). This ignores the dynamics of the lighter fluid phase (in collaboration with Manpreet Singh, M.Tech Thesis, IIT Bombay, 2018) solving only for the denser phase, explicitly accounting for free surface boundary conditions using the approach presented in Malan et al. (2015). This code was developed by an M.Tech student in the lab Manpreet Singh (Singh et al., 2019b), utilising the VoF and the Navier-Stokes module described earlier in Chapter 2. The details of implementation of the FSS is available at Singh

	(Benjamin and Ursell, 1954; Debnath, 1994) and we further extend these results
	in the present study. Analytical results presented for the circular base state
	(Cauchy-Poisson problem and Faraday waves on a two-dimensional, circular
	drop) are presented in Singh et al. (2019a).

et al. (2019a).

and Ursell (1954); Kumar and Tuckerman (1994) that

  of this as a function of time. Within a linearised approximation, it was shown by Benjamin and Ursell (1954) that a(t) is governed by the Mathieu equation, the stability chart of equation 6.13 on the forcing strength (h) versus wavenumber (k) plane, shows tongue like curves. Choosing a pair of coordinates (h, k) to be inside the tongue implies an unstable, exponentially growing response while choosing parameters outside, leads to a bounded, oscillatory response of the interface. This implies that, at any driving frequency W, a wavenumber k can be de-stabilised through sufficiently high forcing strength Benjamin

	2019b),
	f I (x, y, t) = k 1 cos(kx)	✓	exp [ky] + exp [ ky 2kH 1 ] 1 exp [ 2kH 1 ]	a(t),
	f O (x, y, t) = k 1 cos(kx)	✓	exp [ky] + exp [ ky + 2kH 2 ] 1 exp [2kH 2 ]	a(t).	(6.14)
						(6.15)
	Equation 6.13 has periodic coefficients and can be analysed using Floquet theory
	Benjamin and Ursell (1954); Kumar (1996); Kumar and Tuckerman (1994). It is
	well-known Benjamin
	.. a(t) + k		Tk 2 + (r I r O )h cos(Wt) r I coth(kH 1 ) + r O coth(kH 2 )	a(t) = 0	(6.13)
	where dots on top indicate differentiation. With a(t) known from the solution of
	equation 6.13, the velocity potential in both fluids are obtained as (Singh et al.,

◆ .

◆ .

and Ursell (1954); Kumar

  .1 with symmetry boundary conditions (zero normal derivative for pressure, volume fraction and for the tangential component of velocity, and zero normal component of velocity). As a comparison of the curvature calculation module in the FSS versus OFS, we show results from the static circle test. As discussed earlier in Chapter 2, this test (cf. Tryggvason et al. (2011), sec. 7.3.1 for details) involves demonstrating that at large time, the pressure jump due to surface tension i.e. p I p O = T/R 0 is produced inside the static droplet (circle of radius R

0 in two dimensions as seen in figure 6.3b.). As seen from the contours in figure 6.3a, the pressure inside and outside are seen to be uniform at t = 20 s with a jump of the expected magnitude. Parasitic currents in figure

6

.3c (Tryggvason

  Table 6.2 for We values.

	Case	64 ⇥ 64 FSS OFS	128 ⇥ 128 FSS OFS	256 ⇥ 256 FSS OFS	t
	8	0.11	0.15	0.12	0.13	0.11	0.22	5T 0
	9	0.25	0.15	0.23	0.15	0.22	0.19	5T 0
	10	0.39	0.16	0.37	0.18	0.35	0.18	5T 0
	11	1.37	0.36	1.42	0.31	1.41	0.31	5T 0
	Table 6.3: L • (Dx, Dt) Spatial Convergence
			CFL	FFS	OFS
			0.125 0.117 0.158
			0.25 0.117 0.158
			0.5	0.117 0.158

Table 6 . 4

 64

: L • Temporal convergence for Case 8 (64 ⇥ 64) Dt

Table 6 .5: Computational time(s)

 6

	2	981 1895	1.9
	8	317 1680	5.3
	14	1164 7232	6.2
	18	313 6798	21.7
		Grid -(64 ⇥ 64)

Chirco et al., 2019; Karnakov et al., 2019).

  For capillary flows, spurious currents can arise at long times and at high values of Laplace number which can create numerical instability in flows[START_REF] Popinet | Numerical models of surface tension[END_REF]. Accurate estimation of curvature within a VoF framework can be challenging. Additionally formulating surface tension as a conservative body force is an important challenge (Popinet, 2018). Substantial ongoing research is currently aimed at improving these, leading to better algorithms for surface tension in VoF based methods (We have proposed some new inviscid test cases with linearised analytical solutions in chapter 6. These may be used as a part of a repository of stringent test cases to aid in benchmarking new algorithms. There are research gaps in the analytical theory of waves with Marangoni forces and involving non-newtonian fluids. Here new theories and analytical solutions may be found and used as benchmark tests for the numerical solvers, see Shen et al. (2018) for recent analytical results.

  r U p xB + I(r L kB + r U p xA))

	= a + bI	(I.15)
	x(s) exp[lz] s 2 + sL + w 2 = a cos[ = (a cos[ p xz] b sin[ p xz] + a sin[ p	p xz]) + (b cos[ xz]I + b cos[ p p xz] + a sin[ xz]I b sin[ p p xz])I xz]
	We finally obtain

= ⇢ x(s) exp[lz] s 2 + sL + w 2 = b cos[ p xz] + a sin[ p xz] (I

.16)

I.3 Consistency check for the dispersion relation for hollow gas filament, chapter 4

  A consistency check on 4.38 is to demonstrate that it reduces to the viscous dispersion relation in Cartesian geometry (see Lamb (1932)) in the limit m !•, R 0 ! • such that m/R 0 ! k where k is the wavenumber in Cartesian ((x, z))coordinates (i.e. h(x, t) = a 0 cos(kx)). The following asymptotic formulae (as n ! •) Abramowitz and Segun (1970) are necessary for further analysis,. Consider the limit m ! •, R 0 ! • with m/R 0 = k. Defining k ⌘ kp n s and using I.37, I.38 we obtain the following asymptotic expressions.Using Taylor expansion of ln(1 + e) for e << 1, we obtain from I.41h 2,1

	⇠	✓	1 +	1 k2	◆ 1/2	0 @ 1 ⌥		m	⇣	1 k2 + 1	⌘	1 A + ln	2 6 6 6 4	k +	⇣	k2 + 1	1 ⌥ ⌘ 1/2	1 m ⇢ 1 ⌥	m	⇣	1 k2 +1	⌘	3 7 7 7 5	.
																								(I.41)
						K n (nz) ⇠	⇣ p 2n	⌘ 1/2 exp [ nh] (1 + z 2 ) 1/4		1 + O	✓	1 n	◆	,	(I.37)
						K 0 n (nz) ⇠					⇣ p 2n	⌘ 1/2 1 + z 2 1/4 z	exp [ nh]		1 + O	✓	1 n	◆	,	(I.38)
																			"				#
	with h ⌘ 1 + z 2 1/2 + ln			z 1 + (1 + z 2 )	1/2
	K m 1 ( R0 ) K 0 m 1 ( R0 )	⇠			⇣	1 k2 + 1	⌘ 1/2 + O	✓	1 m	◆	,	(I.39)
	K m+1 ( R0 ) K 0 m 1 ( R0 )	⇠			✓	m 1 m + 1	◆ 1/2	exp	" (	(m + 1) h 2 (m 1) h 1	)#	⇢	1 + z 2 1	2 1 + z 2 z 1	1/4
							+ O	✓	1 m	◆	.					(I.40)
																								2	3
	Here h 2,1 = 1 k ✓ 1 ⌥ ◆ 1 m . Asymptotic expressions for h 2,1 may be obtained as follows. ⇣ 1 + z 2 2,1 ⌘ 1/2 + ln 6 4 z 2,1 1 + ⇣ 1 + z 2 2,1 ⌘ 1/2 7 5 and z 2,1 ⌘ R0 m ± 1	=
	h 2,1 ⇠		1 +	1 k2	✓	1 ⌥	2 m	◆ 1/2	+ ln	2 6 6 6 6 4	1 +	⇢	1 k ✓ 1 + 1 k2 1 ± ✓ 1 ⌥ 1 m	◆ m 2	◆ 1/2	3 7 7 7 7 5

  Re-writing the dispersion relation 4.38 for a hollow jet, in terms of the frequency s and using the asymptotic expressions I.44, we obtains 2 + nk 2

	I.45 may be further written as s 2 " ⇣ k2 + 1 ⌘ 1/2 s 2 " ⇣ k 2 n s + 1 ⌘ 1/2 ( ⇣ k 2 n s + 1 ⌘ 1/2 ( ⇣ k +nk 2 " 1 k ( ⇣ k2 + 1 ⌘ 1/2 k)( ⇣ n )# ⌘ 1/2 s k2 + 1 ⌘ 1/2 k)# 1 + 4 k ⇣ + nk 2 " 1 k n 1/2 ( ⇣ k 2 n s + 1 ⌘ 1/2 k ⇣ n s ⌘ 1/2 )( 1 + 4k ⇣ n s ⌘ 1/2 ⇣ k2 + 1 k 2 n ⌘ 1/2 s + 1 s	⌘ 1/2 )#	s )#	s +	Tk 3 r (I.47) = 0
	+ I.47 may be written as, Tk 3 r Using k = k n s s 2  k 2 n s + 1 k = 0 r 1/2 , we obtain from I.46 n s ⇣ k 2 n s + 1 ⌘ 1/2 + 4k 2 n s + 4 1 4k ⇣ n s ⌘ 1/2 ✓ nk 2 + s + nk 2 s ◆ 1/2 # " ✓ s + nk 2 + s nk 2 Tk 3 r	◆ 1/2 = 0	(I.46) (I.48)
	After cancellations I.48 simplifies to				2
	Using I.42, we obtain at leading order, ⇣ s + 2nk 2 ⌘ 2 + Tk 3 r = 4n 2 k 3	✓	nk 2 + s n	m ◆ 1/2	k	(I.42) (I.49)
	exp	" (	(m + 1) h 2 (m 1) h 1	)#	⇠ 1 + 2	k2 + 2	k ⇣	k2 + 1	⌘ 1/2	(I.43)
	Equations I.39,I.40 and I.43 implies the following leading order behaviour
	K m+1 ( R0 ) K 0 m 1 ( R0 )	⇠	1 + 2	k2 + 2 ⇣ k2 + 1 k ⇣ ⌘ 1/2 k2 + 1	⌘ 1/2	,	G( R0 ) R0 K 0 m 1 ( R0 )	⇠	1 + 2 2 k ⇣ k ⇣ k2 + 1 k2 + 1 ⌘ 1/2 ⌘ 1/2
											(I.44)
			2 6 4	1 + 4 k ⇣ k ⇣ k2 + 1 k2 + 1 ⌘ 1/2 ⌘ 1/2	3 7 5 s +	Tk 3 r	2 6 4	⇣	k2 + 1 ⇣ k2 + 1 ⌘ 1/2 ⌘ 1/2 +	k	3 7 5 = 0	(I.45)

which is the dispersion relation provided in Lamb (1932).

  1 2 [K m 1 (z) + K m+1 (z)], we have K m 1 (z) z KFor small viscosity n, we pose an expansions = s 0 + n 1/2 s 1 + ns 2 + . . .Solving these we obtain a complex conjugate pair of series with the following coefficients s 0 = ±Iw 0 , s 1 = 0, s 2 =

	obtain	0 m 1 (z)	⇠	1 z		and	z K m+1 (z) K 0 m 1 (z)	⇠ (z + 2m). Using these in I.51, we
	s 2 +	2m(m + 1)n R 2 0	✓	1 +	1 2m	( R0 + 2m) +	m 1 R0	◆	s + w 2 0	✓	1 +	m + 1 R0	◆	= 0
																					(I.52)
	which maybe written in dimensional form as,
	s 2 +	2m(m + 1)n R 2 0	✓	2 +	R 0 2m	r	s n	+	m 1 R 0	r	n s	◆	s + w 2 0	✓	1 +	m + 1 R 0	r	n s	◆	= 0
																					(I.53)
																					(I.54)
	and obtain at various orders
	O(1) : O(n 1/2 ) : O(n) :	s 2 0 + w 2 0 = 0, 2s 0 s 1 + m + 1 R 0 2s 0 s 2 + 4m (m + 1) s 0 s 3/2 0 + w 2 0 (m + 1) R 0 R 2 0 + s 2 1 + 3 (m + 1) s 1/2 0 = 0, 2R 0 s 1 s 1/2 0	✓	m + 1 2R 0	◆	w 2 0 s 1 s 3/2 0	= 0
																				2m(m + 1) 0 R 2	. . .	(I.55)

we use the word filament and jet interchangeably, the only difference being the absence or presence of an axial velocity profile respectively, assumed here to be uniform (Chandrasekhar, 1981; Eggers and Villermaux,

2008),99

<(z) and =(z) denote the real and imaginary part of z

Acknowledgements

I begin by thanking my supervisor Professor Ratul Dasgupta whose guidance and support have enabled me to enjoy the pleasures of research. I thank the members of my research progress committee, Professors Mahesh Tirumkudulu, Partha S. Goswami and Prof. Atul Sharma for their feedbacks and suggestions in annual progress seminars. Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Case

Comments the OFS agrees better with analytical predictions. Finally we present the effect of larger values of e and e az at high density ratio r r = 10 3 in figures 6.14a and 6.14b. Like the case of free oscillations, we find here as well that for larger values of e (e az ) the FSS shows better agreement than the OFS, to analytical predictions even when the density ratio is quite high. The OFS is found to display higher numerical dissipation in both circular and planar cases in this parameter range.

Grid

Summary and Conclusions

A brief point-wise summary of significant results obtained in this thesis are summarised below.

• In Chapter 1, we review the literature on viscous, free & forced interfacial waves highlighting the lacunae that exist. This description defines the problems studied in this thesis.

• In Chapter 2, we discuss the numerical implementation and benchmarking of in-house solver for solving the incompressible, two dimensional Navier-Stokes equations with an interface.

• In Chapter 3, the problem of axisymmetric interfacial capillary-gravity waves on a cylindrical deep pool is analysed using the initial value problem approach. We show the presence of continuous spectrum in this geometry.

We also demonstrate the energy transfer to higher modes in the large amplitude regime. Analytical results obtained are verified by our in-house developed solver.

• In Chapter 4, capillary waves on an azimuthally perturbed gaseous filament surrounded by viscous fluid is studied using normal modes as well as the initial value problem approach. We demonstrate the existence of the continuous spectrum using careful normal mode analysis and follow it up with solution to the initial value problem. Our analytical results are verified by the direct numerical simulations using our in-house solver.