⇦Chinese Physics Letters, 2020, Vol. 37, No. 7, Article code 075201 Interface Width Effect on the Weakly Nonlinear Rayleigh–Taylor Instability in Spherical Geometry Yun-Peng Yang (杨云鹏)1,2, Jing Zhang (张靖)3, Zhi-Yuan Li (李志远)3, Li-Feng Wang (王立锋)2,3, Jun-Feng Wu (吴俊峰)3, Wun-Hua Ye (叶文华)2,3*, and Xian-Tu He (贺贤土)2,3 Affiliations 1School of Physics, Peking University, Beijing 100871, China 2Center for Applied Physics and Technology, HEDPS, and College of Engineering, Peking University, Beijing 100871, China 3Institute of Applied Physics and Computational Mathematics, Beijing 100094, China Received 17 February 2020; accepted 15 May 2020; published online 21 June 2020 Supported by the National Natural Science Foundation of China (Grant Nos. 11575033, 11675026, and 11975053).
*Corresponding author. Email: ye_wenhua@iapcm.ac.cn Citation Text: Yang Y P, Zhang J, Li Z Y, Wang L F and Wu J F et al. 2020 Chin. Phys. Lett. 37 075201    Abstract Interface width effect on the spherical Rayleigh–Taylor instability in the weakly nonlinear regime is studied by numerical simulations. For Legendre perturbation mode $P_n$ with wave number $k_n$ and interface half-width $L$, the commonly adopted empirical linear growth rate formula $\gamma_n^{\rm em}(L)=\gamma_n/\sqrt{1+k_nL}$ is found to be sufficient in spherical geometry. At the weakly nonlinear stage, the interface width affects the mode coupling processes. The development of the mode $P_{2n}$ is substantially influenced by the interface width. Moreover, the nonlinear saturation amplitude increases with the interface width. DOI:10.1088/0256-307X/37/7/075201 PACS:52.57.Fg, 47.20.Ma, 52.35.Py © 2020 Chinese Physics Society Article Text The Rayleigh–Taylor instability (RTI)[1,2] plays a key role in many astrophysical phenomena and inertial confinement fusion (ICF)[3–7]. When a heavy fluid is supported or accelerated by a light fluid, the interface separating the two fluids is unstable to the seeded perturbations. In the linear regime, a small perturbation on the interface grows exponentially. In planar cases, the growth rate is given by $\gamma=\sqrt{A_{\rm T}kg}$, where $A_{\rm T}$ is the Atwood number, $k$ is the wave number, and $g$ is the acceleration. When the amplitude of the perturbation is comparable to its wavelength, the RTI grows into a weakly nonlinear (WN) stage. A WN model based on the parameter expansion method is used to describe the WN RTI growth[8,9]. Harmonics are generated through mode coupling processes among lower-order modes. The growth of the fundamental mode is slowed down by the third-order feedback. As the perturbation continues to grow, the RTI eventually enters the fully nonlinear stage. In ICF and astrophysics, many processes happen in spherical geometry. The spherical RTI is more complicated than its planar counterpart. Besides the geometry, the RTI is also affected by many real physical conditions, such as compressibility, ablation,[10,11] and thin shell.[12,13] Particularly, the density profile is usually continuous. In other words, there is a density transition layer at the interface between the two fluids. Therefore, the interface width effect, i.e., the density gradient effect, is non-negligible. In ICF, the fusion target consists of multiple shells. The interface width normally exists at the shell interfaces, with a magnitude of 10 $\mathrm{µ m}$. Wang studied the interface width effect on planar RTI.[14] It is found that the interface width affects the linear growth and the mode coupling processes remarkably. In this letter, we present a simulation investigation on the interface width effect in 2D spherical geometry. The math framework and the numerical setup are introduced, and then the results are shown and discussed. For a perturbed spherical interface separating two irrotational and incompressible ideal fluids, the perturbation can be decomposed into Legendre polynomials in 2D cases, or into spherical harmonics in 3D cases. In 2D cases, the interface is initially located at $r(\theta,t=0)=R_{0}+\eta(\theta,t=0)$, where $R_{0}$ is the unperturbed interface radius, and $\eta$ is the perturbation. For a single-mode perturbation, at the time $t=0$, $\eta$ contains a single Legendre mode $\eta(\theta,t=0)=\eta_0 P_n(\cos\theta)$, where $\eta_0$ is the initial perturbation amplitude, and $n$ is the initial perturbation mode number. When the perturbation driven by the external gravity or acceleration grows, several modes emerge as a result of mode coupling, i.e., $\eta(\theta,t)=\sum_{l} \eta_l(t) P_l(\cos\theta)$. Zhang et al.[15] extended the third-order WN model to spherical geometry and studied the mode coupling processes. In this model, the perturbation and the velocity potentials are expanded into the power series of a formal parameter $\epsilon$ for each mode $P_l(\cos\theta)$, $$\begin{align} & \eta(\theta, t) = \displaystyle\sum^{\infty}_{l=0}\displaystyle\sum^{\infty}_{j=1}\epsilon^j a_{l}^{(j)}(t) P_{l}(\cos\theta),~~ \tag {1} \end{align} $$ $$\begin{align} & \phi^{\mathrm{ex}}(r,\theta,t) = \displaystyle\sum^{\infty}_{l=0}\displaystyle\sum^{\infty}_{j=1}\epsilon^j b^{\mathrm{ex}(j)}_{l}(t)\left(\frac{r}{R_{0}}\right)^{-(l+1)} P_{l}(\cos\theta), \;~~ \tag {2} \end{align} $$ $$\begin{align} & \phi^{\mathrm{in}}(r,\theta,t) = \displaystyle\sum^{\infty}_{l=0}\displaystyle\sum^{\infty}_{j=1}\epsilon^j b^{\mathrm{in}(j)}_{l}(t)\left(\frac{r}{R_{0}}\right)^{l} P_{l}(\cos\theta),~~ \tag {3} \end{align} $$ where $\phi$ is the velocity potential defined as ${\boldsymbol v}=\nabla \phi$, $a_{l}^{(j)}$ is the temporal amplitude of the perturbation at order-$\epsilon^j$ for each Legendre mode $l$, and $b_{l}^{(j)}$ is the temporal amplitude of the velocity potential at order-$\epsilon^j$ for each Legendre mode $l$. The superscript 'in' denotes the physical variables of the interior fluid within the interface, and the superscript 'ex' denotes the physical variables of the exterior fluid. The equation series of $\epsilon^j$ can be solved from the lower order to the higher order.
cpl-37-7-075201-fig1.png
Fig. 1. Density pseudocolor images for $P_n(\cos\theta)$ at $\gamma_n t=6.0$ with (a) $n=6,k_nL=0.2$, (b) $n=6,k_nL=1.2$, (c) $n=9,k_nL=0.2$, and (d) $n=9,k_nL=1.2$. The $y$ axis is the polar direction ($\theta=0$).
For the single-mode RTI, the solution at the first order of $\epsilon$ gives the linear growth of the fundamental mode $P_n$. No other $P_l$ emerges in the first-order solution. The growth rate of $P_n$ is $$\begin{alignat}{1} \gamma_{n}=\sqrt{A_{\rm T}^{\prime} g k_{n}} =\sqrt{\frac{n(n+1)\left(\rho^{\mathrm{ex}}-\rho^{\mathrm{in}}\right) g} {\left[(n+1) \rho^{\mathrm{in}}+n \rho^{\mathrm{ex}}\right] R_{0}}},~~ \tag {4} \end{alignat} $$ where $\rho$ is the density, ${\boldsymbol g}=-g {\boldsymbol e}_{r}$ is the external gravity, and $k_{n}=2\pi/\lambda_{n}=\sqrt{n(n+1)}/R_{0}$ is the wave number defined to be the square root of the Laplacian eigenvalue of Legendre polynomial $P_n$. The equivalent Atwood number $A_{\rm T}'=(\rho^{\mathrm{ex}}-\rho^{\mathrm{in}}) /(\sqrt{(n+1)/n}\rho^{\mathrm{in}}+\sqrt{n/(n+1)}\rho^{\mathrm{ex}})$ is used here, other than the traditional Atwood number $A_{\rm T}=(\rho^{\mathrm{ex}}-\rho^{\mathrm{in}}) /(\rho^{\mathrm{in}}+\rho^{\mathrm{ex}})$. The solutions within the third order[15] reveal the second-order and third-order mode coupling processes. The mode $P_0$ and the mode $P_{2n}$ are mainly generated through the second-order coupling. The feedback to the fundamental mode is mainly the results of the third-order coupling. The math framework introduced above does not contain interface width effects. We investigate the interface width effects in spherical geometry by numerical simulation (NS). The NS is run with the hydrodynamical code FLASH.[16] $A$ directionally split, piecewise-parabolic method solver is employed. For a perturbation mode $P_n$, the radius direction $r$ takes a range of $50\,\mathrm{µ m} < r < 950 \,\mathrm{µ m}$, and occupies $192n$ zones. For odd perturbation modes, the poloidal direction $\theta$ takes a range of $0 < \theta < \pi$, and occupies $512n$ zones. For even perturbation modes, which are symmetric about the equatorial plane, the poloidal direction takes a range of $0 < \theta < \pi/2$, and occupies $256n$ zones. Reflective boundary conditions are set for the poloidal boundary and the inner radial boundary. A hydrostatic boundary condition is set for the outer radial boundary. Initially, the interface is placed at $R_0 = 410\,\mathrm{µ m}$. The initial amplitude of the perturbation is set to be $\eta_0 = 0.002\lambda_n$. The gravity points toward the sphere center. The ideal gas equation of state is used, with the adiabatic index $\gamma_{h}$ set to $5.0$ to reduce compressibility. To induce the interface width effects, a density transition layer is employed at the interface. The initial density profile alongside the radial direction is of the form $\rho_0(r')=(\rho_{\mathrm{in}}+\rho_{\mathrm{ex}}) / 2 -[(\rho_{\mathrm{in}}-\rho_{\mathrm{ex}}) / 2] \tanh (r' / L)$, where $r'$ is the radial distance to the interface, and $L$ is the half width of the transition layer which satisfies $\min[\rho_0/(\mathrm{d}\rho_0/\mathrm{d}r')] = L/A_{\rm T}$. The interface is determined by density contour $\rho=(\rho^\mathrm{in}+\rho^\mathrm{ex})/2$ from the NS data. For even perturbation modes, we mirror the data about the equatorial plane to extend the interface range from $0 < \theta < \pi/2$ to $0 < \theta < \pi$. The amplitude $\eta_l$ for each mode $P_l$ is extracted from the interface through Legendre polynomial fitting. Figure 1 shows the density pseudocolor images at time $\gamma_n t=6.0$ for perturbation modes $P_6/P_9$ with Atwood number $0.5$.
cpl-37-7-075201-fig2.png
Fig. 2. Growth rate versus $k_nL$ for (a) $n=6,A_{\rm T}=0.2$, (b) $n=9,A_{\rm T}=0.2$, (c) $n=6,A_{\rm T}=0.5$, (d) $n=9,A_{\rm T}=0.5$, (e) $n=6,A_{\rm T}=0.8$, and (f) $n=9,A_{\rm T}=0.8$. The blue dots are the fitted growth rates $\gamma_n'$. The yellow lines are the empirical growth rates $\gamma_n^{\rm em}$ with an interface width effect. They are both normalized by the theoretical growth rates $\gamma_n$ for discontinuous density profile.
For perturbation mode $P_n$ with interface half width $k_nL$, the linear growth rate of the fundamental mode can be fitted from the NS data, which is denoted as $\gamma_n'$. Figure 2 plots $\gamma_n'/\gamma_n$ versus $k_nL$. As can be expected, when $k_nL$ is small, $\gamma_n'$ approaches $\gamma_n$, i.e., $\gamma_n'/\gamma_n$ approaches $1$. Similar to planar cases,[14] $\gamma_n'$ decreases with $k_nL$. This means that the interface width stabilizes the spherical RTI in the linear regime. In planar geometry with interface width effect taken into account, a widely used empirical formula of the growth rate is $$\begin{align} \begin{aligned} \gamma^{\rm em}=\gamma/\sqrt{1+kL}. \end{aligned}~~ \tag {5} \end{align} $$ In spherical geometry, a reasonable formula with interface width effect is $$\begin{align} \begin{aligned} \gamma_n^{\rm em}=\gamma_n/\sqrt{1+k_nL}. \end{aligned}~~ \tag {6} \end{align} $$ Here $\gamma_n^{\rm em}/\gamma_n$ is plotted in Fig. 2 as yellow lines. As is shown, for $P_6$ and $P_9$, when the Atwood number is small, $\gamma_n'$ agrees quite well with Eq. (6). When the Atwood number is large, Eq. (6) is still an appropriate estimation, but there may be an error up to $12\%$.
cpl-37-7-075201-fig3.png
Fig. 3. The zeroth mode amplitude versus the square of the fundamental mode amplitude with different $k_nL$ values for (a) $n=6$, $A_{\rm T}=0.2$, (b) $n=9$, $A_{\rm T}=0.2$, (c) $n=6$, $A_{\rm T}=0.5$, (d) $n=9$, $A_{\rm T}=0.5$, (e) $n=6$, $A_{\rm T}=0.8$, and (f) $n=9$, $A_{\rm T}=0.8$. The dotted lines are located at $\eta_0/\lambda _n =0$.
According to the third-order WN model,[15] the mode $P_0$ and the mode $P_{2n}$ mainly emerge from the second-order mode coupling processes. Figure 3 plots the $P_0$ amplitude $\eta_0$ versus the square of the fundamental mode amplitude with different $k_nL$ values. Figure 4 plots the $P_{2n}$ amplitude $\eta_{2n}$ versus the square of the fundamental mode amplitude with different $k_nL$ values. The lines in the two figures are approximately linear, indicating that the mode $P_0$ and the mode $P_{2n}$ are dominated by the second-order coupling, which is consistent with the WN model. From Fig. 3, one can see that $\eta_0$ is not sensitive to $k_nL$. In fact, the amplitude of the zeroth mode reflects the deviation of the average interface radius from $R_0$. It is a consequence of mass conservation in convergent geometry. Figure 4 shows that interface width affects $\eta_{2n}$ substantially. In subplots (c), (d), (e) and (f) of Fig. 4, $|\eta_{2n}|$ tends to decrease with $k_nL$, which indicates suppression of the mode coupling. However, in subplot (a), for $P_6$ with a small Atwood number, $\eta_{2n}$ changes its sign and increase with $k_nL$ when $k_nL$ is large. This result shows a dependent influence of the interface width effect on the mode coupling processes, which cannot be concluded as suppression or enhancement simply.
cpl-37-7-075201-fig4.png
Fig. 4. The $2n$th mode amplitude versus the square of the fundamental mode amplitude with different $k_nL$ values for (a) $n=6$, $A_{\rm T}=0.2$, (b) $n=9$, $A_{\rm T}=0.2$, (c) $n=6$, $A_{\rm T}=0.5$, (d) $n=9$, $A_{\rm T}=0.5$, (e) $n=6$, $A_{\rm T}=0.8$, and (f) $n=9$, $A_{\rm T}=0.8$. The dotted lines are located at $\eta_{2n}/\lambda _n =0$.
At the WN stage of RTI, high-order coupling processes will feed back to the fundamental mode. This feedback slows down the growth of the fundamental mode. In the scenario of the WN model, the feedback is dominated by third-order mode coupling. Based on the amplitude of the fundamental mode $\eta_n$ in the NS, its exponential growth in the linear stage can be fitted out as $\eta_{\rm L}=\eta_0'e^{\gamma_n' t}$, where $\eta_0'$ is the fitted initial perturbation amplitude. There is only a very small discrepancy between $\eta_0$ and $\eta_0'$, due to the initial conditions in the NS and the fitting error. The fundamental feedback can be calculated as $\eta_{\rm fd}=\eta_n-\eta_{\rm L}$. In most cases, $\eta_{\rm fd}$ is negative, and causes the fundamental mode to saturate. The nonlinear saturation time $t_{\rm s}$ is defined by $\eta_n(t_{\rm s})=0.9\eta_{\rm L}(t_{\rm s})$. The nonlinear saturation amplitude $\eta_{\rm sat}$ is defined as $\eta_{\rm sat}=\eta_{\rm L}(t_{\rm s})$. As $P_n(\cos 0)=1$, the amplitude of a Legendre mode equals the height of the bubble/spike at the pole ($\theta=0$). And the bubble/spike nearest to the equator is more similar to that in 2D planar cases. Thus for the convenience of comparison, we multiply $\eta_{\rm sat}$ by the height ratio of the initial bubble/spike nearest to the equator to the bubble/spike at the pole. That is, we use $\eta_{\rm s}=\alpha_n\eta_{\rm sat}$ instead of $\eta_{\rm sat}$ to represent the nonlinear saturation amplitude (NSA), where $\alpha_n$ is the height ratio of Legendre polynomial $P_n$. Figure 5 plots the NSA versus $k_nL$. The NSA tends to increase with $k_nL$, which means the high-order feedback to the fundamental mode is suppressed. However, for $P_6$ and $P_9$, when the Atwood number is small, this tendency is not evident and $\eta_{\rm s}$ can even decrease with $kL$. A possible explanation is that this is due to the reduction of the effective Atwood number. The Atwood number represents the density contrast between the two sides of the interface. The transition layer results in a reduced contrast, which can be regarded as a lowered 'effective Atwood number'. In planar cases, NSA always decreases with the Atwood number,[17] while in spherical geometry, NSA does not vary monotonically with the Atwood number. As shown in Fig. 6, for $P_6/P_9$ with $k_nL=0.1$, when the Atwood number is small, NSA increases with the Atwood number. Also in Fig. 5, it is seen that for a large $k_nL$ (i.e., small effective Atwood number), NSA increases with the Atwood number. Since one of the effects of the interface width is the reduction of the effective Atwood number, it is reasonable that $\eta_{\rm s}$ decreases with $k_nL$ for a small Atwood number.
cpl-37-7-075201-fig5.png
Fig. 5. Normalized NSA versus $k_nL$ for perturbation modes $P_6/P_9$ with Atwood numbers 0.2, 0.5, 0.8.
cpl-37-7-075201-fig6.png
Fig. 6. Normalized NSA versus $A_{\rm T}$ for perturbation mode $P_6/P_9$.
One may notice that the WN growth of $P_6$ is quite different from that of $P_9$. Reasons from two aspects can contribute to this difference. Firstly, $P_6$ has a lower mode number, which means that it is more affected by geometry effects. Secondly, the perturbation mode number $n$ is even for $P_6$, while it is odd for $P_9$. The mode coupling processes for even $n$ and odd $n$ are different.[15] For even $n$, the second-order coupling also generates the fundamental feedback, although it is the third-order coupling that dominates the feedback. In fact, for even $n$, all the modes generated by the second-order coupling can also arise from the third-order coupling. While for odd $n$, the second-order coupling and the third-order coupling generate different modes. In conclusion, we have studied interface width effect on the spherical Rayleigh–Taylor instability in the weakly nonlinear regime by numerical simulations. At the linear stage, the RTI growth rate decreases with the interface width. The commonly adopted empirical growth rate formula $\gamma^{\rm em}(L)=\gamma/\sqrt{1+kL}$, which takes the interface width into account and fits quite well in planar geometry, is found to be sufficient and slightly inaccurate in spherical cases. At the weakly nonlinear stage, the interface width affects the mode coupling processes. The development of the mode $P_0$, which is a result of mass conservation, is not sensitive to the interface width. While the development of the mode $P_{2n}$ is substantially influenced by the interface width. The nonlinear saturation amplitude, which is decided by high-order feedback, increases with the interface width. However, when the Atwood number is small, the relation between the NSA and the interface width is not evident. This can be explained due to the reduction of effective Atwood number. Our research takes a 2D spherical geometry, while real physical processes often take a 3D geometry. Future 3D spherical investigations can be conducted based on our 2D work. References The Form of Standing Waves on the Surface of Running WaterThe Instability of Liquid Surfaces when Accelerated in a Direction Perpendicular to their Planes. IThe physics basis for ignition using indirect-drive targets on the National Ignition FacilityA hybrid-drive nonisobaric-ignition scheme for inertial confinement fusionMain drive optimization of a high-foot pulse shape in inertial confinement fusion implosionsA scheme for reducing deceleration-phase Rayleigh–Taylor growth in inertial confinement fusion implosionsTheoretical and simulation research of hydrodynamic instabilities in inertial-confinement fusion implosionsThree-dimensional Rayleigh-Taylor instability Part 1. Weakly nonlinear theoryTwo-Dimensional Rayleigh?Taylor Instability in Incompressible Fluids at Arbitrary Atwood NumbersSelf-consistent growth rate of the Rayleigh–Taylor instability in an ablatively accelerating plasmaJet-Like Long Spike in Nonlinear Evolution of Ablative Rayleigh-Taylor InstabilityTime evolution of density perturbations in accelerating stratified fluidsPhase Effects of Long-Wavelength Rayleigh–Taylor Instability on the Thin ShellInterface width effect on the classical Rayleigh–Taylor instability in the weakly nonlinear regimeWeakly nonlinear incompressible Rayleigh-Taylor instability in spherical geometryFLASH: An Adaptive Mesh Hydrodynamics Code for Modeling Astrophysical Thermonuclear FlashesNonlinear saturation amplitude in the Rayleigh-Taylor instability at arbitrary Atwood numbers with continuous profiles
[1] Rayleigh L 1883 Proc. London Math. Soc. s1-15 69
[2] Taylor G 1950 Proc. R. Soc. London Ser. A 201 192
[3] Lindl J D et al. 2004 Phys. Plasmas 11 339
[4] He X T et al. 2016 Phys. Plasmas 23 082706
[5] Wang L F et al. 2016 Phys. Plasmas 23 122702
[6] Wang L F et al. 2016 Phys. Plasmas 23 052713
[7] Wang L F et al. 2017 Sci. Chin.-Phys. Mech. Astron. 60 055201
[8] Jacobs J W et al. 1988 J. Fluid Mech. 187 329
[9] Wang L F et al. 2010 Chin. Phys. Lett. 27 025203
[10] Takabe H et al. 1985 Phys. Fluids 28 3676
[11] Ye W H et al. 2010 Chin. Phys. Lett. 27 125203
[12] Mikaelian K O 1983 Phys. Rev. A 28 1637
[13] Li Z Y et al. 2020 Chin. Phys. Lett. 37 025201
[14] Wang L F et al. 2010 Phys. Plasmas 17 052305
[15] Zhang J et al. 2017 Phys. Plasmas 24 062703
[16] Fryxell B et al. 2000 Astrophys. J. Suppl. Ser. 131 273
[17] Wang L F et al. 2010 Europhys. Lett. 90 15001