The Rayleigh-Taylor instability (RTI) in an infinite slab where a constant density lower fluid is initially separated from an upper stratified fluid is discussed in linear regime. The upper fluid is of increasing exponential density and surface tension is considered between both of them. It was found useful to study stability by using the initial value problem approach (IVP), so that we ensure the inclusion of certain continuum modes, otherwise neglected. This methodology includes the branch cut in the complex plane, consequently, in addition to discrete modes (surface RTI modes), a set of continuum modes (internal RTI modes) also appears. As a result, the usual information given by the normal mode method is now complete. Furthermore, a new role is found for surface tension: to transform surface RTI modes (discrete spectrum) into internal RTI modes belonging to a continuous spectrum at a critical wavenumber. As a consequence, the cut-off wavenumber disappears: i.e. the growth rate of the RTI surface mode does not decay to zero at the cut-off wavenumber, as previous researchers used to believe. Finally, we found that, due to the continuum, the asymptotic behavior of the perturbation with respect to time is slower than the exponential when only the continuous spectrum exists.

The Rayleigh-Taylor instability, which occurs due to the gravitational instability of a heavy fluid overlying a lighter fluid1,2 is of fundamental importance in a multitude of applications ranging many fields, such as atmospheric (mushroom clouds in volcanic eruptions and atmospheric nuclear explosions..), astrophysics (supernova explosions), chemistry (reaction-diffusion systems and immiscible fluids in porous media)...

The interface between two different liquids that do not mix is characterized by a surface energy density that is related to the asymmetry of the intermolecular forces acting on the molecules in the interface. When the interfacial area is increased, external energy must be supplied to bring additional molecules into the interface. An equivalent measure of the interfacial energy density is by means of a surface tension, which acts in a direction tangent to the surface so as to minimize its area.3 

Chandrasekhar4 generalized Taylor’s problem by including surface tension. In the absence of surface tension and viscosity, the growth rate increases indefinitely with the wavenumber. According to Chandrasekhar, the presence of surface tension tends to inhibit the growth of the instability. Moreover, for the incompressible case there is a cut-off wavenumber kc, so that the arrangement is stable for k>kc, i.e., the presence of surface tension stabilizes perturbations for the wavenumber larger than a cut-off value. Reid5 has also discussed the effects of surface tension and viscosity on the stability of two superposed fluids. He has discussed the most detailed and considerable effects arising from surface tension. In addition to this, Daly6 made a numerical study of the effect of surface tension on interface instability. He showed how surface tension affects the growth of the Rayleigh-Taylor spike and provides the mechanism for drop separation from the spike.

Vaghela and Chhajlani7 investigated the Rayleigh-Taylor instability of the plane interface separating the two partially-ionized superposed fluids through porous medium. Here, the effect of variable horizontal magnetic field, surface tension and rotation along the vertical axis are also incorporated. They concluded that the system is unstable only for those positive wavenumbers which are less than certain critical value in case of an adverse density gradient.

Mikaelian8 has studied the RTI and Richtmyer-Meshkov instabilities in multilayer fluids with surface tension. He analyzed the general case of N fluids with arbitrary densities and surface tensions and derived the eigenvalue equation determining the growth rate of the perturbations. Besides, he studied turbulence at the interface between two fluids with surface tension and presented specific predictions for the turbulent energy as function of the surface tension. Al Ansary et al.9 discussed theoretically the effects of surface tension and uniform rotation of a vertical axis for three-fluid systems. The fluids are considered to be incompressible with varying density. Numerical results were obtained for two-fluid and three-fluid systems. He noticed that, when both separating surfaces have non-vanishing surface tensions, the short wavelength perturbations will become completely stabilized.

Recently, Chertkov et al.10 have investigated the effect of surface tension on immiscible Rayleigh-Taylor turbulence. They presented phenomenology describing the internal structure of a turbulent zone, produced as the result of the push of a heavy fluid into a light one, for the case of immiscible fluids. Sharma et al.11 have analyzed the RTI in case of fully ionized superposed fluid incorporating the effect of surface tension.

Obied Allah12 has discussed the problem of Rayleigh-Taylor instability in the presence of surface tension, stratifications and two rigid planes. Besides, he analyzed some particular cases: only one rigid plane and rigid planes at infinite. He found that the stratification has a destabilizing effect. Besides, he showed that the presence of one rigid plane or two rigid planes has a stabilizing effect.

From all these discussions, we also conclude that various problems of RTI in the presence of surface tension are under current investigation, but none of the authors has studied how the surface tension affects the RTI surface mode in the presence of a continuous spectrum seeded by a density gradient.

Several authors have purposed solutions of the initial value problem through the use of a Laplace Transform in time. The principal finding of the initial value problem (IVP) approach is that, in addition to the discrete eigenvalues linked to the normal modes, there exists a continuous spectrum of eigenvalues, so, the modal approach cannot provide a complete solution. Thus, Case13 using IVP showed that for appropriately posed problems, the solutions of the linearized Navier-Stokes equations approach those of the linearized Euler equations as the viscosity tends to zero. Moreover, Case14 investigated the stability of an idealized inverted atmosphere. He showed that an initial perturbation grows in a similar way to that found by Taylor for the instability of accelerated fluid surfaces.

