⇦Chinese Physics Letters, 2017, Vol. 34, No. 4, Article code 045201 Weakly Nonlinear Rayleigh–Taylor Instability in Incompressible Fluids with Surface Tension * Hong-Yu Guo(郭宏宇)1,2, Li-Feng Wang(王立锋)2,3, Wen-Hua Ye(叶文华)2,3**, Jun-Feng Wu(吴俊峰)2, Wei-Yan Zhang(张维岩)2 Affiliations 1Graduate School, China Academy of Engineering Physics, Beijing 100088 2Institute of Applied Physics and Computational Mathematics, Beijing 100094 3HEDPS, Center for Applied Physics and Technology, Peking University, Beijing 100871 Received 12 November 2016 *Supported by the National Natural Science Foundation of China under Grant Nos 11275031, 11475034, 11575033 and 11274026, and the National Basic Research Program of China under Grant No 2013CB834100.
**Corresponding author. Email: ye_wenhua@iapcm.ac.cn Citation Text: Guo H Y, Wang L F, Ye W H, Wu J F and Zhang W Y 2017 Chin. Phys. Lett. 34 045201 Abstract A weakly nonlinear model is established for incompressible Rayleigh–Taylor instability with surface tension. The temporal evolution of a perturbed interface is explored analytically via the third-order solution. The dependence of the first three harmonics on the surface tension is discussed. The amplitudes of bubble and spike are greatly affected by surface tension. The saturation amplitude of the fundamental mode versus the Atwood number $A$ is investigated with surface tension into consideration. The saturation amplitude decreases with increasing $A$. Surface tension exhibits a stabilizing phenomenon. It is shown that the asymmetrical development of the perturbed interface occurs much later for large surface tension effect. DOI:10.1088/0256-307X/34/4/045201 PACS:52.57.Fg, 47.20.Ma, 52.35.Py © 2017 Chinese Physics Society Article Text When a lighter fluid accelerates a heavier one, the perturbed interface between the two fluids is subject to the Rayleigh–Taylor instability (RTI).[1,2] RTI plays an important role in several areas such as successful inertial confinement fusion (ICF)[3] and supernova explosions.[4,5] Many issues important to ICF, such as density–gradient stabilization effect, are investigated in laboratory 'water tank' experiments[16] where the surface tension effect is appreciable and the effect of surface tension should be considered. The classical RTI (CRTI) growth case [1,2] refers to two semi-infinite fluids with constant density $\rho^{\rm a}$ and $\rho^{\rm b}$ $(\rho^{\rm a}>\rho^{\rm b})$ undergoing a constant acceleration $g$. Initially a single-mode cosinusoidal perturbation with wavelength $\lambda$ on the interface will grow exponentially in time, $\eta(t)=\eta_0 e^{\gamma t}$, with the linear growth rate $\gamma=\sqrt{Akg}$, where $k=2 \pi/\lambda$ is the wave number, $\eta_0$ is the amplitude of perturbation at the beginning, and $A=(\rho^{\rm a}-\rho^{\rm b})/(\rho^{\rm a}+\rho^{\rm b})$ is the Atwood number. Based on the CRTI model, many studies have been carried out theoretically [6-9] and numerically [10-12] that take different effects into consideration. Bellman et al.[13] first extended the CRTI model to include surface tension $\sigma_{\rm s}$. Chandrasekhar[14] further investigated the effect of viscosity on the shear flow for the perturbation growth. Mikaelian[15] considered RTI in arbitrary stratified density profiles including surface tension. All these linear theories indicate that the significant effect of surface tension on the RTI is the existence of a cutoff wave number $k_{\rm c}$,[13-15] i.e., beyond which perturbations on the interface are stable, $$\begin{align} k_{\rm c}=[(\rho^{\rm a}-\rho^{\rm b})g /\sigma_{\rm s}]^{1/2}=2\pi/\lambda_{\rm c}.~~ \tag {1} \end{align} $$ As the amplitude of perturbation is comparable with the perturbation wavelength in the later time, the nonlinear effect becomes significant and the symmetry of the interface between two fluids is lost. The wave troughs tend to steepen into 'spikes' and wave crests tend to broaden to form 'bubbles'. There exists a weakly nonlinear (WN) regime of perturbation growth before a strong nonlinear regime. The WN growth of perturbation has been extensively investigated through the method of successive approximations.[16-20] In this study, the transition from linear to nonlinear RTI growth with surface tension is studied analytically through a third-order WN solution. We consider two inviscid, irrotational and incompressible fluids in a gravitational field $-g {\boldsymbol e}_y$. A Cartesian coordinate system is introduced, where $x$ and $y$ are along and normal to the unperturbed interface ($y=0$), respectively. In the following analysis, the physical quantities of the fluid above the fluid interface ($y>0$) shall be denoted by the superscript a and fluid below the interface ($y < 0$) by b. The perturbed interface is expressed as $y=0+\eta(x,t)$, where $\eta(x,y)$ is a small perturbation displacement. The governing equations for two incompressible fluids in the presence of surface tension are $$\begin{alignat}{1} &\nabla ^2\phi^{\rm a}=\nabla ^2\phi^{\rm b}=0,~{\rm in~two~fluids},~~ \tag {2} \end{alignat} $$ $$\begin{alignat}{1} &\partial_t\eta +\partial_x\eta \cdot \partial_x\phi^{\rm a}-\partial_y\phi^{\rm a}=0,~{\rm at}~ y=\eta(x,t),~~ \tag {3} \end{alignat} $$ $$\begin{alignat}{1} &\partial_t\eta +\partial_x\eta \cdot \partial_x\phi^{\rm b}-\partial_y\phi^{\rm b}=0,~{\rm at}~ y=\eta(x,t),~~ \tag {4} \end{alignat} $$ $$\begin{alignat}{1} &\rho^{\rm a}\Big[\partial_t\phi^{\rm a}+\frac{1}{2}(\partial_x \phi^{\rm a}+\partial_y \phi^{\rm a})^2+g \eta \Big]\\ &-\rho^{\rm b}\Big[\partial_t\phi^{\rm b}+\frac{1}{2}(\partial_x \phi^{\rm b}+\partial_y \phi^{\rm b})^2+g \eta\Big]\\ =\,&-\sigma_{\rm s} \cdot (1+\partial _x\eta )^{-3/2}\partial_{xx}\eta,~{\rm at}~ y=\eta(x,t),~~ \tag {5} \end{alignat} $$ where $\sigma_{\rm s}$ is the surface tension, $\rho^i$ and $\phi^i$ are respectively density and velocity potentials of the two fluids with $i$ denoting a or b. Equations (3) and (4) are kinematic boundary conditions that the normal velocity is continuous. Equation (5) represents that the pressure is continuous across the fluid interface. Considering an initially single-mode cosinusoidal perturbation $\eta(x,0)=\eta_0\cos(x)$ with $\eta_0$ being the amplitude of perturbation, the interface perturbation $\eta(x,t)$ and velocity potential $\phi^i(x,y,t)$ can be expanded into a power series in a formal parameter $\varepsilon$, $$\begin{align} \eta(x,t)=\,&\sum_{n=1}^{\infty}\varepsilon^{n}\sum_{m=0}^{\lfloor\frac{n}{2}\rfloor} \eta_{n,n-2m}(t)\cos[(n-2m)k x],~~ \tag {6} \end{align} $$ $$\begin{align} \phi^{\rm a}(x,y,t)=\,&\sum_{n=1}^{\infty}\varepsilon^{n}\sum_{m=0} ^{\lfloor\frac{n}{2}\rfloor} [\phi_{n,n-2m}^{\rm a} (t)e^{-k y}\\ &\times\cos[(n-2m)k x]],~~ \tag {7} \end{align} $$ $$\begin{align} \phi^{\rm b}(x,y,t)=\,&\sum_{n=1}^{\infty}\varepsilon^{n}\sum_{m=0} ^{\lfloor\frac{n}{2}\rfloor} [\phi_{n,n-2m}^{\rm b} (t)e^{k y}\\ &\times\cos[(n-2m)k x]],~~ \tag {8} \end{align} $$ where $\eta_{n,n-2m}$ ($\phi^i_{n,n-2m}$) is the amplitude of interface displacement (velocity potential) of the $(n-2m)$th Fourier mode in $n$th-order of $\varepsilon$ ($n=1,2,3,\ldots$; $m=0,\ldots,\lfloor n/2\rfloor$). Gauss's symbol $\lfloor n/2\rfloor$ indicates that the maximum integer is less than or equal to $n/2$. The velocity potential $\phi^i(x,y,t)$ has satisfied the Laplace Eq. (2) and condition $\nabla\phi^{\rm a}|_{y\rightarrow\infty}=0$ and $\nabla\phi^{\rm b}|_{y\rightarrow-\infty}=0$. Substituting Eqs. (6)-(8) into Eqs. (3)-(5), the resulting boundary conditions at the interface $y=\eta(x,t)$ can be obtained. Then based on the Taylor series, expanding every term of equations obtained before at the unperturbed interface $y=0$. The $n$th-order governing equations including $\varepsilon^n$ are derived. The $n$th-order equations are the second order ordinary differential equations that can be easily solved. At the first order $O(\varepsilon)$, the linear solution is $$\begin{align} \eta _{1,1}=\eta_0 \cosh (\gamma_1 t),~~ \tag {9} \end{align} $$ where $\gamma_1=\sqrt{Akg (1-\alpha^2)}$ is the linear growth rate of the perturbation,[14,15] with $A=(\rho^{\rm a}-\rho^{\rm b})/(\rho^{\rm a}+\rho^{\rm b})$ being the Atwood number, and $\alpha=k/k_{\rm c}$ being a normalized parameter which represents the intensity of surface tension, and $k_{\rm c}$ is the cutoff wave number defined in Eq. (1). Note that the perturbed interface is linearly stable when $\alpha>1 (k>k_{\rm c})$. At the second order $O(\varepsilon^2)$, the second order solution is obtained as $$\begin{alignat}{1} \eta _{2,2}=\,&\frac{\eta_0^2}{2}\frac{A k(1-\alpha^2)^2}{1-2\alpha^2-8\alpha^4}\Big\{\cosh (\gamma _2t)\\ &-\cosh ^2(\gamma _1t)+\frac{3 \alpha^2}{1-\alpha^2}\sinh^2(\gamma _1t)\Big\},~~ \tag {10} \end{alignat} $$ where $\gamma_2=\sqrt{2Akg(1-4\alpha^2)}$ is the intrinsic growth rate of the second harmonic. The last two terms of Eq. (10) come from the mode coupling effect of the fundamental mode. At the third order $O(\varepsilon^3)$, the third-order solutions are expressed as $$\begin{align} \eta _{3,3}=\,&\frac{\eta _0^3k^2}{2}\Big\{a_1\cosh ^3(\gamma _1t)+(a_3-a_2-a_1)\cosh (\gamma_3t)\\ &+a_2\cosh (\gamma _1t)- a_3\cosh (\gamma _1t) \cosh (\gamma _2t)\\ &+a_3\cdot \frac{9 \alpha^2}{1-\alpha^2}\frac{\gamma _1}{\gamma _2}\sinh (\gamma _1t) \sinh (\gamma _2t) \Big\},~~ \tag {11} \end{align} $$ $$\begin{align} \eta _{3,1}=\,&-\frac{\eta _0^3k^2}{4}\Big\{b_1\cosh ^3(\gamma _1t)+b_2\cdot\gamma _1 t\cdot \sinh (\gamma _1t)\\ &-(b_1-b_3)\cosh (\gamma _1t)-b_3\cosh (\gamma _1t) \cosh (\gamma _2t)\\ &-b_3\cdot \frac{6 \alpha^2}{1-\alpha^2}\frac{\gamma _1}{\gamma _2}\sinh (\gamma _1t)\sinh (\gamma _2t)\Big\},~~ \tag {12} \end{align} $$ where $\gamma_3=\sqrt{3Akg(1-9\alpha^3)}$ is the intrinsic growth rate of the third harmonic and the coefficients $a_i$ and $b_i$ ($i=1,2,3$) are $$\begin{alignat}{1} a_1=\,&A^2\frac{1-2\alpha^2+\alpha^4}{1+5\alpha^2+6\alpha^4} -\frac{1}{8}\frac{2+\alpha^2}{1+3\alpha^2},\\ a_2=\,&\frac{3\alpha^2 (1-\alpha^2)}{1-10 \alpha^2-39 \alpha^4}\Big(4A^2\frac{ 1-\alpha^2}{1+2 \alpha^2}-\frac{5}{8}\Big),~~ \tag {13} \end{alignat} $$ $$\begin{alignat}{1} a_3=\,&3 A^2 (1-\alpha^2)^3/(2-6 \alpha^2-93 \alpha^4-146 \alpha^6),\\ b_1=\,&3 A^2\frac{ 1-\alpha^2}{1+2 \alpha^2}+\frac{1}{4}\frac{4-7 \alpha^2}{1-\alpha^2},\\ b_2=\,&A^2\frac{ 1+7 \alpha^2-8 \alpha^4}{1-2 \alpha^2-8 \alpha^4}+\frac{1}{4}\frac{4-13 \alpha^2}{1-\alpha^2},\\ b_3=\,&4A^2(1-\alpha^2 )^3/(1-12\alpha^4-16\alpha^6),~~ \tag {14} \end{alignat} $$ where explicit $A^2$, $k^2$ and $\alpha$ dependencies of $\eta_{3,3}(\eta_{3,1})$ are exhibited clearly. As can be seen from Eqs. (9)-(14), if we set $\alpha=0$ (surface tension $\sigma_{\rm s}=0$), the WN results of CRTI[6,21] can be recovered. Comparing Eqs. (9)-(14) with Eqs. (17)-(23) in Ref. [21], additional mode coupling terms exist due to the modification of surface tension effect. Thus surface tension plays a remarkable role in the process of mode coupling, which will be discussed in the following. The influence of surface tension (parameter $\alpha$) on the development of perturbed interface for $A=1.0$ and 0.5 at normalized time $(Akg)^{1/2}t=5.1$ is plotted in Fig. 1. As can be seen in Fig. 1, when the amplitude of interface perturbation is comparable with the perturbation wavelength, the interface ($\alpha=0$) exhibits significant bubble-spike asymmetry. The interface asymmetry is obvious especially for a large Atwood number. With increasing $\alpha$, the interface between the two fluids tends to become stable. This means that surface tension has a definitely stabilizing effect that depends on the Atwood number $A$ on the interface.
cpl-34-4-045201-fig1.png
Fig. 1. The dependence of instability interface on surface tension (represented by parameter $\alpha$) for the Atwood number $A=1.0$ (a) and $A=0.5$ (b) at normalized time $(Akg)^{1/2} t=5.1$. The initial amplitude of perturbation is $\eta_0=0.001 \lambda$.
Temporal evolution of normalized amplitudes of bubble $\eta_{\rm b}/\lambda$ and spike $\eta_{\rm s}/\lambda$ for different parameters $\alpha$ at $A=1.0$ are shown in Fig. 2. With time proceeding, the amplitudes of spike and bubble grow rapidly. The amplitude of spike is greater than that of bubble for fixed time. As the surface tension $\alpha$ increases, the amplitudes of spike and bubble tend to grow equally.
cpl-34-4-045201-fig2.png
Fig. 2. Temporal evolution of the normalized amplitude of bubble $\eta_{\rm b}/\lambda$ (solid line) and spike $\eta_{\rm s}/\lambda$ (dashed line) for $\alpha=0$, 0.3 and 0.6 at the Atwood number $A=1.0$. The initial amplitude of perturbation is $\eta_0=0.001 \lambda$.
In Fig. 3, the temporal evolution of the normalized fundamental mode ($\eta_1=\eta_{1,1}+\eta_{3,1}$), the second harmonic ($\eta_2=\eta_{2,2}$) and the third harmonic ($\eta_3=\eta_{3,3}$) with different $\alpha$ at $A=1.0$ are illustrated explicitly. In Fig. 3(a), the amplitude of the fundamental mode grows quickly with time until the maximum amplitude is reached due to the negative feedback effect of the third order. Since the decreasing amplitude that occurs after the maximum amplitude is invalid in the theoretical model, the applicable scope of expansion theory is within the maximum amplitude of the fundamental mode. The larger the $\alpha$, the larger the maximum amplitude for fixed $A$. The second harmonic plays an important role in the formation of the bubble-spike structure. As can be seen in Fig. 3(b), the smaller the surface tension is, the larger the amplitude of the second harmonic is. This means that the RTI with smaller surface tension will lead to a broader bubble, but a narrower spike which is consistent with Fig. 1. The phase of the third harmonic is greatly affected by surface tension for fixed $A$ which can be seen in Fig. 3(c). The phase of the third harmonic is identical to that of the fundamental mode when $\alpha$ is small. With increasing $\alpha$, the amplitude of the third harmonic decreases to zero then grows rapidly with phase changed by $\pi$. When $\alpha$ approximates to $\sim0.56$, the amplitude of the third harmonic will not grow with time for $A=1.0$.
cpl-34-4-045201-fig3.png
Fig. 3. The temporal evolution of the normalized amplitude of perturbation $\eta_1/\lambda$ (a), $\eta_2/\lambda$ (b) and $\eta_3/\lambda$ (c) for $\alpha=0$, 0.3 and 0.6. The parameters are the same as those of Fig. 2.
As mentioned in Fig. 3(a), the fundamental mode grows gradually reaching a maximum amplitude. One can realize a saturation in the temporal evolution of fundamental mode because of the third-order negative feedback effect. If we define the transition regime from linear to nonlinear growth to happen when the amplitude of fundamental mode $\eta_1$ is reduced by $10\%$ in comparison with its linear growth $\eta_{1,1}$,[19-21] namely, $(\eta_{1,1}-\eta_1)/\eta_{1,1}=0.1$, then its corresponding time is called the nonlinear saturation time (NST) and the amplitude of fundamental mode is called the nonlinear saturation amplitude (NSA). If we drop the decaying terms and only keep the dominant terms (i.e., $e^{\gamma_1 t}, e^{2\gamma_1 t}, e^{3\gamma_1 t}$), the normalized NSA of the fundamental mode is obtained as $$\begin{alignat}{1} \!\!\!\!\!\!\!\!\frac{\eta_1^{\rm (s)}}{\lambda }=\frac{1}{\pi } \sqrt{\frac{8}{5}} \sqrt{\frac{1+\alpha^2-2 \alpha^4}{4+\alpha^2-14 \alpha^4+12 A^2(1-\alpha^2)^2}}.~~ \tag {15} \end{alignat} $$ When the surface tension effect is ignored ($\alpha=0$), Eq. (15) reduces to $$\begin{align} \frac{\eta_1 ^{\rm (s)}}{\lambda }=\frac{1}{\pi}\sqrt{\frac{2}{5}}\frac{1}{\sqrt{1+3A^2}},~~ \tag {16} \end{align} $$ where the result in Ref. [6] is recovered. The dependence of normalized NSA $\eta_1^{\rm (s)}/\lambda$ on the Atwood number $A$ for different $\alpha$ is displayed in Fig. 4. It can be seen that $\eta_1^{\rm (s)}/\lambda$ decreases with increasing the Atwood number $A$ for different $\alpha$. For fixed $A$, the normalized NSA increases with $\alpha$. For $\alpha=0.2$, 0.4 and 0.6 at $A=1.0$, $\eta_1^{\rm (s)}$ is respectively approximated to $\sim0.106\lambda$, $0.121\lambda$ and $0.154\lambda$, which are all larger than the classical NSA of $0.1\lambda$. Noting that the larger NSA means the later nonlinear growth. That is to say, the surface tension effect will mitigate the growth of perturbation.
cpl-34-4-045201-fig4.png
Fig. 4. The dependence of normalized nonlinear saturation amplitude $\eta^{\rm (s)}/\lambda$ on the Atwood number $A$ for different surface tension effects $\alpha=0.2$, 0.4 and 0.6. The classical NSA ($\alpha=0$) is also plotted for comparison.
In conclusion, a third-order WN model for the RTI growth considering surface tension effect has been proposed. The nonlinear solutions up to the third-order correction for initially single-mode interface perturbation are obtained analytically. Our WN results can depict the transition from linear regime to early stage of nonlinear regime where the interface asymmetry is very obvious due to the formation of bubble–spike structure. The effects of surface tension on the interface shape as well as the amplitude of bubble ($\eta_{\rm b}/\lambda$) and spike ($\eta_{\rm s}/\lambda$) are discussed in some detail. The higher the surface tension is, the more stable the perturbed interface is for fixed $A$. The larger $\alpha$ will tend to make the amplitudes of bubble and spike grow equally. The temporal evolution of the first three harmonics is illustrated explicitly. The amplitudes of the first three harmonics are greatly dependent on the surface tension. The phase of the third harmonic is also affected by surface tension. The NSA $\eta_1^{\rm (s)}/\lambda$ of the fundamental mode is strongly dependent on the Atwood number $A$ and surface tension $\alpha$. The $\eta_1^{\rm (s)}/\lambda$ decreases with increasing $A$. For a fixed Atwood number $A$, the normalized NSA increases with $\alpha$. References The Instability of Liquid Surfaces when Accelerated in a Direction Perpendicular to their Planes. ILaser Compression of Matter to Super-High Densities: Thermonuclear (CTR) ApplicationsThermonuclear Supernovae: Simulations of the Deflagration Stage and Their ImplicationsExperimental astrophysics with high power lasers and Z pinchesTwo-Dimensional Rayleigh–Taylor Instability in Incompressible Fluids at Arbitrary Atwood NumbersOn the Second Harmonic Generation through Bell—Plesset Effects in Cylindrical GeometryMagnetic field gradient effects on Rayleigh-Taylor instability with continuous magnetic field and density profilesRayleigh–Taylor instability of steady ablation fronts: The discontinuity model revisitedStabilization of ablative Rayleigh-Taylor instability due to change of the Atwood numberNumerical Simulation of Anisotropic Preheating Ablative Rayleigh–Taylor InstabilityFormation of jet-like spikes from the ablative Rayleigh-Taylor instabilityEffects of surface tension and viscosity on Taylor instabilityRayleigh-Taylor and Richtmyer-Meshkov instabilities in multilayer fluids with surface tensionThree-dimensional Rayleigh-Taylor instability Part 1. Weakly nonlinear theoryNonlinear Perturbation Theory of the Incompressible Richtmyer-Meshkov InstabilityWeakly nonlinear Bell-Plesset effects for a uniformly converging cylinderNonlinear saturation amplitudes in classical Rayleigh-Taylor instability at arbitrary Atwood numbersWeakly nonlinear Rayleigh-Taylor instability of a finite-thickness fluid layerCoupling between interface and velocity perturbations in the weakly nonlinear Rayleigh-Taylor instability
[1]Rayleigh L 1883 Proc. London Math. Soc. 14 170
[2] Taylor G 1950 Proc. R. Soc. London Ser. A 201 192
[3] Nuckolls J, Wood L, Thiessen A and Zimmerman G 1972 Nature 239 139
[4] Gamezo V N, Khokhlov A M, Oran E S, Chtchelkanova A Y and Rosenberg R O 2003 Science 299 77
[5] Remington B A, Drake R P and Ryutov D D 2006 Rev. Mod. Phys. 78 755
[6] Wang L F, Ye W H and Li Y J 2010 Chin. Phys. Lett. 27 025203
[7] Guo H Y et al 2014 Chin. Phys. Lett. 31 044702
[8] Yang B L, Wang L F, Ye W H and Xue C 2011 Phys. Plasmas 18 072111
[9] Piriz A R et al 1997 Phys. Plasmas 4 1117
[10] Ye W H, Zhang W Y and He X T 2002 Phys. Rev. E 65 057401
[11] Wang L F, Ye W H and Li Y J 2010 Chin. Phys. Lett. 27 025202
[12] Wang L F et al 2012 Phys. Plasmas 19 100701
[13] Bellman R and Pennington R H 1954 Quart. Appl. Math. 12 151
[14]Chandrasekhar S 1961 Hydrodynamic and Hydromagnetic Stability (London: Oxford University) chap X
[15] Mikaelian K O 1990 Phys. Rev. A 42 7211
[16] Jacobs J W and Catton I 1988 J. Fluid Mech. 187 329
[17] Velikovich A L and Dimonte G 1996 Phys. Rev. Lett. 76 3112
[18] Wang L F, Wu J F, Guo H Y, Ye W H, Liu J, Zhang W Y and He X T 2015 Phys. Plasmas 22 082702
[19] Liu W H, Wang L F, Ye W H and He X T 2012 Phys. Plasmas 19 042705
[20] Wang L F, Guo H Y, Wu J F, Ye W H, Liu J, Zhang W Y and He X T 2014 Phys. Plasmas 21 122710
[21] Wang L F et al 2012 Phys. Plasmas 19 112706