Interface coupling effects of weakly nonlinear Rayleigh–Taylor instability with double interfaces
Li Zhiyuan1, Wang Lifeng1, 2, Wu Junfeng2, Ye Wenhua1, 2, †
Institute of Applied Physics and Computational Mathematics, Beijing 100094, China
Center for Applied Physics and Technology, HEDPS, and College of Engineering, Peking University, Beijing 100871, China

 

† Corresponding author. E-mail: ye_wenhua@iapcm.ac.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11575033, 11675026, and 11975053) and the Science Foundation from China Academy of Engineering Physics (Grant No. CX2019033).

Abstract

Taking the Rayleigh–Taylor instability with double interfaces as the research object, the interface coupling effects in the weakly nonlinear regime are studied numerically. The variation of Atwood numbers on the two interfaces and the variation of the thickness between them are taken into consideration. It is shown that, when the Atwood number on the lower interface is small, the amplitude of perturbation growth on the lower interface is positively related with the Atwood number on the upper interface. However, it is negatively related when the Atwood number on the lower interface is large. The above phenomenon is quantitatively studied using an analytical formula and the underlying physical mechanism is presented.

1. Introduction

The Rayleigh–Taylor instability (RTI)[1,2] exists widely in inertial confinement fusion (ICF) and supernova explosions. Understanding the mechanism of the RTI is important to the success of the ignition of ICF[3] and to explain the nonlinear evolution of supernova explosions.[4] Generally, it will happen when a heavy fluid is accelerated by a light one, if some perturbations are on the interface between two fluid layers. Focusing on this kind of RTI, many studies have been done in its linear,[1,2] nonlinear,[5,6] and turbulence mixing[7,8] regimes. However, in some practical applications, the RTI with multiple material interfaces is more attractive. For examples, there are an ablator layer, a deuterium–tritium (DT) fuel layer filled with a low density DT gas in the one-shell ICF targets.[9] Alternatively, double-shell targets, which consist of more material layers, were suggested by Amendt[10] and Canaud.[11] Especially, in the implosion process of ICF, the ablation effect makes the material distribution in the capsule more complicated.[1216] The classic RTI models which only include a single material interface are not good enough for depicting the perturbation growth of the RTI with multiple interfaces. However, to understand the mechanism of the RTI with multiple interfaces, some efforts have been done by researchers. The first series of studies were made by Mikaelian. The author presented an analytical model to describe the temporal evolution of the RTI in the linear growth regime at the interfaces of any number of stratified fluids forming an arbitrary density profile.[17,18] Using this model, some meaningful works have been done by him. For example, the density gradient stabilization effect[19] and the way that how to adjust the perturbation amplitude on the two interfaces to kill the exponentially growing mode.[20] Besides, the author extended the model in the plane geometry to that in the spherical[21] and the cylindrical geometry.[22] Following Mikaelian’s works, the weakly nonlinear growth regime of the RTI was taken into consideration by Wang et al.[23] They presented a weakly nonlinear model of the RTI with a finite-thickness fluid layer. They showed that the weakly nonlinear effect and interface thickness effects are important to understand the flow phenomenons of ICF and supernova explosions. It should be pointed out that this weakly nonlinear model only considers the situation that the Atwood number on the lower interface is 1 and the Atwood number on the upper interface is −1. Focusing on the plane and cylindrical geometry, recently, Guo et al.[24,25] used the potential model to study the linear RTI with multiple fluid layers. Making a summary about the existing studies, one can see that most of studies about the RTI with multiple fluid layers are focussed on in the linear regime. The weakly nonlinear RTI with double interfaces studied by Wang et al.[23] did not cover the situation with varied Atwood numbers. In this study, the situations, the Atwood number on the lower interface is from 0 to 1 and the Atwood number on the upper interface is from −1 to 0, are taken into consideration. The interface coupling effects of the RTI in the weakly nonlinear regime are studied. Euler numerical simulations along with analytical studies are used to do the research in this paper. The results presented in this study are helpful to understand the physical mechanism of the interface coupling effects of the RTI with double interfaces.

The remainder of this article is organized as follows. Firstly, the perturbation interactions between two interfaces are studied by numerical simulations. Secondly, the variation rule of the perturbation amplitude on the lower interface with the Atwood number on the upper interface is studied by the analytical model. Meanwhile, a suppression strategy of the perturbation growth of the RTI on the lower interface is proposed. Thirdly, the underlying mechanisms are explained. Finally, the major conclusions are summarized.

2. Statement of the problem