On the other hand, Case15 analyzed the stability of an idealized atmosphere with constant shear flow and exponentially decreasing density, finding that no initial perturbation increases in time. Consequently, the decrease is, in general, not faster than a small fractional power of the time. Burger16 solved a general initial value problem for wave perturbations of small amplitude on a zonal atmospheric flow with constant vertical shear and vanishing temperature lapse rate, using frictionless, hydrostatic, adiabatic and geostrophic theory and beta-plane dynamics. His analysis showed that, for the exceptional wavelengths for which no such unstable normal mode exists, the continuous spectrum contribution is unstable.

Sedlácêk17 investigated small amplitude electrostatic oscillations in a cold plasma with a continuously varying density. He constructed the Green’s function of the differential equation of the problem to show that there were branch-point singularities on the real axis of the complex frequency-plane. In addition, he found an infinity of new singularities (simple poles) of the analytic continuation of the Green’s function into the lower half of the complex frequency-plane.

Ott and Russell18 investigated smooth decreasing density profiles in the presence of a gravitational field. This survey showed that the Rayleigh-Taylor instability exhibits essentially different behavior above and below a certain critical wavenumber kc. Menikoff et al.19 studied the IVP problem associated with the development of small amplitude disturbances in Rayleigh-Taylor unstable, viscous, incompressible fluids. In addition to normal modes, a set of continuum modes, not treated explicitly in the literature, makes an important contribution to the development of the fluid motion.

Gustavsson20 analyzed the behavior of small disturbances in a boundary layer flow by IVP. It is found that a disturbance evolves not only as discrete waves but also has a portion described by a continuous spectrum. Similarly, Ott and Russell21 presented a linear theory of the Rayleigh-Taylor instability in the equatorial ionosphere. For a purely exponential density profile, they found that no unstable eigenmode solutions exist. To clarify the situation in the case that eigenmodes do not exist, they solved the initial value problem for the linearized ion equation of motion in the long time asymptotic limit. Finally, Bogdan and Cally22 considered the linear oscillations of a plane-parallel semi-infinite electrically conducting atmosphere with a constant temperature gradient, subjected to an imposed uniform gravitational acceleration and uniform magnetic field. The results show that atmospheric motions with prescribed horizontal variations of the form exp(ikx), with k real, possess both a discrete set of complex eigenfrequencies ωn, n=0,1,2, and a continuous spectrum.

This paper is organised as follows. In Sec. II the stratified fluid model, the governing equations and linear stability analysis are presented. A discussion via the initial value problem (IVP) is then given in Secs. II and III. The connection of this approach with a correct normal mode treatment is indicated in Sec. IV. In Sec. V the temporal evolution of disturbance is analyzed showing that the application of the hydrodynamic-stability theory care must be taken not to omit continuous spectra.

We consider a similar framework to the one used in Ref. 12, where a vertically stratified fluid at the top and an incompressible fluid at the bottom rest in a stationary Cartesian coordinate system x,y,z in the presence of a gravity field g = (0,0,−g). Typical examples of stratified fluids are the Earth’s ocean or its atmosphere. In the unperturbed state we consider both fluids at rest, i.e., uo = 0. In this state the interface separating the two fluids at zo = 0 is an horizontal plane. The suffix zero will be used to characterize the equilibrium quantities. We suppose that the fluids are semi-infinite. The density distribution ρo(z) is given by equation (1) and reported in Fig. 1 

(1)
FIG. 1.

Density profile at the unperturbed state.

FIG. 1.

Density profile at the unperturbed state.

Close modal

where ρ1 and ρ2 are independent constants, being ρ2>ρ1, and β represents the inverse scale height in its corresponding region.

In the static state of a fluid arranged in horizontal layers in a uniform gravitational field g, the density ρo and the pressure po are functions of the coordinate z satisfying the equation of hydrostatic equilibrium

(2)

This means that the pressure at any level is determined by the total weight of the fluid above this level.

To investigate the character of the equilibrium of this system, we shall consider the evolution of a small disturbance. So, we now introduce a small disturbance in the system and linearized around the unperturbed state. Thus, the first order linearized equations of conservation of mass, motion and incompressibility are, respectively

(3)
(4)
(5)

where u(u,v,w), δp and δρ denote the perturbations in the velocity, pressure and density, respectively. The location of the perturbed interface is described by zs + ξs(x, y, t), where the function ξs(x, y, t) with ξs/t=ws is the perturbation of the interfacial surface and the subscript s indicates the interfacial surface.

Introducing Fourier transforms23 with respect to x and y,

(6)
(7)
(8)
(9)
(10)
(11)

and taking account the influence of the surface tension,4 we obtain

(12)
(13)
(14)
(15)
(16)

Here, δ(zzs) denotes the Dirac function and Ts is the surface tension of interfacial surface. Moreover, k2=kx2+ky2 and ω are both the square wavenumber and the growth rate of the harmonic perturbations.

After eliminating δpk, δρk, uk and vk from the above equations, an equation for wk is obtained as

(17)

Next, we consider the IVP and introduce the Laplace transform.23 So we define the forward Laplace transform in the complex ω -plane as

(18)

Here we will use the second derivative property of the Laplace Transform

(19)

being F1 and F2 the initial conditions that are imposed at t=0

(20)

So, application of this transformation to equation (17) gives

(21)

where H(z; ω, k) term is given by

(22)

and the Γz;ω,k term involves all initial conditions being

(23)

This differential equation (21) is now to be solved for w¯k, which represents the general relation of RTI of two superposed fluids incorporating the effect of surface tension. For simplicity, and according to Ott and Russell,21 we choose an initial condition corresponding to the Green’s function of equation (21) namely

