Pure single-mode Rayleigh-Taylor instability for arbitrary Atwood numbers

The validity of theoretical investigation on Rayleigh-Taylor instability (RTI) with nonlinearity is quite important, especially for the simplest and the commonest case of a pure single-mode RTI, while its previous explicit solution in weakly nonlinear scheme is found to have several defections. In this paper, this RTI is strictly solved by the method of the potential functions up to the third order at the weakly nonlinear stage for arbitrary Atwood numbers. It is found that the potential solution includes terms of both the stimulating and inhibiting RTI, while the terms of the decreasing RTI are omitted in the classical solution of the weakly nonlinear scheme, resulting in a big difference between these two results. For the pure single-mode cosine perturbation, comparisons among the classical result, the present potential result and numerical simulations, in which the two dimensional Euler equations are used, are carefully performed. Our result is in a better agreement with the numerical simulations than the classical one before the saturation time. To avoid the tedious expressions and improve a larger valid range of the solution, the method of the Taylor expansion is employed and the velocities of the bubble and spike are, respectively, obtained. Comparisons between the improved and the simulation results show that the improved theory can better predict the evolution of the interface from the linear to weakly nonlinear, even to later of the nonlinear stages.

For a simple case of a single cosine-mode perturbation at the interface between the heavier fluid overlapping the lighter one, i.e., with ε the amplitude, k = 2π∕λ wave number, λ the wavelength of the initial perturbation and ε ≪ λ, this interface is prone to RTI. This interface initially develops with time at growth rate γ = Agk , where A = (ρ h − ρ l )∕(ρ h + ρ l ) is Atwood number, g is gravitational acceleration and ρ h (ρ l ) is density of the heavier (lighter) fluid. That is, the perturbed interface evolves in the form classical L in which η ε γ = t exp( ) L is called as the linear amplitude of the interface. This is known as the linear theory, see both the Rayleigh and Taylor papers 1,2 . After a short time, this evolution interface includes the second, third harmonics, and so on, because of the mode-coupling effect in a nonlinear system. The appearance of higher harmonics results in the different amplitude growth of the interface in the heavier and lighter fluids. The portion of the interface entering the heavier (lighter) fluid is named as a bubble (spike). Generally, the spike moves faster than the bubble, especially for the larger Atwood number. The higher harmonics play an important role for the amplitude difference between the bubble and spike. The evolution amplitudes of the first three harmonics are given in weakly nonlinear theory [22][23][24][25] , which are as below: where the amplitude of the fundamental mode, classical 1 η , includes the linear one η L and the third-order feedback/ correction. Within the third-order weakly nonlinear framework, the second and the third harmonics have no feedbacks from higher orders. In the recent work 26 , by employing the tracking technology of the mode coupling up to the fifth order, the effect of the mode-coupling branches/channels on harmonic amplitudes in single-mode classical Rayleigh-Taylor instability for arbitrary Atwood numbers in the weakly nonlinear stage is investigated.
Because the bubble (spike) is the structure which moves the fastest to the heavier (lighter) fluid, its evolution is greatly focused on. For the asymptotical behavior, the expansion method of the potential function was first employed by Layzer to bubbles 27 , and then bubbles and spikes [28][29][30] . Based on the weakly nonlinear solution, the perturbation solutions can lead to a perfect Padé approximation by matching the linear solution and the asymptotic solution [31][32][33] .
However, when the above theory of the third-order weakly nonlinearity is deeply investigated, deficiencies are found, resulting in the corresponding prediction largely deviating results from the numerical simulations or experiments. Within the third-order weakly nonlinear framework, the interface equation should be When the amplitudes of the first three harmonics are replaced with the classical ones (3a-c) the two deficiencies are found as below.
(i) the initial velocity of the whole interface is not equal to zero, while the general RTI starts from the initial interface with the rest; (ii) the interface is not in the form of the pure single-mode perturbation, but multi-mode perturbations, and then the initial interface condition (1) is not strictly satisfied.
When the linear amplitude is carefully taken into account, it should be which was also obtained by Forbes 34 . In Eq. (5), these two terms play a different role on RTI: the former stimulates the RTI, while the latter inhibits the RTI. Generally speaking, just the stimulating term is focused on, while the inhibiting term is neglected. As a result, the prediction of the linear amplitude without the inhibiting term is always higher than that with this term. For a nonlinear system, due to the effect of the nonlinear mode coupling, the neglect of the inhibiting term in the first order can make not only the inhibiting terms but also the stimulating ones disappear. That is, in the weakly nonlinear stage, the disappearance terms including the inhibiting and stimulating ones can make the prediction have a larger difference from the related results of the numerical simulations or experiments. Meanwhile, this classical solution of the weakly nonlinear stage can mislead the readers and the researchers for better understanding the weakly nonlinearity, even the whole nonlinearity. In this paper, the evolution equations of the interface are strictly solved, and the solutions satisfying the condition of the initial velocity, together with the pure single-mode perturbation, are obtained order by order. Therefore, this work raised the linearized solution, while the previous works neglected the role related to the inhibiting terms ( γ − t exp[ ]), resulting to a big difference between the third theory of the weakly nonlinearity and the numerical simulations. Based on the pure single-mode RTI, our work presents the third-order theory of the weakly nonlinearity with the role of the whole terms, and predicts the better results from the initial linear stage to the later stage of the nonlinearity.

