Early- and late-time evolution of Rayleigh-Taylor instability in a finite-sized domain by means of group theory analysisAnnie Naveh

TL;DR

This paper presents a group theory-based analysis of the late-time evolution of Rayleigh-Taylor instability in finite domains, revealing how domain size influences flow dynamics and solution multiplicity.

Contribution

It introduces a systematic theoretical framework using group theory to analyze nonlinear Rayleigh-Taylor instability in finite-sized domains, including asymptotic solutions and domain size effects.

Findings

Flow in finite domains is slower than in infinite domains.

Domain size affects the flow's diagnostic parameters.

Inter-facial shear function links to solution multiplicity.

Abstract

We have developed a theoretical analysis to systematically study the late-time evolution of the Rayleigh-Taylor instability in a finite-sized spatial domain. The nonlinear dynamics of fluids with similar and contrasting densities are considered for two-dimensional flows driven by sustained acceleration. The flows are periodic in the plane normal to the direction of acceleration and have no external mass sources. Group theory analysis is applied to accurately account for the mode coupling. Asymptotic nonlinear solutions are found to describe the inter-facial dynamics far from and near the boundaries. The influence of the size of the domain on the diagnostic parameters of the flow is identified. In particular, it is shown that in a finite-sized domain the flow is slower compared to the spatially extended case. The direct link between the multiplicity of solutions and the inter-facial shear function is explored. It is suggested that the inter-facial shear function acts as a natural parameter to the family of analytic solutions.

Equations

$$v_A=\sqrt{\frac{g}{3k}}\text{and}{\zeta }_A=-\frac{k}{6}$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0$$

$$\frac{\partial \mathbf{v}}{\partial t}+\left(\mathbf{v}\cdot \nabla \right)\mathbf{v}+\frac{\nabla P}{\rho }=0$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left(E+P\right)\mathbf{v}=0$$

$$\left(\rho ,\mathbf{v},P,E\right)={\left(\rho ,\mathbf{v},P,E\right)}_{h}H(\theta )+{\left(\rho ,\mathbf{v},P,E\right)}_{l}H(-\theta )$$

$$\rho \left(\frac{1}{\nabla \theta }\frac{\partial \theta }{\partial t}+\mathbf{v}\cdot \mathbf{n}\right)=0$$

$$\big{[}\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}\big{]}=0\quad,\quad\big{[}\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\tau}\big{]}=\textrm{arbitrary}\quad,\quad[P]=0\quad,\quad[W]=\textrm{arbitrary}$$

$${\mathbf{v}}_h=0\text{at}z=Z,{\mathbf{v}}_l=0\text{at}z=-Z$$

$$\mathbf{v}=\nabla \Phi +\nabla \times \mathbf{\varphi }$$

$$\mathbf{u}=\nabla \Phi$$

$${\nabla }^2\Phi =0$$

$$\rho \left(\frac{\partial \Phi }{\partial t}+\frac{\nabla {\Phi }^2}{2}\right)+P=0$$

$${\rho }_h\left(\nabla {\Phi }_h\cdot \mathbf{n}+\frac{\dot{\theta }}{|\nabla \theta |}\right)={\rho }_l\left(\nabla {\Phi }_l\cdot \mathbf{n}+\frac{\dot{\theta }}{|\nabla \theta |}\right)=0$$

$$\nabla {\Phi }_h\cdot \mathbf{\tau }-\nabla {\Phi }_l\cdot \mathbf{\tau }=\text{arbitrary}$$

$${\rho }_h\left[\frac{\partial {\Phi }_h}{\partial t}+\frac{|\nabla {\Phi }_h{|}^2}{2}+\left(g+\frac{dv}{dt}\right)z\right]={\rho }_l\left[\frac{\partial {\Phi }_l}{\partial t}+\frac{|\nabla {\Phi }_l{|}^2}{2}+\left(g+\frac{dv}{dt}\right)z\right]$$

$${\frac{\partial {\Phi }_h}{\partial z}|}_{z=Z}=-v(t),{\frac{\partial {\Phi }_l}{\partial z}|}_{z=-Z}=-v(t)$$

$${\Phi }_h=\sum _{m=1}^{\infty }{\Phi }_m(t)\left[z\sinh (mkZ)+\frac{\cos (mkx)}{mk}\cosh \left[mk(z-Z)\right]\right]+f_h(t)$$

$${\Phi }_l=\sum _{m=1}^{\infty }{\tilde{\Phi }}_m(t)\left[-z\sinh (mkZ)+\frac{\cos (mkx)}{mk}\cosh \left[mk(z+Z)\right]\right]+f_l(t)$$

$$z^{\ast }(x,t)=\sum _{i=0}^N{\zeta }_i(t)x^{2i}$$

$$\begin{array}{r}{M}_n=\sum _{m=1}^{\infty }{\Phi }_m(t)(mk{)}^n\sinh (mkZ)\\ {\tilde{M}}_n=\sum _{m=1}^{\infty }{\tilde{\Phi }}_m(t)(mk{)}^n\sinh (mkZ)\\ N_n=\sum _{m=1}^{\infty }{\Phi }_m(t)(mk{)}^n\cosh (mkZ)\\ {\tilde{N}}_n=\sum _{m=1}^{\infty }{\tilde{\Phi }}_m(t)(mk{)}^n\cosh (mkZ)\end{array}$$

$$\dot{{\zeta }_1}-3N_1{\zeta }_1-\frac{M_2}{2}=\dot{{\zeta }_1}-3{\tilde{N}}_1{\zeta }_1+\frac{{\tilde{M}}_2}{2}=0$$

$$(1+A)\left(\frac{\dot{N_1}}{2}+{\zeta }_1\dot{M_0}-\frac{N_1^2}{2}-g{\zeta }_1\right)=(1-A)\left(\frac{\dot{\widetilde{N_1}}}{2}-{\zeta }_1\dot{\widetilde{M_0}}-\frac{{\tilde{N}}_1^2}{2}-g{\zeta }_1\right)$$

$$M_0=-v(t),{\tilde{M}}_0=v(t),M_0=-{\tilde{M}}_0$$

$$N_1-{\tilde{N}}_1=\text{arbitrary}$$

$${\Phi }_m(t),{\tilde{\Phi }}_m=0\forall m>1$$

$$\dot{{\zeta }_1}+\frac{v{k}^2}{2}=0$$

$$\dot{v}k+2Ag{\zeta }_1\tanh (kZ)=0$$

$$v(t)\sim \exp \left(t\sqrt{\tanh (kZ)}/\tau \right)$$

$$|{\zeta }_1(t)|\sim \exp \left(t\sqrt{\tanh (kZ})/\tau \right)$$