(24)

where δ(z + z1) is the Dirac delta function centered at z = −z1 and Γo is a constant.

To proceed further, the solutions of equation (21) must satisfy certain boundary conditions. Since we are interested in modes whose behavior is dominated by the discontinuity in the density at the interface separating the fluids, a suitable boundary condition is that w¯k0 as z. Additional conditions must be satisfied on the horizontal plane (z = zs) on which there are density discontinuities. Evidently, one condition is that w¯k is continuous on the interface. Following to Chandrasekhar4 and Sharma et al.,11 a second condition is obtained by integrating equation (21) across the boundary z = zs between zsε and zs + ε and then letting ε0. In virtue of the continuity of w¯k across z = zs and the boundedness of ρo this limiting process leads to (25)

(25)

where Δs(f) = f(zs + 0) − f(zs − 0) is the jump, that a function f experiences at the interface z = zs. The subscript s distinguishes the value of a continuous quantity known at the interface z = zs.

Concerning to the initial condition given by equation (24), w¯k(z,z1;ω,k) must be continuous at z = −z1. Besides, the first derivative of w¯k(z,z1;ω,k) is discontinuous by a unit amount at z = −z1, so that the second derivative then yields a delta function singularity. Thus, we obtain the following jump condition at z = −z1

(26)

Taking account the above conditions, the general solution of (21) is a linear combination of two linearly independent solutions. Thus, w¯k(z,z1;ω,k) becomes

(27)

where α is given by equation being

(28)

where Re(α)>0 and the constants A, B, C and D are independent from z. Those constant measure the amplitude of the perturbation and they are determined by the boundary and matching conditions as

(29)
(30)
(31)
(32)

For simplicity, we use the solution for z>0

(33)

where the dispersion function Δ(k, ω) is given by

(34)

The zeros of dispersion function Δ(k, ω) are calculated by the dispersion relation Δ(k, ω) = 0, and represent the discrete spectrum. This expression, can be solved to produce a relation of the type ω = ω(k). It is not difficult to reduce equation (34) to a polynomial expression

(35)
(36)
(37)
(38)

Consequently, w¯k(z,z1;ω,k) has two poles ±ωp , i.e. the zeros of the dispersion function Δ(k, ω). The observation of this result is congruous with the Obied Allah’s findings.12 He obtained the same dispersion relation by using a normal-mode analysis. Before proceeding with a detailed discussion of the results of these computations, we shall consider the asymptotic behavior of the solutions of equation (35) for small β. Thus, as the fluid of exponentially varying density becomes constant, i.e., in the limit of a β0, the dispersion relation is

(39)

Consequently, we can appreciate that the relation derived by Chandrasekhar is recovered.4 

The continuous spectrum arises as a result of the multi-valuedness of w¯k. Thus, expression (33) is the basis of our investigation. Following to Ott and Russell,18 due to the fact that the function w¯k does not turn out to be an even function of α, this function has branch points at ω=±ωb (α = 0) on the real axis associated with a continuous spectrum, where

(40)

If the ω -plane is cut along the interval connecting - ωb and ωb, w¯k is a single-valued function of ω. This cut [−ωb, +ωb] is the interval of the continuous spectrum.

The complex inversion formula is a very powerful technique for computing the inverse of a Laplace transform. The technique is based on the methods of contour integration and requires that we consider our parameter ω to be a complex variable.

To proceed further, we must invert w¯k(z,ω) to obtain wk(z, t) as function of time. The inverse Laplace transform of w¯k(z,ω) is given by

(41)

The expression (41) is known as the complex (or Fourier-Mellin) inversion formula. The Bromwich contour σ is an appropriate integration path (see figure 2.a) which excludes the singularities of w¯kz,ω in the region Reωσ. This contour runs parallel to the imaginary axis and the constant σ has to been chosen such the path of integration in the complex ω -plane is on the right of all singularities of w¯k.The evaluation of the integral (41) is done with the method of residues, and the integration path γ is, therefore, closed to the left ω-half-plane so that all the singularities of w¯kz,ω are enclosed (see figure 2.b) by a semicircle denoted by Γ. Now, we deform the contour of integration without affecting the value of the integral as long as we do not cross any singularities of the integrand during the deformation. Therefore, w¯k(z,ω) is then a single-valued analytic function in the region inside the lines Γ + γ, C1, C2 and Cb. In accordance with the boundary condition w¯k±,z1;ω,k=0, the square root given by (28) is defined, so that Re(α)>0 for Re(ω)>σ. Thus, this integral (41) in the ω- plane can now be written as

(42)
FIG. 2.

a) The Bromwich integration path. The dots on the real axis refer to branch points, and the wiggly lines refer to branch cuts, representing the appropriate Riemann sheet. Simple poles are indicated by crosses. b) Example of lines used for the integral evaluation (41).

FIG. 2.

a) The Bromwich integration path. The dots on the real axis refer to branch points, and the wiggly lines refer to branch cuts, representing the appropriate Riemann sheet. Simple poles are indicated by crosses. b) Example of lines used for the integral evaluation (41).

Close modal

where Res means the residues values at the poles ωn and the summation runs over all poles enclosed (two poles) which represent the discrete spectrum.

The integral

(43)

can be expressed as a linear combination of terms of the type eωt where ω is real and runs between ±ωb. These values clearly represent a continuum of eigenvalues. The branch cut along the interval connecting − ωb and ωb is associated with this unstable continuous spectrum. Consequently, we can consider that (43) is the integral on the continuous spectrum. Since boundaries Γ and γ both extend to infinity and the integrals on C1 and C2 cancel each other, the first four terms in equation (42) disappear, then