Theory Framework and Explicit Solutions
When two incompressible, inviscid and irrotational fluids with arbitrary density rations (i.e., arbitrary Atwood numbers) in two dimensions are taken into account, the governing equations are When the heavier fluid is overlapping the lighter one in the gravitational field of acceleration g ge y → = − → , any perturbations at the interface can arouse the rest interface. Here, we consider the pure single-mode cosine perturbation in the form of According to the mode coupling rule, the evolution interface, within the three-order framework, is given by Eq. ((4)) where the amplitude of the fundamental mode η 1 (t) = η 1,1 (t) + η 3,1 (t). The velocity potentials of the fluids, satisfying their Laplace equations and infinity conditions, can be expressed as The time coefficients of the evolution interface need to be solved at last. The solving steps are listed as: (i) substitute Eqs. (4), (9) and (10)  (v) construct an ordinary differential equation for the first-order, second-order, and third-order terms, respectively, consisting of the coefficients of ε, ε 2 , and ε 3 ; (vi) eliminate the unknown coefficients of the velocity potential, and obtain the corresponding deference equations on amplitude coefficients of the evolution interface.
A second-order ordinary differential equation for the first-order terms together with the initial conditions is given as below: where the double dots denote the second derivative of the time and the single dot denotes the first derivative of time. The solution is obtained as 1,1 As to differential equations for the second and the third terms, the initial values of the amplitude and velocity, The initial conditions are taken into account, these equations can be solved in turn. The amplitudes of the first three harmonics are expressed as bellow.
classical p resent  As mentioned above, from the results one finds that there are several terms which either stimulate or inhibit the RTI, especially for the solutions from the third order. These solutions consist of hyperbolic sine and cosine functions including stimulating and inhibiting terms. That is to say, if the inhibiting term in the linear solution is deleted, many terms will disappear.