$$v=-\frac{12\sqrt{-Ag{\zeta }_1}{\sinh }^2(kZ)\left[2{\zeta }_1\sinh (kZ)+k\cosh (kZ)\right]}{k^2\left[3\sinh (kZ)+\sinh (3kZ)\right]\sqrt{\frac{A\left(k^2-4{\zeta }_1^2\right)+A\left(4{\zeta }_1^2+k^2\right)\cosh (2kZ)-4{\zeta }_{1}k\sinh (2kZ)}{{\left[k\cosh (kZ)-2{\zeta }_1\sinh (kZ)\right]}^2}}}$$

Full text

Early- and late-time evolution of Rayleigh-Taylor instability in a finite-sized domain by means of group theory analysis

Annie Naveh

Miccal T. Matthews

Snezhana I. Abarzhi Corresponding author, email snezhana.abarzhi@gmail.com School of Mathematics and Statistics, The University of Western Australia, Perth, Western Australia 6009, Australia

We have developed a theoretical analysis to systematically study the late-time evolution of the Rayleigh-Taylor instability in a finite-sized spatial domain. The nonlinear dynamics of fluids with similar and contrasting densities are considered for two-dimensional flows driven by sustained acceleration. The flows are periodic in the plane normal to the direction of acceleration and have no external mass sources. Group theory analysis is applied to accurately account for the mode coupling. Asymptotic nonlinear solutions are found to describe the inter-facial dynamics far from and near the boundaries. The influence of the size of the domain on the diagnostic parameters of the flow is identified. In particular, it is shown that in a finite-sized domain the flow is slower compared to the spatially extended case. The direct link between the multiplicity of solutions and the inter-facial shear function is explored. It is suggested that the inter-facial shear function acts as a natural parameter to the family of analytic solutions.

1 Introduction

In 1883 Lord Rayleigh questioned the outcome of supporting a denser layer of fluid on top of a lighter fluid layer under the effect of gravity, and with this seemly straightforward idea the challenge of understanding the Rayleigh-Taylor instability was born [22]. An instability is generally identified as a component within a system that grows without bound. The Rayleigh-Taylor instability (RTI) is a well-known example of an instability in fluid dynamics. It occurs in fluids, plasmas or (in extreme circumstances) materials. RTI develops on the interface of two fluids (or plasmas or solids) of different densities as they are accelerated against their density gradients; that is, the acceleration is directed from the heavy to the light fluid. While the fluid interface is kept perfectly planar and normal to the external acceleration it is in an unstable state of equilibrium and the instability will not form. However, even a small disturbance of this interface is sufficient to initiate the unbounded growth of RTI. For example we observe RTI whenever a denser, heavier fluid such as water is placed on top of a lighter fluid such as oil and the external acceleration in this case is due to gravity. RTI is initiated as potential energy is released and finger-like protrusions are observed as large coherent bubble structures of the lighter fluid penetrates the heavier fluid and spikes of the heavy fluid penetrates the light. A shear function across the interface results in vortical small scale structures (Kelvin-Helmholtz instabilities) on the side of the evolving spikes. This periodic array of bubbles and spikes as the two fluids mix in time is what defines RTI. An exact reliable description of this evolution is the intellectually rich challenge set by Lord Rayleigh decades ago.

The challenge to understand the growth of RTI remains relevant today for its role in natural phenomena ranging from astrophysical to microscales; in inertial confinement fusion, laser ablation, combustion, atmospheric flows, as well as in supernovae explosion and the formation of the universe [9][21][23][15]. Further the ability to understand and control this instability is vital for industrial applications such as in free-space optical telecommunication, laser micromachining and in aeronautics [18][10].

Previous theoretical studies have considered RTI in an infinitely extended spatial domain in the direction of acceleration. More often however, real world physical phenomena and applications in industry exist in a bounded space. This creates a need to translate these infinitely extended mathematical models to more useful physical models that may be applied to real world systems. An important question must be asked: how does finite spatial boundaries influence the characteristic evolution of the instabilities? Answering this question not only reveals the particular effects of imposing a finite domain on the evolution of the instabilities, but also the convergence to the infinitely extended case.

1.1 Review of approaches

Lord Rayleigh [22] was the first to theoretically study RTI in the incompressible, inviscid fluid case. In doing so, he distinguished himself as not only the first to recognise the significance of hydrodynamic instabilities, but also for his contributions to the theoretical understanding of the initial evolution. Rayleigh found that RTI initially develops faster for smaller values of the wavelength of the initial perturbation, $\displaystyle\lambda$.

Since Rayleigh’s paper, our understanding of the initial stages of RTI has been enhanced significantly and the idealisations of the inviscid, incompressible fluid overcome. The early time growth rate of the instability has been found for more involved cases taking into account effects such as surface tension, viscosity and compressibility in non-ideal fluids [7]. However, finding a description of the instability evolution in the later, nonlinear time regime presents a more formidable challenge explicit in the long history of the development of solution methods.

1.1.1 First estimations

G.I. Taylor [24] was the first to observe the nonlinear evolution of the instabilities experimentally, working with fluids of contrasting densities (specifically water or ethanol and air) in cylindrical tubes. In his work the physical importance of RTI in accelerated fluid layers could first be demonstrated. Taylor considered three-dimensional cylindrical symmetric dynamics in the approximation of a potential flow and driven by a constant acceleration. Taylor numerically estimated a solution to the steady state velocity of RTI and used the tube radius $\displaystyle\lambda/2$ and gravity $\displaystyle g$ as the scales of the problem. He found the velocity of the bubble to be $\displaystyle v\approx 0.49\sqrt{g(\lambda/2)}$ which corresponded with his experimental results. Indeed he lends his name to the instability for his contributions and predominantly for his cautious experimental observations.

1.1.2 Layzer-approach and a unique solution

A major advancement towards a trustworthy description of the nonlinear instability evolution was from Layzer in 1955 [17]. Layzer assumed the flow of a single fluid system to be a potential flow. For incompressible flows a velocity potential satisfies the Laplace equation. Accounting for the conservation laws of mass and momentum, Layzer expanded the potential flow using Bessel functions (for the three-dimensional case) and Fourier series (for the two-dimensional case). He retained only the lowest-order terms (first harmonic) and expanded the conservation laws in the vicinity of the bubble tip. An important difference of the Layzer approach to the method presented in this work, as part of a group theory approach, is that only a single Fourier harmonic is retained. Although the asymptotic solutions of the Layzer approach was in good agreement with the results of Taylor, it does not account for the interplay of multiple harmonic modes and fails when the harmonic exceeds unity. A multiple harmonic description must be considered as it influences the evolution of the diagnostic parameters.

1.1.3 Parameter family of solutions

Garabedian (1957) developed a much deeper analysis that considered multiple harmonics using conformal mapping. He was the first to suggest that there is not one singular steady state solution but a family of solutions governed by a parameter [12]. This implies that nonlinear solutions are in fact not unique but belong to a finite family bounded by the symmetries of the system. He identified the physically significant solution in the family as the fastest one. This solution, however, differed from the Layzers solution by approximately 10 %.

1.1.4 Moments equations