(44)

Let us examine the impact of the continuous spectrum on the discrete spectrum.

It is always useful to convert the expression (33) of wk for z>0 to a dimensionless form. So, the most relevant dimensionless parameters become

(45)

On the other hand, we define the Atwood number as

(46)

where r=ρ2ρ1 is the density jump. Hence, w¯k for z̃>0 may be written in the normalized form

(47)

where the denominator of (47) is the normalized dispersion relation Δ(k̃,ω̃)

(48)

Fig. 3 shows the normalized square growth rate ω  ̃2 against the normalized wavenumber k̃ obtained when Δ(k̃,ω̃)=0 in equation (48) for At = 0.4 and At = 0.7 and τ = 0.01 and τ = 0.02. The boundary of the region of the (k̃,ω̃2) -plane from which surface modes are excluded, see equation (40), is shown by a dashed line (normalized square growth rate of the branch points named as ω̃b2) in each plot. The solid lines stand for the normalized square growth rate of the poles (discrete spectrum) ω̃p2. As it can be appreciated, in every branch, there are two terminal points with critical wavenumbers k̃cr and k̃cr+, respectively, that limit the dispersion relation domain. As can be observed, for k̃>k̃cr+ or k̃<k̃cr the ω̃p2 (poles) move towards ω̃b2 (branch points). Due to role of surface tension here, surface RTI modes (discrete spectrum) have a maximum and transform into internal RTI modes at k̃cr+. This is possible when the continuous spectrum, neglected by Obied,12 is now taken into account. Thus, for k̃=k̃cr±, then ω̃p2=ω̃b2 and the poles move to another sheet of α (the surface modes disappear). In the absence of continuous spectrum, the discrete spectrum domain would be described by the equation (39). However, the intersection of the discrete and continuous spectra at k̃cr+ is a consequence of the combined effect of the surface tension and the existence of internal modes. We can also observe in Fig. 3, that for a fixed Atwood number At, as τ increases the domain [k̃cr,k̃cr+] where the discrete modes appear will be reduced. On the other hand, when the value of At is increased for a fixed dimensionless surface tension parameter τ, the discrete mode domain will be larger. Moreover, here it is shown how the critical wave numbers k̃cr and k̃cr+ depend on the density jump, indicating the sensitivity of surface mode to At changes.

FIG. 3.

The dependence of the square of the dimensionless growth rate on dimensionless wavenumber. The results are obtained for τ = 0.01 and τ = 0.02 and different Atwood numbers At = 0.4 and At = 0.7. Solid lines correspond to our numerical results (discrete modes). Dashed lines correspond to branch point. For the sake of clarity, the k̃cr and k̃cr+ values are only indicated for τ = 0.02.

FIG. 3.

The dependence of the square of the dimensionless growth rate on dimensionless wavenumber. The results are obtained for τ = 0.01 and τ = 0.02 and different Atwood numbers At = 0.4 and At = 0.7. Solid lines correspond to our numerical results (discrete modes). Dashed lines correspond to branch point. For the sake of clarity, the k̃cr and k̃cr+ values are only indicated for τ = 0.02.

Close modal

At the secondary level, we estimate numerical values for k̃cr, which are easily derived by combining equations (35) and (40). As k̃cr depends on At and τ, we studied this dependence with At for a fixed τ. Thus, in Fig. 4 we have represented the critical wavenumber k̃cr in terms of the dimensionless parameter At for τ = 0.01 and τ = 0.02. The plot presents two branches: the lower and the upper branch for k̃cr and k̃cr+, respectively. Also, it can be shown in Fig 4 that when At is decreased for a fixed τ value, the range in which the surface mode exists also decreases. It is worth mentioning that for a fixed τ, as the value of At is further decreased, the two terminal points approach two threshold limit values k̃cr* and k̃cr+* on each branch, respectively, for a terminal Atwood number At*. These findings provide evidence that suggests that the necessary and sufficient condition for the existence of the discrete spectrum for a fixed τ is that At>At*. Thus, for instance, for τ = 0.01 the terminal wavenumbers are k̃cr*3.5395 and k̃cr+*3.5423, being At*0.208942 the terminal Atwood number.

FIG. 4.

Dimensionless critical wave number k̃cr as a function of the dimensionless parameter At for different τ dimensionless parameter. The lower and the upper branch correspond to k̃cr and k̃cr+ respectively.

FIG. 4.

Dimensionless critical wave number k̃cr as a function of the dimensionless parameter At for different τ dimensionless parameter. The lower and the upper branch correspond to k̃cr and k̃cr+ respectively.

Close modal

Taking into account the previous analysis, the square of the dimensionless growth rate against the dimensionless wavenumber close to At* was plotted for τ = 0.01 and two different Atwood numbers, see Fig. 5.

FIG. 5.

Dependence of the square of the dimensionless growth rate on dimensionless wavenumber for τ = 0.01 and different Atwood numbers At = 0.22 and At = 0.2089422. Solid lines correspond to our numerical results (discrete modes). Dashed lines correspond to branch point.

FIG. 5.

Dependence of the square of the dimensionless growth rate on dimensionless wavenumber for τ = 0.01 and different Atwood numbers At = 0.22 and At = 0.2089422. Solid lines correspond to our numerical results (discrete modes). Dashed lines correspond to branch point.