A two-dimensional (2D) problem in the plane geometry, as shown in Fig. 1, is concerned in this study. Along the direction of acceleration, a finite-thickness fluid layer with density ρ2 is located between two semi-infinite fluid layers with densities ρ3 and ρ1 respectively. These three fluids are all incompressible and immiscible. Apparently, there are two interfaces among these three fluids. Initially, the upper one between fluid layers with densities ρ3 and ρ2 is at the position x = 0. The lower one between fluid layers with densities ρ2 and ρ1 is at the position x = −d. A2 and A1, which are defined as (ρ3ρ2)/(ρ3 + ρ2) and ρ2ρ1)/(ρ2 + ρ1), are the Atwood numbers on the upper and lower interfaces. The situations that A1 > 0 and A2 < 0 are discussed in this study. A perturbation is set on the lower interface. It is a single-mode cosinusoidal perturbation defined as η(x, t = 0) = η0cos(kx), where η0 is the amplitude of the initial perturbation; k = 2π/L is the wavenumber; and L is the wavelength.

Fig. 1. Schematic drawing about the initial condition of the RTI with double interfaces.
3. Numerical simulation

Euler simulations are made in an L × 10L domain containing fluids with three different densities. The domain length in the acceleration direction is 10L. The value of L is 0.01π cm. The acceleration g is 1000 cm/μs2. ρ2 is fixed as 1 g/cm3. In this way, A1 and A2 on the interfaces can be changed by adjusting the values of ρ1 and ρ3. The boundaries at the ends of the y direction are walls. At the ends of x direction, the periodic boundary condition is used. At the initial time, the wavelength on the lower interface is L, and its amplitude is 0.01L. There are 256 uniform grids in the x direction and 2000 grids in the y direction. Most of grids are located in the evolution region of the perturbation. The Euler simulations are done by a well-tested second-order Godunov hydrodynamic simulation code. This code is recently developed in house. The finite volume method and the massively parallel computation are used. The second order HLLE approximation Riemann solver[26] is used. The slope limit is moncen.[27]

To study the interface coupling effects, the temporal evolution of the perturbation growth amplitude on the lower interface (η1) is depicted with some combinations of A1 and A2. Ten situations depicted in Fig. 2 are chosen to show this phenomenon. They are the conditions that A1 = 0.2 and A1 = 0.9 with five different A2. The thickness between two interfaces is 0.004 cm, which means that the normalized thickness kd is 0.5. The computational results are represented by dotted lines. To judge when the perturbation growth reaches the weakly nonlinear regime, the analytical results of the linear RTI model, which is presented by Mikaelian,[18] are also depicted by solid lines in the figure. It is shown that, early in the temporal evolution of the perturbation, the analytical results are closed to the simulation results. This implies that this period of time is the linear RTI regime. At the time when the perturbation growth amplitude η1 reaches 0.2 time of the wavelength, the simulation results deviate from the linear analytical results greatly. This means that the RTI has already reached the weakly nonlinear regime. One can see that, in the weakly nonlinear regime, η1 is positively related with A2 when A1 is small. However, it is negatively related when A1 is large. More Euler simulations have shown that this phenomenon exists in situations with kd ⩾ 0.4. The reason why situations with kd < 0.4 are not concerned is that the single medium Euler simulation can not accurately capture the interfaces when kd is too small. The following studies are made under the condition kd ⩾ 0.4.

Fig. 2. The temporal evolution of the perturbation amplitude calculated by Euler simulation (dotted lines) and linear growth analytical model (solid lines); (a) A1 = 0.2; (b) A1 = 0.9.
4. Variation rule of η1 with A2