In 1992 Inogamov advanced the idea of a one parameter family of solutions and introduced moment equations (weighted sums of Fourier amplitudes) to solve the nonlinear system to higher orders. He used potential flow and stream functions to reach a one parameter family of solutions [16].

1.1.5 Group theory

Results from a group theory approach, first applied to the RTI problem by Abarzhi in 1995 [6], ties together and makes sense of the previous results of Layzer (1955), Garabedian (1957) and Inogamov (1992) (all studied for a single fluid system). Using group theory, Abarzhi showed that in the two-dimensional case there is a one parameter family of solutions. In the first-order approximation the fastest solution in the family has a bubble curvature and velocity in agreement with Layzer’s solution [6]. In higher orders the velocity value of the fastest solution increases and agrees with Garabedian’s solution. The discrepancy between Layzer and Garabedian’s results is due to an order of approximation. Further, there is no Layzer-type solution, only a Layzer first order approximation. The necessity of considering higher-order interactions for obtaining a reliable description is one of the indications of the non-local character of the instability evolution.

Finally, the group theory approach finds a stability region for the curvature of the bubble tip (the family parameter) and the corresponding velocity [5]. Inogamov’s solution sits outside of this stability region. The results obtained from a group theory approach are confirmed in advanced numerical simulations and in more accurate experiments [25] [13].

1.2 Previous results of group theory for an infinite domain

An analysis based on group theory is applicable for many RTI-related problems in two and three dimensions with various symmetries for $\displaystyle 0<A\leq 1$, where $\displaystyle A$ is the system’s Atwood number [19] [20] [14] [8]. Local properties of the evolution are accounted for in local spatial expansions and global properties are accounted for in Fourier series expansions [11]. A group theory approach reliably accounts for the asymptotic dynamics for nonlinear RTI bubbles in an infinite spatial domain. A stability analysis suggests that the obtained RTI families are complete as they involve all the possible solutions allowed by the spatial symmetry of the flow [5]. For asymptotic time $\displaystyle t/\tau\rightarrow\infty$ (where $\displaystyle\tau$ is a suitable time scale) to the first-order RTI dynamics approaches a steady-state bubble velocity and curvature [2].

For contrasting fluid densities, $\displaystyle A=1$, and at the first-order of approximation the fastest solution in the parameter family of solutions corresponds with a Layzer-approach approximation. This yields a maximum velocity and corresponding curvature of

[TABLE]

where $\displaystyle k$ is the wave number of the Fourier expansion. The critical bubble curvature, $\displaystyle\zeta_{cr}$, is defined as the maximum curvature supported by the system. Independent of Atwood number it is found in an infinite spatial domain as $\displaystyle\zeta_{cr}=-k/2$ and a bubble curvature outside of the region $\displaystyle 0\leq|\zeta_{1}|\leq|\zeta_{cr}|$ is not supported by the convergence of the solutions.

2 Theoretical modelling approach

2.1 System description

Fluid flow problems in general are very difficult and often impossible to solve analytically without appropriate assumptions set in place. This is due to the equations governing the fluid flow; namely the celebrated Euler or (compressible) Navier-Stokes equations. The set of equations span three dimensions, are time dependent, have pressure gradient terms and include viscous stress. Indeed the specific fluid problem of modelling the nonlinear evolution of acceleration-driven hydrodynamic instabilities is undeniably complex [5]. As in many fluid flow descriptions, to successfully navigate this problem and obtain a rigorous description of the diagnostic parameters certain theoretical approximations must be set in place.

2.1.1 Fluid idealisation

Firstly, in this work the effects of compressibility, viscous stress and surface tension are neglected. This approximation proves important as finding the proper choice of boundary conditions for compressible or miscible fluids is a very difficult and, to a large extent, an unsolved problem [5].

2.1.2 Two-dimensional

The flow is assumed to be two-dimensional; periodic in the $\displaystyle x-$direction, motionless in the $\displaystyle y-$direction and spatially bounded in the $\displaystyle z-$direction. The latter is set as the direction of acceleration. Fixed boundaries are imposed for all values of $\displaystyle x$ at positions $\displaystyle z=\pm Z$. The set-up of the problem is illustrated in Figure 1. The extended periodicity of the coherent structures in the $\displaystyle x-$direction is an important symmetry property enabling a group theory consideration. It is for this reason that the finite boundaries are imposed only in the direction of acceleration and the periodic structures repeat infinitely in the $\displaystyle x-$direction. To expedite the calculations, an analysis is performed in the non-inertial frame of reference moving with velocity $\displaystyle v(t)$ in the $\displaystyle z-$direction, where $\displaystyle v(t)$ is the bubble velocity at the bubble tip in the laboratory frame of reference. Solutions are translated back into the laboratory reference frame where necessary.

2.1.3 Foundational scales

The evolution of the instability should be found such that it is true at all scales, requiring an identification of the appropriate scales of the problem. The spatial period of the initial perturbation $\displaystyle\lambda$, $\displaystyle[\lambda]=m$ and the uniform acceleration $\displaystyle g$, $\displaystyle[g]=ms^{-2}$ are the foundational scales for RTI. The length scale may be represented as the wave number $\displaystyle k$ where $\displaystyle k=2\pi/\lambda$ and $\displaystyle[k]=m^{-1}$. These basic scales define the time scale of the flow $\displaystyle\tau\sim\sqrt{gk}^{-1}$, the initial bubble growth rate $\displaystyle v_{0}\sim\sqrt{g/k}$ and the frequency $\displaystyle\omega\sim\sqrt{gk}$. For a realistic system of non-ideal fluids the characteristic length scale $\displaystyle\lambda$ is set by the mode of fastest growth, whereas for this analysis, regarding ideal fluids, the spatial period $\displaystyle\lambda$ is set by the initial perturbation.

2.1.4 Scale separation

For a large density ratio of the two fluids the nonlinear regime dynamics of the RT flow can be categorised into two relatively independent scales [1]. The large scale set by the initial perturbation wavelength ($\displaystyle\sim\lambda$) includes the coherent structures of bubbles and spikes, while the small-scale set by the amplitude of the initial perturbation ($\displaystyle\ll\lambda$) consists of the vortical structures caused by the shear function. It should be noted that this theoretical separation of scales is not applicable for fluids with very close densities (that is, for small Atwood numbers $\displaystyle A$). This is because as the difference in fluid densities approach zero the vortical structures become large and the scale separation is no longer valid. The small scale vortical structures are a consequence of the non-linearities and secondary instabilities from the full Navier-Stokes equations [3]. These small scale dynamics and the interaction between the scales results in a randomness that convolutes the mixing as a complex process. Fortunately, the multi-scale mixing process does maintain certain features of coherence and order associated primarily with the dynamics of the large scales [14]. To theoretically solve this problem we therefore focus our solutions to the large-scale coherent dynamics; the bubble and spike structures.

2.1.5 Localised to the bubble-tip