Close modal

It is remarkable that, although the problem analyzed here is a special case considered by Obied Allah,12 this author did find no evidence of a continuous spectrum. On the other hand, our results seem to refute the common belief that, due to the existence of surface tension, the growth instability for the RTI surface mode goes to zero for the cut-off wavenumber kc. This result is not surprising, because the cut-off wavenumber does not exist at this stage, due to the existence of the continuous spectrum. Moreover, Obied Allah12 did not make any explicit comment about the nonexistence of this cut-off wavenumber.

Let us illustrate how the normalized w̃k(z̃,z̃1;ω̃,k̃) behaves as a function of both normalized real part and normalized imaginary part growth rate, respectively. To study graphically the existence of continuous spectrum, we consider a contour plot of the imaginary part of w̃k(z̃,z̃1;ω̃,k̃) for At = 0.7, τ = 0.01, k̃=6, z̃1=0.5, z̃=0 and Γ̃o=1, see Fig. 6.

FIG. 6.

Contour plot of normalized imaginary part of w̃k(z̃,z̃1;ω̃,k̃) versus normalized growth rate, for At = 0.7, k̃=6, z̃1=0.5, z̃=0, τ = 0.01 and Γ̃o=1. It illustrates the branch cut running along real axis between Re(ω̃) =-1 and Re(ω̃) =1. Each pair of lobes represents a pole.

FIG. 6.

Contour plot of normalized imaginary part of w̃k(z̃,z̃1;ω̃,k̃) versus normalized growth rate, for At = 0.7, k̃=6, z̃1=0.5, z̃=0, τ = 0.01 and Γ̃o=1. It illustrates the branch cut running along real axis between Re(ω̃) =-1 and Re(ω̃) =1. Each pair of lobes represents a pole.

Close modal

This clearly illustrates the continuous spectrum through the branch cut running along real axis between Re(ω̃)=1 and Re(ω̃)=1. On the other hand, poles of w̃k(z̃,z̃1;ω̃,k̃) function can be appreciated due to the existence of pairs of lobes at Re(ω̃)=±1.626.

Now, we plot the contour plot of the imaginary part of w̃k(z̃,z̃1;ω̃,k̃) versus normalized growth rate for At = 0.7, τ = 0.01, k̃=20, z̃1=0.5, Γ̃o=1 at z̃=0, (see Fig. 7). In this case k̃cr+8.9, as k̃>k̃cr+, there is not a discrete spectrum and only the branch cut Re(ω̃)=1 and Re(ω̃)=1 appears.

FIG. 7.

Contour plot of normalized imaginary part of w̃k(z̃,z̃1;ω̃,k̃) function versus normalized growth rate for At = 0.7, k̃=20, z̃1=0.5, z̃=0, τ = 0.01 and Γ̃o=1.

FIG. 7.

Contour plot of normalized imaginary part of w̃k(z̃,z̃1;ω̃,k̃) function versus normalized growth rate for At = 0.7, k̃=20, z̃1=0.5, z̃=0, τ = 0.01 and Γ̃o=1.

Close modal

Our goal in this section is to analyze the time evolution of a perturbation and observe the transient response of the system.

Before proceeding further, it may be of some interest to treat the nature of the unstable growth in those cases where no solutions are obtained from the eigenvalue problem and only continuum exists. The growth of the response to an initial perturbation will be studied, firstly, for the cases k<<kcr+ and k>>kcr+ and, secondly, for k=kcr and k=kcr+. Consequently, the only contribution to w¯kz,ω comes from the branch point at ω = ωb. To calculate the asymptotic behavior of wkz,t we must find an approximate expression of this perturbation valid at large time values.

To transform the problem into the time domain, we take the inverse transform of w¯kz,ω by integrating along the contour σ. To find the asymptotic behavior of the integral (41), we deform the contour towards the left until the first singularity is encountered, which is the branch point at ωb, see Fig. 8. Note that σ is larger than the real parts of all singularities of w¯k, except the one at ω = ωb.

FIG. 8.

The contour σ used to invert w̃kz,ω.

FIG. 8.

The contour σ used to invert w̃kz,ω.

Close modal

To obtain the time-asymptotic solution wkz,t, we require t to be so large that the integral along Re(ω)=σ can be neglected when compared to the integral obtained in the proximities of ωb. Since we may take σ to be arbitrarily close to but greater than zero, the error we make in neglecting the integral along the regions of Re(ω)=σ far from ωb decrease as exp(ωbt).The asymptotic solution can be expressed as

(49)

Using the expressions (28),(40) and (33) we replace ωωb by reiθ in the integrand of (49). In the vicinity of the branch point, we replace ω by ωb in the equation (49) and we obtain for z>0

(50)

where Δ(r, θ) can be expressed as

(51)

Since t is large (t>>1), the integral from the expression (50) along Re(ω) = σ is exponentially small and it can be approximated by

(52)

where the constants are

(53)
(54)
(55)
(56)

Taking r = ρ2 and dr = 2ρdρ, we can write (52) as

(57)

using a new variable x = ρC, the former expression can be written as

(58)

where using expressions (40), (53) and (55) we can define

(59)

1. k<<kcr+ and k>>kcr+ cases

Let us consider the limits k<<kcr+ and k>>kcr+. Taking account (59), in these cases it is easy to show that B2+(Ax/C)2B2, due to the fact that (A/C)20 as k0 and k. In this case the expression (58) changes to

