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.
I. INTRODUCTION
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 , 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 , with k real, possess both a discrete set of complex eigenfrequencies ωn, , 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.
II. LINEAR STABILITY ANALYSIS
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
where ρ1 and ρ2 are independent constants, being , 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
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
where u, δ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 is the perturbation of the interfacial surface and the subscript s indicates the interfacial surface.
Introducing Fourier transforms23 with respect to x and y,
and taking account the influence of the surface tension,4 we obtain
Here, δ(z − zs) denotes the Dirac function and Ts is the surface tension of interfacial surface. Moreover, and ω are both the square wavenumber and the growth rate of the harmonic perturbations.
After eliminating δpk, δρk, uk and from the above equations, an equation for is obtained as
Next, we consider the IVP and introduce the Laplace transform.23 So we define the forward Laplace transform in the complex ω -plane as
Here we will use the second derivative property of the Laplace Transform
being F1 and F2 the initial conditions that are imposed at t=0
So, application of this transformation to equation (17) gives
where H(z; ω, k) term is given by
and the term involves all initial conditions being
This differential equation (21) is now to be solved for , 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
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 as . Additional conditions must be satisfied on the horizontal plane (z = zs) on which there are density discontinuities. Evidently, one condition is that 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 . In virtue of the continuity of across z = zs and the boundedness of ρo this limiting process leads to (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.
III. DISPERSION RELATION
Concerning to the initial condition given by equation (24), must be continuous at z = −z1. Besides, the first derivative of 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
Taking account the above conditions, the general solution of (21) is a linear combination of two linearly independent solutions. Thus, becomes
where α is given by equation being
where 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
For simplicity, we use the solution for
where the dispersion function Δ(k, ω) is given by
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
Consequently, has two poles , 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 , the dispersion relation is
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 . Thus, expression (33) is the basis of our investigation. Following to Ott and Russell,18 due to the fact that the function does not turn out to be an even function of α, this function has branch points at (α = 0) on the real axis associated with a continuous spectrum, where
If the ω -plane is cut along the interval connecting - ωb and ωb, 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 to obtain (z, t) as function of time. The inverse Laplace transform of is given by
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 in the region . 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 .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 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, is then a single-valued analytic function in the region inside the lines Γ + γ, C1, C2 and Cb. In accordance with the boundary condition , the square root given by (28) is defined, so that for . Thus, this integral (41) in the ω- plane can now be written as
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).
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).
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
can be expressed as a linear combination of terms of the type eωt where ω is real and runs between . 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
Let us examine the impact of the continuous spectrum on the discrete spectrum.
IV. NUMERICAL RESULTS
It is always useful to convert the expression (33) of for to a dimensionless form. So, the most relevant dimensionless parameters become
On the other hand, we define the Atwood number as
where is the density jump. Hence, for may be written in the normalized form
where the denominator of (47) is the normalized dispersion relation
Fig. 3 shows the normalized square growth rate against the normalized wavenumber obtained when in equation (48) for At = 0.4 and At = 0.7 and τ = 0.01 and τ = 0.02. The boundary of the region of the -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 ) in each plot. The solid lines stand for the normalized square growth rate of the poles (discrete spectrum) . As it can be appreciated, in every branch, there are two terminal points with critical wavenumbers and , respectively, that limit the dispersion relation domain. As can be observed, for or the (poles) move towards (branch points). Due to role of surface tension here, surface RTI modes (discrete spectrum) have a maximum and transform into internal RTI modes at . This is possible when the continuous spectrum, neglected by Obied,12 is now taken into account. Thus, for , then 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 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 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 and depend on the density jump, indicating the sensitivity of surface mode to At changes.
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 and values are only indicated for τ = 0.02.
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 and values are only indicated for τ = 0.02.
At the secondary level, we estimate numerical values for , which are easily derived by combining equations (35) and (40). As depends on At and τ, we studied this dependence with At for a fixed τ. Thus, in Fig. 4 we have represented the critical wavenumber 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 and , 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 and 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 . Thus, for instance, for τ = 0.01 the terminal wavenumbers are and , being the terminal Atwood number.
Dimensionless critical wave number as a function of the dimensionless parameter At for different τ dimensionless parameter. The lower and the upper branch correspond to and respectively.
Dimensionless critical wave number as a function of the dimensionless parameter At for different τ dimensionless parameter. The lower and the upper branch correspond to and respectively.
Taking into account the previous analysis, the square of the dimensionless growth rate against the dimensionless wavenumber close to was plotted for τ = 0.01 and two different Atwood numbers, see 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.
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.
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 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 for At = 0.7, τ = 0.01, , , and , see Fig. 6.
Contour plot of normalized imaginary part of versus normalized growth rate, for At = 0.7, , , , τ = 0.01 and . It illustrates the branch cut running along real axis between Re =-1 and Re =1. Each pair of lobes represents a pole.
Contour plot of normalized imaginary part of versus normalized growth rate, for At = 0.7, , , , τ = 0.01 and . It illustrates the branch cut running along real axis between Re =-1 and Re =1. Each pair of lobes represents a pole.
This clearly illustrates the continuous spectrum through the branch cut running along real axis between and . On the other hand, poles of function can be appreciated due to the existence of pairs of lobes at
Now, we plot the contour plot of the imaginary part of versus normalized growth rate for At = 0.7, τ = 0.01, , , at , (see Fig. 7). In this case , as , there is not a discrete spectrum and only the branch cut and appears.
Contour plot of normalized imaginary part of function versus normalized growth rate for At = 0.7, , , , τ = 0.01 and .
Contour plot of normalized imaginary part of function versus normalized growth rate for At = 0.7, , , , τ = 0.01 and .
V. TEMPORAL EVOLUTION OF DISTURBANCE
Our goal in this section is to analyze the time evolution of a perturbation and observe the transient response of the system.
A. Asymptotic analysis for large t (t1)
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 and and, secondly, for and . Consequently, the only contribution to comes from the branch point at ω = ωb. To calculate the asymptotic behavior of 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 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 , except the one at ω = ωb.
To obtain the time-asymptotic solution , we require t to be so large that the integral along 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 far from ωb decrease as .The asymptotic solution can be expressed as
Using the expressions (28),(40) and (33) we replace ω − ωb by in the integrand of (49). In the vicinity of the branch point, we replace ω by ωb in the equation (49) and we obtain for
where Δ(r, ) can be expressed as
Since t is large (), the integral from the expression (50) along Re(ω) = σ′ is exponentially small and it can be approximated by
where the constants are
Taking r = ρ2 and dr = 2ρdρ, we can write (52) as
using a new variable x = ρC, the former expression can be written as
1. and cases
Let us consider the limits and . Taking account (59), in these cases it is easy to show that , due to the fact that as and . In this case the expression (58) changes to
Since we are looking at (z, t) for large t values, equation (60) can be simplified by Taylor expansion of the trigonometric functions at x = 0, and written as
This definite integral may be calculated easily
We conclude that asymptotic behavior of the perturbation for and is not that a pure exponential, for , growing with the time as
2. and cases
Finally, let us consider the asymptotic behavior at the limits and . 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
Consequently, expression (57) may be expressed as
By taking a new variable x = ρC, expression (66) here becomes
Approximating to unity at x = 0, (67) turns into
After evaluating the last integral, we obtain
Therefore, the temporal evolution of the disturbance for , and is given by
B. Numerical solution
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
For convenience, we write the normalized function for as
being the square normalized branch point
1. case
In order to obtain , 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 , according to the equations (72) and (73), as a function of time at for At=0.7, τ = 0.01, , and where the discrete spectrum dominates (see Fig. 3). For , we can approximate expression (72) by its discrete spectrum contribution and fit the signal for large time values with a function such as . The constant determined from the least square fitting is when . We obtain the following fitted curve which gives the contribution of the discrete mode for
Plot of normalized as a function of time for At=0.7, τ = 0.01, , , and . The solid line is a fit of the form . The perturbed velocity is represented on a logarithmic scale. Here is both plotted the continuum and the discrete contributions.
Plot of normalized as a function of time for At=0.7, τ = 0.01, , , and . The solid line is a fit of the form . The perturbed velocity is represented on a logarithmic scale. Here is both plotted the continuum and the discrete contributions.
Taking account of expression (75), it is easy to obtain for all-time range
In order to study how the continuum affects the disturbance evolution, we plot in Fig. 10 the ratio against normalized time for . Thus, we can compare our numerical finding with the disturbance evolution due to the discrete mode . As it can be seen from Fig. 10, the dotted line is the asymptote to which tends as , i.e., in this asymptotic regime, . Firstly, for , we find on the interval being . This result indicates that the full perturbation grows faster than the mode . On this interval, the continuum modes reinforce the effect of the discrete mode giving an amplified response.
Ratio between the full normalized perturbed velocity in and the corresponding to the discrete mode for At=0.7, τ = 0.01, , , and . The dashed line is the asymptote where .
Ratio between the full normalized perturbed velocity in and the corresponding to the discrete mode for At=0.7, τ = 0.01, , , and . The dashed line is the asymptote where .
Secondly, let us consider how the continuum affects to for larger values of . When and 10, due to the action of the continuum, we can see that the complete disturbance grows slower than the discrete mode given by , 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 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. and case
Finally, for the case At=0.7, τ = 0.01, , and , it is illustrative to study the asymptotic perturbation evolution as the continuum only exists and being and . We have fitted the data for with a function of the form where is the normalized branch point. Table I illustrates several values of the fitting parameters a, b and corresponding to the expression for different normalized wave numbers. It is worth noting that for the exponent b increases with . Firstly, we consider and choose . Using expression (72), the evolution of the disturbance due to the continuum for is given by
Fitting parameters a, b and in the expression for different values of , , At = 0.7, τ = 0.01, and .
. | a . | 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.0571974 | 0.505777 | 0.820709 | |
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 |
. | a . | 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.0571974 | 0.505777 | 0.820709 | |
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 .
Secondly, for , for instance, , the fit yields
Thirdly, for , we obtain
And finally, for we can write
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 (here the b exponent in equations (79) and (80) is ).
Finally, in order to highlight the effects of the continuum as , it is advisable to define the following parameter δ as
which is a measurement of the contribution of the continuum to the whole response. Thus, the dependence of δ on for τ = 10−4, , , and different Atwood numbers is shown in Fig.11. As can be observed in this figure, for a given value of At, as increases, δ also increases approaching unity for large values of . As increases, the contribution of the continuous spectrum dominates completely. Note that when is increased, the saturation of δ is reached earlier for small Atwood numbers. However, note that for low values of the contribution of the continuum is negligible, being this effect more accused for increasing At values. We can conclude that, for a given , the effect of the continuous spectrum grows for decreasing values of the At number.
Dependence of δ on for τ = 10−3, , , and different Atwood number At = 0.3, 0.5, 0.7 and 0.9. The dashed line is the asymptote to which .
Dependence of δ on for τ = 10−3, , , and different Atwood number At = 0.3, 0.5, 0.7 and 0.9. The dashed line is the asymptote to which .
It is also of some interest here to examine δ when At, τ and are kept fixed (At=0.7, , τ = 10−3 and ), where δ decreases as is reduced, see Fig. 12. The minimum value of δ is always reached at which corresponds to the maximum value of the pole as expected.
VI. SUMMARY AND CONCLUSIONS
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 and 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 .
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.