The mushroom-type spike gets its shape from the non-deterministic small-scale vortical structures. Further to describe RTI evolution in the vicinity of the spike requires a very high order of approximation to deal with the many singularities [14]. Therefore to obtain a trustworthy deterministic description of the unstable dynamics in the large scales we find solutions in the vicinity of the coherent, large scale bubble front.

2.1.6 Time regimes

To date an accurate theoretical description of the evolution of RTI for all times has not been found [5]. Instead studies of these instabilities (and the work presented here) are predominantly concerned with the evolution in early time $\displaystyle t\ll\tau$ (the linear regime) and latter time in which the bubble and spike structures are formed (the nonlinear regime). For a description in the nonlinear regime the dynamical system is solved for asymptotically large time $\displaystyle t/\tau\rightarrow\infty$. The evolution in this regime reveals the true intellectual richness of this problem and so is the main focus of this work.

2.1.7 Multiple harmonics

A multiple harmonic analysis of the system is considered such that waves that are integer multiples of the initial perturbation wavelength are retained in our analysis. This is crucial as these higher-order harmonics contribute to the diagnostic parameters of the motion and should be accounted for in the description of the evolution [16]. The major assumption of the multiple harmonic analysis is that the dynamics of the flow are governed by a dominant mode. For a realistic system of non-ideal fluids the dominant mode and characteristic length scale is set by the mode of fastest growth [2]. On the other hand, for this analysis regarding ideal fluids, the dominant mode is set by the initial perturbation for the RTI.

2.2 Governing equations

The dynamics of the idealised fluid system considered is governed by the set of incompressible Euler equations; a simplification of the more general Navier-Stokes equations with zero viscosity and a constant density imposed. The set of incompressible Euler equations depict, respectively, the conservation of mass, momentum and energy:

[TABLE]

[TABLE]

[TABLE]

Here $\displaystyle\rho,\boldsymbol{v},P,E$ is the field of density, velocity, pressure and energy and $\displaystyle t$ is time. The energy field can be further expressed as $\displaystyle E=\rho\left(e+\frac{\boldsymbol{v}^{2}}{2}\right)$ with $\displaystyle e$ being a specific internal energy. All the variables are considered to be continuous functions of the spatial coordinates and time. For an ideal incompressible fluid the density and the internal energy of the fluid are constant in the fluid’s bulk.

2.2.1 Inter-facial boundary conditions

For incompressible, immiscible fluids, the fluid interface of the two fluids is discontinuous and this motivates the introduction of a local scalar function $\displaystyle\theta(x,y,z)$ defined such that $\displaystyle\theta=0$ at the interface, the heavy (denser) fluid is located at $\displaystyle\theta>0$ and the light (less denser) fluid at $\displaystyle\theta<0$. It is assumed that the derivatives $\displaystyle\dot{\theta}$ and $\displaystyle\nabla\theta$ exist. Using this simple function the system can be explicitly expressed as

[TABLE]

where $\displaystyle H(x)$ is the Heaviside function and the subscripts $\displaystyle h$ and $\displaystyle l$ refer to the observables in the heavy and light fluid, respectively. Describing the incompressible two fluid system holistically using the $\displaystyle\theta$ function and substituting it into the governing Euler equations (2.2) yields the following inter-facial boundary conditions in the case of zero mass flux across the interface

[TABLE]

[TABLE]

The parenthesis $\displaystyle[\dots]$ denotes a jump of functions across the interface, and $\displaystyle W=e+\frac{P}{\rho}$ is the specific enthalpy. The inter-facial boundary conditions indicate, respectively, that across the interface; mass flux is continuous, normal component of velocity of the fluid is continuous, tangential component of velocity is not continuous and pressure is continuous. As the fluids are considered incompressible we can omit the discontinuity of specific enthalpy $\displaystyle W$ from our consideration. Extra attention should be given to the boundary condition indicating the discontinuity of the tangential component of velocity as this allows a shear function to develop along the interface. Later it is shown that this shear function may act as the natural parameter to the family of solutions.

The boundary conditions at the outside boundaries are simply

[TABLE]

2.2.2 Large-scale dynamics

Any vector field can be written as the sum of a scalar potential function and a curl and so the velocity vector fluid of the fluids can be expressed as

[TABLE]

The large-scale dynamics are assumed to be irrotational therefore $\displaystyle\nabla\times\boldsymbol{\phi}=0$. Indeed in the small scale the flow is not irrotational (due to the interfacial vortical structures). Considering only the large scales and expressing the velocity as a scalar potential is a valid and necessary approximation allowing us to reach a solution analytically. For large scales, therefore, the fluid velocity field is conservative and may be presented by a gradient potential function

[TABLE]

for some scalar potential $\displaystyle\Phi$. Limiting our analysis to the large-scale dynamics, the governing Euler equations (2.2) and derived inter-facial boundary conditions (2.2.1) can be expressed in terms of the scalar potential $\displaystyle\Phi$.

Substituting (7) into the conservation equations (2.2) reduce to

[TABLE]

[TABLE]

Notice that this system is governed by the Laplace equation; this is as expected as the flow is both irrotational and incompressible. To complement the governing conservation equations at large scales the inter-facial boundary conditions (2.2.1) should likewise be expressed in terms of the scalar potential $\displaystyle\Phi$. Substituting (7) into (2.2.1) yields the appropriate inter-facial boundary conditions at large scales given below in the non-inertial reference frame of the bubble-tip:

[TABLE]

[TABLE]

[TABLE]

The set of defining equations is completed by defining the outside boundary conditions for the instability in a finite domain

[TABLE]

2.3 Method of solution

2.3.1 Symmetry groups

Symmetry in mathematics implies invariance of a given system under certain transformations. The periodicity of the large scale coherent motion in the $\displaystyle x-$direction is a natural property of two-dimensional RTI, and is used as the defining symmetry element. This natural periodic structure is invariant with respect to one of the seven one-dimensional crystallographic symmetry groups. Explicitly these seven groups are $\displaystyle p1$, $\displaystyle p1m1$, $\displaystyle p11\texttt{g}$, $\displaystyle p2$, $\displaystyle p2m\texttt{g}$, $\displaystyle p11m$, and $\displaystyle p2mm$ using international classification and Fedorov’s notation. However, not all of the seven invariant groups should be considered as an appropriate portrayal of RTI. This is because the instability flows are essentially anisotropic and the dynamics in the $\displaystyle z-$direction of acceleration differs from that in the other directions. An additional requirement for the correct symmetry group is one with coherent structures that are observable and repeatable – we need to find a group such that the symmetry properties do not change over time and therefore are structurally stable. This is satisfied if the group is a symmorphic group with inversion in the plane [4]. Group $\displaystyle p1m1$, hereafter $\displaystyle pm1$, generators are translation $\displaystyle x+\lambda\rightarrow x$ and mirror reflection in the $\displaystyle(z,y)$ plane $\displaystyle x\rightarrow-x$ and is the symmetry group most suitable for the analysis of the coherent structures of RTI growth.