Comparison between the Results of the Theory and the Numerical Simulations
In order to validate the present theory, the results of the classical and the present theories are compared with the numerical simulation. The code, based on CFD 2 , used here is developed by our research team. For the two dimensional Euler equations of the inviscid and compressible fluid dynamics, fifth-order finite difference WENO schemes with the third-order Runge-kutta time discretization are used. In our numerical simulations, we set up the problem, in which the related parameters step from the ref. 35  For the whole field of the fluids, the initial velocity in the x-or y-direction is zero. Because of the RTI, the interface perturbation will grow with time. The top of the bubble (spike) is located at x = 0(x = 0.125) and the corresponding amplitude or velocity can be obtained by tracking its location at different time.  www.nature.com/scientificreports www.nature.com/scientificreports/ For the whole evolution interface in the RTI, the top of the bubble or the spike is the structure which grows fastest and has the longest amplitude. In our explicit interface (4), the tops of the bubble and the spike are located at x = 0 and x = π∕k, and their amplitudes are, respectively, η bb (t) = η(x = 0, t) and η sp (t) = |η(x = π/k, t)|. Here, the positive amplitudes are taken into our account, and then their velocities are, respectively, v t . Due to the tedious expressions of the amplitude and the velocity of the bubble and spike, they are not explicitly given here.
In order to determine the value of the mesh step h, some results of the mesh convergence analysis are shown in Fig. 1. From comparisons of the dimensionless amplitudes of the bubbles and the spikes among the mesh steps 1/1250, 1/2500, 1/5000 and 1/6000 in Fig. 1, the mesh step is set as the uniform h = 2 × 10 −4 in our following computations.
The amplitude and velocity of the bubble and spike from the classical theory (denoted by circles), the present potential theory (denoted by dots) and numerical simulation (denoted by stars) are, respectively, compared in Figs. 2 and 3 for Atwood number A = 0.2 (left) and A = 0.8 (right). Here, the perturbation wave number k and gravitational acceleration g are used to normalized the time, amplitude and velocity.
From Fig. 2, one finds that in weakly nonlinear stage, the prediction of the present potential theory is in good agreement with the results of the numerical simulation, while the amplitude predicted by the classical theory is always larger than our result. Because of the negative feedback to the fundamental mode from the third order, the amplitude of the the bubble first increases up to its maximum at the saturation time, and then decreases with time for A = 0.2 and A = 0.8. It is obvious that the saturation time predicted by the classical theory is smaller than ours. That is to say, our theory predicts the range of the valid time longer than the classical. For A = 0.8 spike, the feedback to the fundamental mode from the third order is not negative, but positive, resulting in the amplitude of the spike growing all the time. Anyway, it is found that the present potential theory gives a better agreement with the numerical simulation for longer time than the classical.
As stated above, the classical theory did not predict the RTI with the initial velocity being zero, as shown in Fig. 3. It is obvious that the velocity predicted by the classical theory is larger than ours. Our results agree better with the numerical simulations for A = 0.2 and A = 0.8.
Considering the tedious expressions of the amplitude and the velocity of the bubble and spike, together with the larger range of the valid time in the weakly nonlinear stage predicted by present potential theory, the with initial amplitude σ = kε and  = v v k g / . For the sake of checking whether above velocity formulas are effective, the improved velocities predicted by Eqs. (18a) and (18b) (denoted by pentagons), our potential theory (denoted by dots) and the numerical simulations (denoted by stars) are compared in Fig. 4. It is found that the improved theory predicts a better result than our potential theory. On the one hand, the results of these two theories are in good agreement in the early stage of the weakly nonlinearity. On the other hand, our improved theory extends the valid range of our potential theory, either for the bubble or the spike, especially for the larger Atwood number, as can be affirmed by the good agreement of our potential theory and the numerical simulations in Fig. 4.
In fact, these purely potential-flow solutions only exist for a finite time interval, since they are ultimately killed at a finite time, by the curvature singularity at the interface predicted by Moore 36 . In this paper, the valid time of the improved theory is extended due to two factors. One is the the theory of the weakly nonlinearity containing the whole terms [not only the term γt exp( ) but also the one γ − t exp( )]. The role of the inhibiting terms can not be neglected. The other is the application of Taylor expansion. These two aspects play an important role.

conclusion
In this paper, we employ the method of the small parameter expansion of the potential function with nonlinear corrections up to the third order to analytically explore the history evolution of the bubble and the spike for the pure single-mode cosine perturbation at the interface between two different fluids in classical (irrotational, incompressible, and inviscid fluids) RTI. First, we analysis the defects in the classical results: (i) the initial interface is not in the form of the pure single mode, but in the multi mode; (2) the initial velocity of the interface perturbation is not zero. These two defects are resulted from the absence of the time inhibiting term in the linear www.nature.com/scientificreports www.nature.com/scientificreports/ solution. As a result, we keep the inhibiting term in the linear solution, and obtain the solutions satisfying the differential equations of the evolution interface in potential functions. The comparisons of the amplitude and the velocity evolutions of the bubble and the spike are performed among the classical, our potential theories and our numerical simulations. It is found that before the saturation time, the results of the potential theory agree with the numerical simulations better than the classical one. After the saturation time, the results of the classical and our present theories start decreasing rapidly with time due to the negative feedback to the fundamental mode from the third order. For this reason, we expand our potential results of the velocities of the bubble and the spike in Taylor series in time, and obtain the improved theory. Comparisons among the potential theory, the improved theory and numerical simulations show that the improved theory predicts a better result, either for the velocities of the bubble or the spike, and the valid range of the time is also larger than the potential theory, as can be confirmed in Fig. 4.