(60)

Since we are looking at wk(z, t) for large t values, equation (60) can be simplified by Taylor expansion of the trigonometric functions at x = 0, and written as

(61)

This definite integral may be calculated easily

(62)

If we combine equations (56) and (62), we obtain

(63)

We conclude that asymptotic behavior of the perturbation for k<<kcr and k>>kcr+ is not that a pure exponential, for t>>1, growing with the time as

(64)

2. k=kcr and k=kcr+ cases

Finally, let us consider the asymptotic behavior at the limits k=kcr and k=kcr+. Here, due to the fact that at these critical wavenumbers, as said above, the poles and the branch points collapse with each other (ωp = ωb), the dispersion relation (34) becomes

(65)

Consequently, expression (57) may be expressed as

(66)

By taking a new variable x = ρC, expression (66) here becomes

(67)

Approximating cos(x) to unity at x = 0, (67) turns into

(68)

After evaluating the last integral, we obtain

(69)

Then, combining equations (69) and (56) we can write

(70)

Therefore, the temporal evolution of the disturbance for t>>1, k=kcr and k=kcr+ is given by

(71)

A numerical analysis is needed in order to study a more accurate all-time response of the system. Firstly, we use the dimensionless form of the expression (49) as

(72)

For convenience, we write the normalized w¯k(z̃,ω̃) function for z̃>0 as

(73)

being ω̃b2 the square normalized branch point

(74)

1. k̃cr<k̃<k̃cr+ case

In order to obtain w̃k(z̃,t̃), we calculate the integral (72) inverting the Laplace Transform by the Fourier-Euler (FE) inversion method. Now, we are going to analyze the evolution of the perturbation in different cases.

In Fig. 9 we have plotted w̃k, according to the equations (72) and (73), as a function of time at z̃=7 for At=0.7, τ = 0.01, k̃=6, Γ̃o=1 and z1̃=104 where the discrete spectrum dominates (see Fig. 3). For t̃>>1, we can approximate expression (72) by its discrete spectrum contribution and fit the signal for large time values with a function such as w̃k*(z̃=7,t̃)=ãeω̃pt̃. The constant ã determined from the least square fitting is ã=5.4755104 when ωp̃=1.6255. We obtain the following fitted curve which gives the contribution of the discrete mode for t̃>>1

(75)

In Fig. 9, we can observe an excellent matching for t̃>>1 with the expression given by the equation (75).

FIG. 9.

Plot of normalized w̃k as a function of time for At=0.7, τ = 0.01, k̃=6, Γ̃o=1, z̃=7 and z1̃=104. The solid line is a fit of the form Aeω̃pt̃. The perturbed velocity w̃k is represented on a logarithmic scale. Here is both plotted the continuum and the discrete contributions.

FIG. 9.

Plot of normalized w̃k as a function of time for At=0.7, τ = 0.01, k̃=6, Γ̃o=1, z̃=7 and z1̃=104. The solid line is a fit of the form Aeω̃pt̃. The perturbed velocity w̃k is represented on a logarithmic scale. Here is both plotted the continuum and the discrete contributions.

Close modal

Taking account of expression (75), it is easy to obtain w̃k*(z̃=7,t̃) for all-time range

(76)

In order to study how the continuum affects the disturbance evolution, we plot in Fig. 10 the ratio w̃k/w̃k* against normalized time t̃ for z̃=103,1,5,10. Thus, we can compare our numerical finding w̃k with the disturbance evolution due to the discrete mode w̃k*. As it can be seen from Fig. 10, the dotted line is the asymptote to which w̃k/w̃k* tends as t̃, i.e., in this asymptotic regime, w̃kw̃k*. Firstly, for z̃=103, we find w̃k>w̃k* on the interval (0,t̃c) being t̃c=5.99. This result indicates that the full perturbation w̃k(z̃=103,t̃) grows faster than the mode w̃k*(z̃=103,t̃). On this interval, the continuum modes reinforce the effect of the discrete mode giving an amplified response.

FIG. 10.

Ratio between the full normalized perturbed velocity in z̃ and the corresponding to the discrete mode for At=0.7, τ = 0.01, k̃=6, Γ̃o=1, z1̃=104 and z̃=103,1,5,10. The dashed line is the asymptote where w̃kw̃k*.

FIG. 10.

Ratio between the full normalized perturbed velocity in z̃ and the corresponding to the discrete mode for At=0.7, τ = 0.01, k̃=6, Γ̃o=1, z1̃=104 and z̃=103,1,5,10. The dashed line is the asymptote where w̃kw̃k*.

Close modal

Secondly, let us consider how the continuum affects to w̃k(z̃,t̃) for larger values of z̃. When z̃=1,5 and 10, due to the action of the continuum, we can see that the complete disturbance w̃k grows slower than the discrete mode given by w̃k*, i.e. the response has been reduced by 4%, 32% and 43% respectively (average percentage). Thus, according to Fig. 10, we can establish that when z̃>1 the role of the continuous spectrum is to diminish the disturbance response. This result is not surprising if we bear in mind that the contribution from continuous spectrum is always present.

These observations provide evidence that suggests that the continuous spectrum should not be neglected when studying the temporal evolution of the disturbance.

2. k̃k̃cr and k̃k̃cr+ case