2.3.2 Fourier series expansion

Using a group theory approach with a strong consideration for the symmetries in the system inspires the application of a Fourier Series in solving the governing Laplace equation. In the case of continuous translational symmetry $\displaystyle x+\lambda\rightarrow x$ in which the translation parameter $\displaystyle\lambda$ can take any value, the set $\displaystyle\sin{(kx)}$ forms a complete set of irreducible representations of the odd functions and $\displaystyle\cos{(kx)}$ for the even functions. The Fourier series are hence irreducible representations of the group of translations. For the $\displaystyle pm1$ group the fluid potentials are expressed as

[TABLE]

[TABLE]

where $\displaystyle\Phi_{m}(t)$ and $\displaystyle\tilde{\Phi}_{m}(t)$ are, respectively, the Fourier amplitudes of the heavy and light fluid corresponding with the mode $\displaystyle m$, and $\displaystyle f_{h}(t)$ and $\displaystyle f_{l}(t)$ are time-dependent functions.

2.3.3 Spatial expansion

In the vicinity of the bubble tip the interface, set as $\displaystyle z^{*}$, evolution must take into account the $\displaystyle pm1$ generator $\displaystyle x\rightarrow-x$, simplifying the Taylor expansion of the bubble tip to

[TABLE]

where $\displaystyle N$ is the order of approximation. In this work the case of $\displaystyle N=1$ will be examined due to the complexity of the problem. In the frame of reference of the bubble tip $\displaystyle\zeta_{0}=0$ and so to the leading order ($\displaystyle N=1$); $\displaystyle z^{*}(x,t)=\zeta_{1}(t)x^{2}$. As the bubble tip takes a predominantly quadratic form $\displaystyle\zeta_{1}(t)$ is the principal curvature of the bubble tip and as the bubble is concave we look for solutions corresponding to $\displaystyle\zeta_{1}(t)<0$.

For this expansion to be valid we assume the following conditions $\displaystyle\frac{x}{\lambda}<<1$ and $\displaystyle\frac{z-z^{*}}{Z}<<1$ or equivalently for the last constraint $\displaystyle z-z^{*}<<Z$. Hence the expansions and following solutions hold only for a sufficiently large domain. Solutions are being considered for a finite but large domain.

2.3.4 Moments

The weighted sums of the infinite number of Fourier amplitudes are named the moments and expressed as

[TABLE]

where $\displaystyle n=0,1,2\dots$

Using these moment expressions the dynamical system in $\displaystyle N=1$ is derived for the first time to give the conditions for the continuity of mass flux, normal component of velocity and pressure in a finite domain at the interface and at the outside boundaries of the finite domain as

[TABLE]

[TABLE]

[TABLE]

and the condition for the discontinuity of the tangential component at the interface is

[TABLE]

Recall that $\displaystyle A$ is the Atwood number of the two fluid system and defined as $\displaystyle A=\frac{\rho_{h}-\rho_{l}}{\rho_{h}+\rho_{l}}$ and a dot marks a time-derivative.

Expressing the system in terms of infinite sums provides the principal opportunity to derive the regular asymptotic solutions for $\displaystyle t/\tau\rightarrow\infty$ with a desired accuracy that accounts for the higher-order spatial modes and enables a stability analysis of the asymptotic dynamics. To find the solutions for the dynamical system in the linear and nonlinear regimes, we have to solve the closure problem. The closure problem can be addressed using the arguments of symmetry; specifically by considering all the local asymptotic solutions allowed by the spatial symmetry of the flow. For the discrete group $\displaystyle pm1$ the number of family parameters $\displaystyle N_{p}=1$. The expressions for the moments in the defining equations should be expanded such that in every order of approximation $\displaystyle N$, the number of equations $\displaystyle N_{e}$ and the number of variables $\displaystyle N_{t}$ (the Fourier amplitudes $\displaystyle\Phi_{m},\tilde{\Phi}_{m}$ and surface variables $\displaystyle\zeta_{i}$) obey the relation $\displaystyle N_{t}\geq N_{e}+N_{p}$. In this case, $\displaystyle N_{p}=1$ and $\displaystyle N_{e}=4$ indicating that $\displaystyle N_{t}\geq 5$. To find the nonlinear asymptotic solution the proper order of approximation of the moments should be taken to solve the closure problem and establish the proper relations between the moments.

3 Results

The diagnostic parameters of RTI i.e the growth rate and curvature of the bubble tip, can be characterised for early and asymptotically late time. These time regimes are named, respectfully, the linear and nonlinear regime. The strength of the outlined methodology is that it is applicable in both these time regimes and for systems with a range of Atwood numbers. For consistency solutions are given for systems of all Atwood numbers $\displaystyle 0<A\leq 1$. It should be stressed, however, that the potential separation from which the equations are derived are only valid for large Atwood numbers and so we should proceed cautiously when approaching $\displaystyle A\rightarrow 0$. Further these solutions hold only for a finite but large domain size. Although for completion the solutions are given in a form that appears to hold for all boundary heights including $\displaystyle kZ\rightarrow 0$, in this small boundary limit the results should be received tentatively.

Prior to discussing the specific effects of imposing a finite domain on the growth of the instability it is useful to explore the convergence of the finite domain to the infinitely extended case as $\displaystyle kZ\rightarrow\infty$. This convergence is found to be independent of Atwood number and for any fixed curvature a critical scaled boundary height of $\displaystyle kZ=5$ is established. At this domain size the percentage difference between the diagnostic parameters in the finite domain and those in the infinite domain falls below a negligible 0.02 %. In other words, at a boundary height of $\displaystyle kZ=5$ the effects of the finite boundaries on the instability become negligible and therefore is a suitable approximation for an infinite domain.

3.1 Linear regime – unique solution

For early time $\displaystyle t/\tau\ll 1$, namely the linear regime, the perturbation amplitude is small and all harmonics above the first order are disregarded

[TABLE]

as it is assumed that for early time the multiple modes have not had time to form. In the linear time regime we work under the approximation of an initially small bubble curvature and perturbation amplitude, $\displaystyle 0<|\zeta_{1}(t_{0})k|<<1$ and $\displaystyle 0<|(k/g)^{1/2}v(t_{0})|<<1$ where $\displaystyle t_{0}$ represents the initial instant of time and $\displaystyle v(t_{0})$ the initial growth rate. For early time $\displaystyle t\approx t_{0}$ and $\displaystyle v\approx|v(t_{0})|$, therefore enabling a trustworthy simplification of the governing dynamical system to the first-order in $\displaystyle v(t_{0})$ and $\displaystyle t_{0}$ by assuming that the contribution of higher orders is negligible. The governing system (2.3.4), retained to the first order in $\displaystyle v(t)$ and $\displaystyle t$ is reduced to a set of simplified coupled ordinary differential equations

[TABLE]

[TABLE]