To further confirm how the interface coupling effects impact on the correlation between η1 and A2, an analytical research is made in the following. Comparing the simulation results and linear analytical results depicted in Fig. 2, one can see that, at the end of the linear regime, the correlations between η1 and A2 for the two methods are consistent. In this way, we choose the time τa, when the classical single interface RTI with Atwood number A1 reaches the saturation amplitude, to make the analytical research. This time is depicted by blue vertical line in Fig. 2. Results show that, at the time τa, the correlations between τ1 and A2 are consistent with the one in the weakly nonlinear regime. Therefore, the interface coupling effects at the time τa can be used to judge the correlations between τ1 and A2 in the weakly nonlinear regime. According to the linear analytical model,[18] we can obtain the partial derivative ∂η1(τa)/∂A2. Figure 3 depicts its contours with two fixed normalized thicknesses. They are kd = 0.4 and 0.8. All the calculations are done with g = 1000 cm/μs2, L = 0.01 cm, and η0 = 0.01L. From Fig. 3, one can see that the lines ∂η1(τa)/∂A2 = 0 are nearly straight (highlighted by black dotted lines) under these two fixed non-dimensional thicknesses. Below the lines ∂η1(τa)/∂A2 = 0, ∂η1(τa)/∂A2 is almost larger than 0. However, ∂η1(τa)/∂A2 is less than 0 above the lines ∂η1(τa)/∂A2 = 0. More researches with other concerned non-dimensional thicknesses, which are not demonstrated in this study, also show the same phenomenon. It means that under the line ∂η1(τa)/∂A2 = 0, η1 and A2 are positively related. While η1 and A2 are negatively related above the line. The parameter kd only has influences on the variation degree of η1 along with the change of A2. With the analysis of Fig. 2, one has already known that the above interface coupling effects at the time τa are consistent with ones presented in the weakly nonlinear RTI regime. In this way, if the line ∂η1(τa)/∂A2 = 0 is defined, the correlation between η1 and A2 in the weakly nonlinear regime under arbitrary concerned A1 and kd will be clear.

Fig. 3. Contours of ∂η1(τa)/∂A2 with the normalized thickness kd = 0.4 (a) and kd = 0.8 (b).

As the line ∂η1(τa)/∂A2 = 0 is nearly linear at different kd, the line can be fixed if positions of points (−1,ya) and (A2 → 0,yb) are known. These points are labeled in Fig. 3. The line can be expressed as

It should be noted that the values g, η0, and L have no influence on the construction of ∂η1(τa)/∂A2 = 0. As they are eliminated when the value of τa is substituted into ∂η1(τa)/∂A2 = 0. In this way, coefficients a and b in Eq. (1) are only related with kd. Furthermore, the points (−1,ya) and (A2 → 0,yb) are only related with kd. When A2 is fixed, the variation rules between ya or yb with kd could be depicted by using the formula ∂η1(τa)/∂A2 = 0. Their relationship is shown in Fig. 4.

Fig. 4. The relationship between A1 and kd under the condition that A2 = −1 and A2 → 0 when ∂η1(τa)/∂A2 = 0

It is shown that A1 and kd are nearly exponential relationship when A2 = −1. When A2 → 0, A1 almost has no relationship with kd. Calculating by the data fitting, ya and yb can be defined as 0.289 e−3.252kd +0.678 and 0.396, respectively. Substituting the two points into Eq. (1), one can obtain

In this way, a variation rule of η1(τa) with A2 under the condition with an arbitrary concerned A1 and kd can be defined as

Equations (3) present a criteria to judge the correlation between η1 and A2 in the weakly nonlinear regime. It is shown that when the normalized thickness kd is fixed, the position of the line ∂η1(τa)/∂A2 = 0 will be defined. The study in Ref. [23] has shown that, when A1 = 1 and A2 = −1, the exist of the upper interface will accelerate the RTI perturbation growth in the weakly nonlinear regime on the lower interface. Using Eqs. (3), one can also make the same conclusion. From Eqs. (3), it is shown that η1 and A2 are negatively related when A1 = 1. In this way, comparing with the condition A2 → 0, the RTI perturbation on the lower interface with the condition A2 = 1 will be accelerated in the weakly nonlinear regime. Besides the extreme situation studied by Wang et al.,[23] the results presented in this study can be applied to more situations. Equations (3) are maybe useful to give some guidance to the domination of η1 in the weakly nonlinear RTI regime.

5. Underlying physics mechanism of η1 to A2

Following, the underlying physics mechanism is trying to be described. Firstly, the condition that A1 = 0, which means there is only an upper interface in the system, is concerned. According to the classical linear RTI growth model , the perturbation on the upper interface is absolutely stable when A2 < 0. Therefore, there are two kinds of movement trends for the perturbation on the upper interface. As shown in Fig. 5, the initial perturbation will make a downward movement following with an upward movement repeatedly. Besides, the smaller A2 is, the stronger the trend is. For convenience, this movement trend is called the self-excited response of the upper interface.

Fig. 5. Two movement trends of the self-excited response on the upper interface.