Finally, for the case At=0.7, τ = 0.01, Γ̃o=1, z̃=103 and z̃1=104, it is illustrative to study the asymptotic perturbation evolution as the continuum only exists k̃k̃cr+ and k̃k̃cr being k̃cr=0.71822 and k̃cr+=8.93074. We have fitted the data for t̃>>1 with a function of the form at̃beω̃bt̃ where ω̃b is the normalized branch point. Table I illustrates several values of the fitting parameters a, b and ω̃b corresponding to the expression at̃beω̃bt̃ for different normalized wave numbers. It is worth noting that for k̃>k̃cr+ the exponent b increases with k̃. Firstly, we consider k̃>>k̃cr+ and choose k̃=30. Using expression (72), the evolution of the disturbance due to the continuum for t̃>>1 is given by

(77)
TABLE I.

Fitting parameters a, b and ω̃b in the expression w̃k(z̃=103,t̃)=at̃beω̃bt̃ for different values of k̃, Γ̃o=1, At = 0.7, τ = 0.01, z̃=103 and z1̃=104.

k̃abω̃b
0.10000 0.2206610 1.494522 0.196116 
0.20000 0.1818642 1.479673 0.371391 
0.30000 0.1808763 1.434328 0.514496 
0.40000 0.1906886 1.352626 0.624695 
0.50000 0.1780879 1.201882 0.707107 
0.71822(k̃cr) 0.0571974 0.505777 0.820709 
8.93074(k̃cr+) 0.6207830 0.494106 0.998437 
9.00000 1.3890063 0.907706 0.998460 
12.0000 0.0109028 1.424671 0.999133 
15.0000 0.0019315 1.463316 0.999445 
20.0000 0.0003093 1.463883 0.999688 
30.0000 0.0000388 1.476592 0.999861 
60.0000 0.0000029 1.498662 0.999981 
k̃abω̃b
0.10000 0.2206610 1.494522 0.196116 
0.20000 0.1818642 1.479673 0.371391 
0.30000 0.1808763 1.434328 0.514496 
0.40000 0.1906886 1.352626 0.624695 
0.50000 0.1780879 1.201882 0.707107 
0.71822(k̃cr) 0.0571974 0.505777 0.820709 
8.93074(k̃cr+) 0.6207830 0.494106 0.998437 
9.00000 1.3890063 0.907706 0.998460 
12.0000 0.0109028 1.424671 0.999133 
15.0000 0.0019315 1.463316 0.999445 
20.0000 0.0003093 1.463883 0.999688 
30.0000 0.0000388 1.476592 0.999861 
60.0000 0.0000029 1.498662 0.999981 

We can observe that expression (77) has a good agreement with the asymptotic approximation given by Eq. (64) being the exponent b in equation (77) equal to b=1.4993103/2.

Secondly, for k̃<<k̃cr+, for instance, k̃=0.1, the fit yields

(78)

Thirdly, for k̃=k̃cr=0.71822, we obtain

(79)

And finally, for k̃=k̃cr+=8.93074 we can write

(80)

It is remarkable that expressions (79) and (80) are both in good agreement with asymptotic analytical expression (71) due to the dependence of the perturbed velocity in z, approximately, as w̃kt̃0.5 (here the b exponent in equations (79) and (80) is b0.5).

Finally, in order to highlight the effects of the continuum as t̃0, it is advisable to define the following parameter δ as

(81)

which is a measurement of the contribution of the continuum to the whole response. Thus, the dependence of δ on z̃ for τ = 10−4, Γ̃o=1, k̃=12, z1̃=104 and different Atwood numbers is shown in Fig.11. As can be observed in this figure, for a given value of At, as z̃ increases, δ also increases approaching unity for large values of z̃. As z̃ increases, the contribution of the continuous spectrum dominates completely. Note that when z̃ is increased, the saturation of δ is reached earlier for small Atwood numbers. However, note that for low values of z̃ the contribution of the continuum is negligible, being this effect more accused for increasing At values. We can conclude that, for a given z̃, the effect of the continuous spectrum grows for decreasing values of the At number.

FIG. 11.

Dependence of δ on z̃ for τ = 10−3, Γ̃o=1, k̃=12, z1̃=104 and different Atwood number At = 0.3, 0.5, 0.7 and 0.9. The dashed line is the asymptote to which δ1.

FIG. 11.

Dependence of δ on z̃ for τ = 10−3, Γ̃o=1, k̃=12, z1̃=104 and different Atwood number At = 0.3, 0.5, 0.7 and 0.9. The dashed line is the asymptote to which δ1.

Close modal

It is also of some interest here to examine δ when At, τ and z1̃ are kept fixed (At=0.7, Γ̃o=1, τ = 10−3 and z1̃=104), where δ decreases as z̃ is reduced, see Fig. 12. The minimum value of δ is always reached at k̃=16.7156 which corresponds to the maximum value of the pole ω̃p=2.82214 as expected.

FIG. 12.

Dependence of δ on k̃ for τ = 10−3, z1̃=104, At = 0.7, Γ̃o=1 and z̃=1, 5, 10, and 15.

FIG. 12.

Dependence of δ on k̃ for τ = 10−3, z1̃=104, At = 0.7, Γ̃o=1 and z̃=1, 5, 10, and 15.

Close modal

A linear Rayleigh-Taylor instability for two incompressible immiscible fluids with surface tension, evolving with a free interface in the presence of a uniform gravitational field has been studied. The problem keeps the framework presented in Ref. 12, but instead of using a normal-mode analysis in the frequency domain, a Laplace transform approach is employed to solve the corresponding initial-value problem. When initial conditions corresponding to the Green’s function of the differential equation of the problem are chosen, branch-points are found on the real axis of the complex growth rate-plane. This approach has the potential advantage of including a new part of the spectrum due to branch cuts in the complex plane, not treated by Obied Allah,12 and known as the continuum modes.