Solving these equations yields the growth and curvature rate of the initial RTI perturbation at the fluid interface. We find that for RTI bounded at $\displaystyle Z$ and $\displaystyle-Z$ in the direction of acceleration, the early time perturbation amplitude and curvature in the case of constant acceleration increases exponentially with time $\displaystyle t$ at a rate

[TABLE]

[TABLE]

where $\displaystyle\tau=\sqrt{Agk}^{-1}$ is the characteristic time scale. Introducing a finite domain decelerates the exponential growth of the perturbation amplitude for early time by an exponent of $\displaystyle\sqrt{\tanh(kZ)}$. Recall that $\displaystyle kZ$ and $\displaystyle-kZ$ are the positions of the imposed boundaries in the direction of acceleration (given in a dimensionless scaled form). The smaller the domain the greater this exponential growth is hindered. In the limit $\displaystyle kZ\rightarrow\infty$ this result agrees with the growth in an infinitely extended spatial domain [8].

3.2 Nonlinear regime – family of solutions

Identifying the evolution of RTI in the nonlinear regime is a more involved task than that in the linear regime. In RTI with constant acceleration, as $\displaystyle t/\tau\rightarrow\infty$, the velocity, moments and surface variables asymptotically approach steady values to the leading-order in time, with the next order terms decaying exponentially with time [1]. To reflect the steady-state conditions as $\displaystyle t/\tau\rightarrow\infty$ the time derivatives for both the moments and the surface curvature variable, $\displaystyle\zeta_{1}$, in the dynamical system (2.3.4) are set to zero.

The next consideration in solving the dynamical system (2.3.4) in the nonlinear regime is correctly satisfying the closure requirement. Truncating the moment expressions such that $\displaystyle\Phi_{m},\tilde{\Phi}_{m}=0$ $\displaystyle\forall$ $\displaystyle m>2$ and for an order of approximation in the spatial expansion $\displaystyle N=1$, retains the variables $\displaystyle\Phi_{1},\tilde{\Phi_{1}},\Phi_{2},\tilde{\Phi}_{2}\text{ and }\zeta_{1}$. In this case the number of variables $\displaystyle N_{t}=5$ and sufficiently overcomes the closure problem. Therefore, to obtain a reliable description of the large scale coherent dynamics of the bubble front that correctly captures the influences of higher harmonics the dynamical system in (2.3.4) is truncated to two harmonic modes.

Solving this dynamic system finds a family of all possible bubble velocity solutions parametrised by curvature. The solution family is complete and finds all the solutions supported by the prescribed global symmetries.

The analytic expression is

[TABLE]

Figure 2 illustrates the family of bubble velocity solutions for various sized finite domains of height $\displaystyle kZ$.

The solution family is plotted in dimensionless, positive units. Recall that the bubble curvature, $\displaystyle\zeta_{1}$, is a negative value, as the bubble is concave, and constant in time in the nonlinear regime for RTI. The units of velocity come from the fundamental velocity scale $\displaystyle\sqrt{g/k}$.

Two fundamental consequences of imposing a finite domain on the dynamic parameters of RTI are clearly illustrated. In a finite domain the maximum velocity solution is decreased and the greatest amount of curvature supported by the symmetries of the system is increased. These two responses are enhanced the smaller the finite domain imposed. This means that compared to RTI’s evolving in an infinite unbounded spatial domain, RTI bubble fronts growing in a finite domain travel slower with the same bubble front curvature. This result that the boundary hinders motion is physically intuitive.

Travelling at the same fixed velocity, RTI bubble structures are more curved in a finite domain than in an infinite domain. Indeed the maximum bubble front curvature supported by the symmetries of the system, denoted as $\displaystyle\zeta_{cr}$, is a monotone decreasing function of the boundary size:

[TABLE]

Note that for $\displaystyle kZ\rightarrow\infty$, $\displaystyle\zeta_{cr}=-k/2$ and for $\displaystyle kZ\rightarrow 1$, $\displaystyle\zeta_{cr}\approx-0.657k$. Figure 2b shows the effect of a finite boundary on a range of Atwood numbers. We find that the family parameter of solutions is altered uniformly for all Atwood numbers. The values for $\displaystyle kZ=1/2$ are given only for the purposes of completeness.

3.3 Atwood bubble

The parameter family of solutions can be represented by the Atwood bubble. The Atwood bubble is the fastest solution in the family with velocity $\displaystyle v_{A}$ and corresponding curvature $\displaystyle\zeta_{A}$ such that

[TABLE]

Notice that $\displaystyle\zeta_{A}$ is the curvature corresponding to the maximum velocity solution, not the maximum curvature solution which is instead denoted $\displaystyle\zeta_{cr}$. The curvature of the Atwood bubble $\displaystyle\zeta_{A}$ satisfies

[TABLE]

and $\displaystyle C>0$. The system is described using hyperbolic tangents, for convenience $\displaystyle T1$ and $\displaystyle T2$ are henceforth used as shorthand notation for $\displaystyle\tanh(kZ)$ and $\displaystyle\tanh(2kZ)$, respectively.

Transformation scalings are introduced to bridge between solutions in a finite and infinitely extended spatial domain. The fundamental scales of a fluid system in a finite domain may be expressed as a function of the infinitely extended domain.

[TABLE]

[TABLE]

where the superscript ∗ corresponds to the parameter in the infinitely extended spatial domain. These fundamental scales lead to the following additional transformation terms:

[TABLE]

[TABLE]

[TABLE]

The solution governing the local dynamics of the bubble front has a non-trivial dependence on the Atwood number and the size of the domain. The solution retains the first four harmonics but for all values the lowest order harmonics are dominant. The general expressions for $\displaystyle v_{A}$ and $\displaystyle\zeta_{A}$ are cumbersome and hence not presented here. The simplest information to be drawn from these cumbersome expressions is the existence of an invariant property satisfying

[TABLE]

and holding for all Atwood numbers. This is similar to what is found in the infinite domain case [8]. Dependency of the Atwood bubble on $\displaystyle A$ for various boundary sizes is illustrated in Figure 3.

For any fixed Atwood number RTI bubble structures travel faster and with a less curved bubble front as the domain size increases. Figure 4 demonstrates the influence of a finite domain on the diagnostic parameters as a function of the boundary size $\displaystyle kZ$. For any fixed boundary size, the RTI grows faster and with a more curved bubble front for contrasting fluid densities. Convergence is seen for both Figures 3 and 4 and as the domain size is increased the diagnostic parameters coincide with the infinitely extended case. For similar fluid densities, $\displaystyle A\approx 0$, it is well demonstrated that the bubble velocity and curvature are forced to zero tainting the accuracy of the solutions in that vicinity.

Analytic expressions for the Atwood bubble may be simplified in the limiting cases of $\displaystyle A\approx 1$ and $\displaystyle A\approx 0$ truncated to give the first two harmonics. The case of $\displaystyle A\approx 1$ is a natural consideration and the theory is applicable with no limitations. Alternatively the properties of the system for $\displaystyle A\approx 0$ is presented for completeness but it is noted that the solutions may be out of range of applicability of the theory.