For the RTI with double interface studied in this study, the perturbation seed of the upper interface is from the perturbation evolution on the lower interface. Conversely, the perturbation evolution on the upper interface will impact the perturbation evolution on the lower interface. The downward movement trend of the self-excited response on the upper interface will accelerate the perturbation evolution on the lower interface. The upward movement trend of the self-excited response on the upper interface will suppress the perturbation evolution on the lower interface. It is clear that the value of A1 has influences on the intensity of the perturbation seed of the upper interface, while the value of kd impacts on the generation time of the perturbation seed on the upper interface. Firstly, the rationality of impacts of A1 is concerned. For a small A1, the perturbation intensity generated on the upper interface is small. When the RTI on the lower interface reaches the weakly nonlinear RTI regime, the self-excited response of the upper interface should present the upward movement trend. This trend will suppress the perturbation evolution on the lower interface. As the smaller A2 is, the stronger the trend is. A2 will be positively related to η1. However, for a large A1, the perturbation intensity generated on the upper interface is larger than the situation with a small A1. When the RTI on the lower interface reaches the weakly nonlinear RTI regime, the self-excited response of the upper interface should present the downward movement trend. This trend will accelerate the perturbation evolution on the lower interface. In this way, A2 will be negatively related to η1 when A1 is large. In terms of kd we concerned, one can see that it influences the gradient of the line ∂η1(τa)/∂A2 = 0. The larger the kd is, the smaller the gradient of this line will be. Following, a situation just above the line ∂η1(τa)/∂A2 = 0 in Fig. 3 is focussed on. When kd is shortened, the correlation between η1 and A2 in this situation will turn from negative relationship to positive relationship. This means that the self-excited response of the upper interface turns from the downward movement trend to the upward movement trend. This is reasonable because the generation time of the perturbation seed on the upper interface will be earlier when kd is shortened. This means that the perturbation evolution time on the upper interface is longer when the RTI on the lower interface reaches the weakly nonlinear regime. Therefore, the self-excited response of the upper interface has already presented the upward movement trend for a smaller kd. In this way, the gradient of ∂η1(τa)/∂A2 = 0 will be decreased when kd is shortened.

6. Conclusion

The interface coupling effects of the weakly nonlinear RTI with double interfaces are studied by taking into consideration of the two situations: the situation that the Atwood number on the lower interface is from 0 to 1 and the Atwood number on the upper interface is from −1 to 0. This is an expending research of Ref. [23]. It is shown that due to the self-excited response of the upper interface, the amplitude of perturbation growth on the lower interface is positively related with the Atwood number on the upper interface when the Atwood number on the lower interface is small. However, it is negatively related when the Atwood number on the lower interface is large.

Reference
[1] Rayleigh L 1883 Proc. London Math. Soc. 14 170
[2] Taylor G 1950 Proc. R. Soc. London 201 192
[3] Wang L F et al. 2017 Sci. China-Phys. Mech. Astron. 60 055201
[4] Gamezo V N et al. 2003 Science 299 77
[5] Tao Y S et al. 2012 Acta Phys. Sin. 61 075207 in Chinese
[6] Zhang J et al. 2018 Phys. Plasmas 25 022701
[7] Youngs D L 1984 Physica D: Nonlinear Phenomena 12 32
[8] Xu Z et al. 2009 Chin. Phys. Lett. 26 084703
[9] Betti R Hurricane O A 2016 Nat. Phys. 12 435
[10] Amendt P et al. 2010 Nucl. Fusion 50 105006
[11] Canaud B et al. 2011 Nucl. Fusion 51 062001
[12] Bodner S 1974 Phys. Rev. Lett. 33 761
[13] Takabe H 1985 Phys. Fluids 28 367
[14] Wang L F et al. 2012 Phys. Plasmas 19 012706
[15] Wang L F et al. 2012 Phys. Plasmas 19 100701
[16] Wang L F et al. 2016 Phys. Plasmas 23 052713
[17] Mikaelian K O 1982 Phys. Rev. 26 2140
[18] Mikaelian K O 1983 Phys. Rev. 28 1637
[19] Mikaelian K O 1990 Phys. Rev. 42 4944
[20] Mikaelian K O 1995 Phys. Fluids 7 888
[21] Mikaelian K O 1990 Phys. Rev. 42 3400
[22] Mikaelian K O 2005 Phys. Fluids 17 094105
[23] Wang L F et al. 2014 Phys. Plasmas 21 122710
[24] Guo H Y et al. 2017 Chin. Phys. Lett. 34 045201
[25] Guo H Y et al. 2017 Chin. Phys. 26 125202
[26] Einfeldt B 1988 SIAM Journal on Numerical Analysis 25 294
[27] Roe P L 1986 Ann. Rev. Fluid Mech. 18 337