In every branch, due to the interaction of surface mode with localized continuous spectrum, two terminal points k̃cr and k̃cr+ determine the interval where the dispersion relation exists. Therefore, due to the existence of the continuous spectrum, the role of surface tension here is to transform surface RTI modes (discrete spectrum) into internal RTI modes at a critical wavenumber k̃cr+.

For the density profile studied, the nonexistence of cut-off wavenumber (wave number corresponding to zero growth rate) is observed. Consequently, our findings seem to refute the common belief that, due to the existence of surface tension, the growth instability for the RTI surface mode goes to zero for the cut-off wavenumber kc, as previous researchers used to believe.

Finally, due to the continuum spectrum, the behavior of the perturbed system is modified when compared to the dominant part of the discrete spectrum. We confirmed that in the absence of discrete spectrum, an initial perturbation grows by the effect of the continuous spectrum as an exponential times a time negative fractional power for large time values.

1.
Lord
Rayleigh
,
Scientific Papers
(
Dover
,
New York
,
1965
), Vol. 2.
2.
G. I.
Taylor
, “
Instability of liquid surfaces when accelerated in a direction perpendicular to their planes
,”
Proc. R. Soc. London, Ser. A
201
,
192
(
1950
).
3.
W. C.
Elmore
and
W. C.
Heald
,
Physics of waves
(
Dover
,
New York
,
1985
).
4.
S.
Chandrasekhar
,
Hydrodynamic and Hydromagnetic Stability
(
Dover
,
New York
,
1981
).
5.
W. H.
Reid
, “
The effect of surface tension and viscosity on the stability of two superposed fluids
,”
Proc. Camb. Philos. Soc.
57
,
415
(
1961
).
6.
B. J.
Daly
, “
Numerical study of the effect of surface tension on interface instability
,”
Phys. Fluids
12
,
1340
(
1969
).
7.
D. S.
Vaghela
and
R. K.
Chhajlani
, “
Rayleigh-Taylor instability of magnetized partially-ionized superposed fluids with rotation and surface tension in porous medium
,”
Astrophys. Space Sci.
149
,
301
(
1988
).
8.
K. O.
Mikaelian
, “
Rayleigh-Taylor and Richtmyer-Meshkov instabilities in multilayer fluids with surface tension
,”
Phys. Rev. A
42
,
7211
(
1990
).
9.
N. F.
El-Ansary
,
G. A.
Hoshoudy
,
A. S.
Abd-Elrady
, and
A. H. A.
Ayyad
,
Phys. Chem. Chem. Phys.
4
,
1464
(
2002
).
10.
M.
Chertkov
,
I.
Kolokolov
, and
V.
Lebedev
,
Phys. Rev. E
71
,
055301
(
2005
).
11.
P. K.
Sharma
,
R. P.
Prajapati
, and
R. K.
Chhajlani
, “
Effect of surface tension and rotation on Rayleigh-Taylor instability of two superposed fluids with suspended particles
,”
Acta Physica Polonica A
118
,
576
(
2010
).
12.
M. H.
Obied Allah
, “
Effects of stratification, surface tension and rigid planes on Rayleigh-Taylor instability
,”
Ind. J. Pure Appl. Math.
32
(
3
),
303
(
2001
).
13.
K. M.
Case
, “
Hydrodynamic stability and the inviscid limit
,”
J. of Fluid Mech.
10
,
420
(
1961
).
14.
K. M.
Case
, “
Taylor instability of an inverted atmosphere
,”
Phys. Fluids
3
,
366
(
1960
).
15.
K. M.
Case
, “
Stability of an idealized atmosphere
,”
Phys. Fluids
3
,
149
(
1960
).
16.
A.
Burger
, “
Instability associated with the continuous spectrum in a baroclinic flow
,”
J. Atmos. Sci.
23
(
1966
).
17.
Z.
Sedlácêk
, “
Electrostatic oscillations in cold inhomogeneous plasma I. Differential equation approach
,”
J. of Plasma Phys.
5
,
239
(
1971
).
18.
E.
Ott
and
D. A.
Russell
, “
Diffuse-boundary Rayleigh-Taylor instability
,”
Phys. Rev. Lett.
41
,
1048
(
1978
).
19.
R.
Menikoff
,
R. C.
Mjolsness
,
D. H.
Sharp
,
C.
Zemach
, and
B. J.
Doyle
, “
Initial value problem for Rayleigh-Taylor instability of viscous fluids
,”
Phys. Fluids
21
,
1674
(
1978
).
20.
L. H.
Gustavsson
, “
Initial value problem for boundary layer flows
,”
Phys. Fluids
22
,
1602
(
1979
).
21.
D. A.
Russell
and
E.
Ott
, “
The linear theory of the Rayleigh-Taylor instability in the equatorial ionosphere
,”
J. Geophys. Res.
84
,
6573
, (
1979
).
22.
T. J.
Bogdan
and
P. S.
Cally
, “
Waves in magnetized polytropes
,”
Proc. R. Soc. Lond. A, Math. Phys. Eng. Sci.
453
(
1960
),
943
(
1997
).
23.
F. B.
Hildebrand
(
Dover
,
New York
,
1992
).