For highly contrasting fluids, $\displaystyle A\approx 1$, we express the diagnostic parameters to the first-order as

[TABLE]

[TABLE]

For similar density fluids, $\displaystyle A\approx 0$, the fastest solution variables are explicitly

[TABLE]

3.4 Taylor’s bubble

Taylor’s bubble is named here as the particular solution for the Atwood bubble in a one fluid system, $\displaystyle A=1$. In this case the problem variables have the form

[TABLE]

Notice that due to the $\displaystyle(-T1+T2)$ term, the second harmonic of the heavier fluid goes to zero $\displaystyle\Phi_{2}\rightarrow 0$ as the boundary height, and hence domain size, increases $\displaystyle Z\rightarrow\infty$.

3.5 Qualitative velocity fields

In the approximation of a potential flow the velocity of the bubble is given by

[TABLE]

Taking the partial derivatives of the potential function, $\displaystyle\frac{\partial\Phi}{\partial x}$ and $\displaystyle\frac{\partial\Phi}{\partial z}$, we may examine the tangential ($\displaystyle x-$direction) and normal ($\displaystyle z-$direction) components of the velocity as a function of $\displaystyle kz$, the dimensionless displacement from the interface. Here we consider only small distances from the interface, $\displaystyle z\approx z^{*}$. The system is described such that the heavy (light) fluid is situated at $\displaystyle kz>0$ ($\displaystyle kz<0$). This enables a specific examination of the influence of a finite boundary on the vector components of the fluid motion, see Figure 5.

The results are given in the inertial, laboratory frame of reference and the physically significant (fastest) solution in the family at $\displaystyle A=1$ is examined. We find that the normal component of velocity is greatest and continuous at the interface. Indeed, this continuity condition was imposed as one of the interface conditions, that is, the condition of no mass flux across the interface. The tangential component of velocity is discontinuous, implying that a shear function develops at the interface. This presence of a shear function leads to the Kelvin-Helmholtz small scale vortical structures. The value of the tangential component of velocity in Figure 5a is low at a short horizontal distance from the bubble tip, $\displaystyle kx=10^{-3}$. For both tangential and normal components, anisotropy is observed in the direction of acceleration as the velocity profile is asymmetric about the origin. We observe that this anisotropy is enhanced in a smaller finite domain. Further enforcing a smaller spatial domain increases the magnitude of the tangential component of velocity while the normal component decreases. This indicates that the shear function is greater (and hence we can expect larger vortical structures) in RTI evolving in a smaller, finite domain when compared to a larger or infinite spatial domain.

Figure 6 shows the qualitative velocity fields of RT flows both in the laboratory and non-inertial bubble tip frame of reference. The dashed line represents the interface of the fluids. As described in the system description, the heavy fluid layer is located above the interface and the light fluid below. The constant externally imposed acceleration is downwards, from heavy to light.

In the laboratory frame of reference, Figure 6a, RTI is characterised by intense motion of the fluids in the vicinity of the interface and effectively no motion away from the interface. A shear function is present at the interface and is greater further from the bubble tip which may cause vortical structures. There is no shear function exactly at the bubble tip, $\displaystyle x=0$. The method of solution applied to the large-scale dynamics of the instabilities finds the large scale ‘envelope’ properties of the small scale vortical structures [8]. According to the velocity field in Figure 6a the vortical structures rotate ‘inwards from light to heavy’.

In a finite domain the fastest solution travels with a greater curvature and the inter-facial vortical structures appear more pronounced when compared to the infinite domain evolution [8], in agreement with the results of Figure 5a. At the upper and lower boundaries the normal component of velocity vanishes and only a tangential component remains as set by the boundary conditions

This velocity field agrees qualitatively with accurate experiments and simulations and suggests that linear and nonlinear RT dynamics in incompressible immiscible fluids (1) is essentially inter-facial, with intense motion of the fluids near the interface and effectively no motion of the fluid away from the interface; (2) has potential flow in the bulk and vortical structures at the interface. Our theoretical results are valid in the vicinity of the bubble tip. To fully describe the mushroom-type shape of the spike, further investigations are required, including the non-ideal effects to be done in the future.

4 Shear function analysis

Previous studies of RTI consider a bubble curvature parametrisation to the family of solutions. This means that the diagnostic parameters are given as a function of curvature. A parametrisation of curvature is favoured as it is observable and hence comparisons with experiments are simplified. A closer analysis, however, indicates that the shear function may be a more natural and suitable parameter to the solutions. The shear function is defined here as a quantification of the discontinuity of the tangential component of velocity. Indeed, it is this discontinuity that is responsible for the multiplicity of solutions and is believed to significantly influence the nonlinear dynamics [8]. Mathematically the shear function, $\displaystyle\Gamma$, is defined as

[TABLE]

where $\displaystyle v_{h}(x,0)$ and $\displaystyle v_{l}(x,0)$ are the tangential components of the velocity field in a vicinity of the interface between the heavy and light fluid, respectively. The shear function is an important dynamical parameter as it drives the Kelvin-Helmholtz vortical structures. The theoretical approach used in this work allows the quantification of the shear function at the interface in the vicinity of the bubble tip. Although the methodology limits our analysis to the large-scale dynamics, considering the shear function enables us to determine the large scale qualitative properties of the vortical structures such as the relative strength and rotation direction. In the linear regime there is no shear function in the vicinity of perturbation front, this result fits with observations as the vortical structures are not visible for early time.

In the nonlinear regime the shear function is clearly present due to the noticeable Kelvin-Helmholtz vortical structures. To leading-order in terms of the moment equations for RTI the shear function is given, in terms of moments, as

[TABLE]

4.1 Shear function as a diagnostic parameter

Solving (31) to the first-order and taking the first two harmonic modes finds that the shear function is dependent on boundary size, bubble curvature and Atwood number. The explicit analytic expression is cumbersome and so not presented here, instead Figure 7 demonstrates the dependencies of the shear function. The shear function is a non-monotone function of curvature and has a maximum corresponding with the critical curvature defined previously $\displaystyle|\zeta_{cr}/k|=|\coth(kZ)/2|$. Truncated to the region $\displaystyle 0\leq|\zeta_{1}|\leq|\zeta_{cr}|$ the shear function is in fact directly related and a one-to-one function of curvature. This range consists of all the possible shear function values supported by the system.

For any fixed curvature value a decreasing domain size corresponds to an increasing shear function. RTI bubble structures evolving in a finite domain have a greater amount of shear function (and hence more distinct vortical structures) at the interface when compared to their evolution in a larger or infinite spatial domain. Irregardless of the size of the domain, instabilities growing in a highly contrasting fluid system $\displaystyle A\approx 1$ have a larger contribution from inter-facial shear function than those with similar fluid densities $\displaystyle A\approx 0$. Further, instabilities at all Atwood numbers converge at the same rate with increasing boundary size.

Particular attention should be focused on the solution with maximum shear function. This solution is denoted the critical bubble with the maximum shear function value given as $\displaystyle\Gamma_{S}$ and corresponding curvature as $\displaystyle\zeta_{S}$. The critical bubble satisfies the equation

[TABLE]

Solving this finds the analytic expression for the maximum shear function value and corresponding bubble front curvature respectively as

[TABLE]

The dimension of the shear function, $\displaystyle s^{-1}$, comes from the foundational scale $\displaystyle(gk)^{-1/2}$. The maximum shear function value increases for contrasting fluid densities and for smaller domain sizes. The bubble curvature that maximises shear function is independent of Atwood number and is a function of the boundary size. RTI bubbles evolving in a smaller domain experience maximum inter-facial shear function for more curved bubble fronts. In a finite domain $\displaystyle|\zeta_{S}|<|\zeta_{cr}|$, whereas in the limit of an infinite domain with $\displaystyle kZ\rightarrow\infty$ the values are

[TABLE]

4.2 Shear function as a parameter to the family of solutions

The direct link between the multiplicity of solutions and the shear function is explored by considering the shear function as an alternate parameter to the family of solutions. Figure 8 illustrates the family of bubble velocity solutions in the nonlinear regime as parametrised by the shear function.

For each Atwood number the range of velocity values is conserved whether using a curvature or a shear function parametrisation. An interesting property of Figure 8 is that the bubble velocity appears independent of Atwood number until the vicinity of maximum shear function (which is dependent on Atwood number) in which the solution falls quickly to zero. RTI evolution in a smaller finite domain is characterised by a lower maximum velocity and a greater range of possible inter-facial shear function values.

Figure 9 illustrates the family of bubble curvature solutions as parametrised by the shear function. A direct, monotone relationship of the two diagnostic parameters is observed. Using the shear function as an alternate parameter has the benefit of tying together the diagnostic parameters in the bubble front vicinity.

5 Study of convergence

The regular asymptotic solutions for the late-time evolution involve multiple harmonics and for all Atwood numbers the lowest order harmonic is dominant [16]. This allows us to study the convergence properties of the solution. To prove a good convergence of the solution, it should be shown that for all harmonic modes $\displaystyle m$, $\displaystyle\Phi_{m+1}>\Phi_{m}$. In this analysis, in which the moments are truncated to retain two harmonic modes, the convergence of the solution is investigated by comparing the magnitude of the first, $\displaystyle\Phi_{1}$, and second, $\displaystyle\Phi_{2}$, Fourier amplitude modes. In the case that $\displaystyle\Phi_{1}>\Phi_{2}$ indicates that the lowest harmonic is dominant and our solution has good convergence. Both Fourier amplitudes are plotted for $\displaystyle 0\leq|\zeta_{1}|\leq|\zeta_{cr}|$ and the corresponding shear function region.

The regions in which $\displaystyle\Phi_{1}$ represented by a continuous line lies below $\displaystyle\Phi_{2}$, the dashed line, indicates regions of poor convergence. In Figure 10 we explore the convergence properties of the heavy fluid for various Atwood numbers at a fixed boundary size $\displaystyle kZ=5$.

6 Conclusion

One of the fascinating things about fluid flow is the incredible richness of the mathematics and physics that appears in the problems. The challenge of describing RTI flows is no exception. We learn the importance of accurate approximations to deal with these flows. In this work, approximations such as neglecting non-idealised fluid effects, a scale separation and approximating a potential flow enables the Navier-Stokes equations to be solved and yields a description of RTI in the vicinity of the bubble tip.

Long-standing problems such as this, challenge us to find innovative techniques and approaches. The group theory approach with a strong consideration for the symmetries of the flow, is one such creative technique that helps find a description of RT evolution. The large scale bubble structure have symmetry properties defined by the $\displaystyle pm1$ symmetry group.

Under the influence of a constant acceleration the initial wave perturbation curvature and growth rate increase exponentially with time. Imposing a finite domain hinders this growth by an exponent of $\displaystyle\tanh(kZ)$ dependent on dimensionless boundary position $\displaystyle kZ$.

In the more involved nonlinear regime, solving the Navier-Stokes equations for the large-scale bubble structures with symmetry group $\displaystyle pm1$ finds a parameter family of solutions supported by the symmetries. The fastest solution in this family is physically significant. For RTI the fastest solution corresponds to a curved bubble. Introducing a finite domain decreases the maximum velocity and increases the possible bubble front curvatures. For completeness a description was found for all Atwood numbers, however as the scale separation breaks down for similar density fluids we should be cautious in this vicinity.

The multiplicity of solutions in the nonlinear regime is a consequence of the governing boundary condition enabling inter-facial shear to develop. In RTI flows the shear function is a monotone increasing function of curvature. Using the shear function as an alternate parameter to the multiplicity of solutions has the benefit of tying together the observables and improving the convergence properties of the harmonic amplitudes.

7 Further work

The work discussed in this paper sets the groundwork for further research. Specifically, we should expand our understanding of the evolution of RTI in a finite domain by considering the three-dimensional case. A group theory approach is applicable in three dimensions and so the method of solution discussed here is appropriate with some alterations. The properties of the diagnostic parameters in the transition between three and two dimensional spatially bounded flows should also be considered. The evolution of the bubble front in a finite domain was described only to a first order approximation, $\displaystyle N=1$. To better the accuracy of our description of RTI evolution we should consider higher orders of approximation.

Our suggestion that the shear function provides a better parametrisation to the family of solutions should also be explored further. The stability of the diagnostic parameters as a function of the shear function should be found by perturbing the solutions or using bifurcation theory. It must also be checked that the dependence of the diagnostic parameters on time in the nonlinear regime remains the same in a shear function parametrisation. For example as a function of curvature, velocity is constant in time in the nonlinear regime, $\displaystyle v\sim O(1)$. However with a different parametrisation this should be tested accordingly.

A closer analysis of inter-facial shear may reveal properties of the Kelvin-Helmholtz vortical structures. Specifically we may identify the influence of non-idealised properties of the fluids, such as viscous stress and surface tension, on the shear function. Numerical simulations should be run including non-idealised effects. The results provided here give a quantitative ‘ideal’ comparison to any future numerical simulations. This may also reveal if the non-idealised effects hinder the shear function near the bubble-tip and hence explain the discrepancy of the vortical structures in the vicinity of the bubble tip.

This paper is ended with a discussion on the possible avenues to expand this work as a motivation for future studies. This motivation stems from the need to understand and control RTI in a range of industrial applications. Although our understanding of RTI flows has greatly enhanced since Lord Rayleigh first defined the fluid problem more than $\displaystyle 100$ years ago, there is still a need for the development of powerful theoretic approaches to facilitate a greater understanding of these fascinating flows. However, describing the influence of a finite domain and exploring the idea that using the shear function as a parameter may lead to more accurate results, brings us one step further.