Ultrasound Triggering of Rayleigh-Taylor Instability: Solution of Compressible Navier-Stokes Equation by a Non-Overlapping Parallel Compact Scheme

The Rayleigh-Taylor instability (RTI) [1, 2] is present in different branches of physics, from hydrodynamic applications [3, 4, 5, 6] to high energy physics of astrophysics [7, 8] and nuclear confinement fusion [9, 10]. Similarly, the baroclinic instability can arise during onset of RTI due to misalignment of density and pressure gradients, a major source of torque that contributes to the instability [11, 12]. In the present research, the compressible flow formulation of the Navier-Stokes equation (NSE) in [12] has been used. Such a formulation obviates the need of making the Boussinesq approximation, as needed for incompressible flow formulations. This has been critiqued by Mikaelian [13] highlighting the need for using compressible formulation. The baroclinic torque term is directly evident from the inviscid vorticity transport equation [11] as,

Buoyancy-induced acceleration due to gravity has been considered in [16, 17, 18, 19] for a tank with heavy fluid resting atop a lighter fluid, initially separated by a horizontal partition, removal of which causes RTI. A schematic of the physical domain along with initial conditions on either side of the interface is shown in Fig. 1. Such experiments allow detailed observation of the evolving mixing layer after the partition is removed at $t=0$, quite unlike those experiments where RTI is triggered by the acceleration of the chamber. In the present simulations, a similar configuration for RTI is studied, with large temperature difference (200K) between the top and bottom fluids causing the Atwood number ($=\frac{\rho_{u}-\rho_{l}}{\rho_{u}+\rho_{l}}$) to be $A_{t}=0.5$. Air is the working fluid that is treated as a perfect gas. In the present case, the box is considered to be thermally insulated, which enables it to be viewed as an universe, without any mass, momentum, and energy transfer with the outside. The ensuing RTI upon removal of the partition enables one to view the process as one of non-equilibrium thermodynamics, as reported for a two-dimensional (2D) flow in [20, 14, 21] and a three-dimensional (3D) flow with $A_{t}=0.2$, due to initial temperature differential of 149K. In Read’s experiments [22], two geometrical setups were considered with the purpose of providing results for 2D and 3D flows. The present research aims to consider only the 3D RTI setup with $A_{t}=0.5$, while in [12], corresponding 2D and 3D setups were computed by solving the 3D geometry with cross-sections of $(0.15m\times 0.025m)$ and $(0.15m\times 0.15m)$, respectively.

As implied in linear instability theory [2, 23, 24], RTI is caused at the interface of two fluids with dissimilar density due to three main parameters: (i) the imposed acceleration that causes the dense fluid (with density, $\rho_{u}$) to penetrate into the lighter fluid ($\rho_{l}$) as a spatio-temporal event; (ii) a non-dimensional density ratio given by the Atwood number ($A_{t}$) and (iii) a length scale given by a real wavenumber ($\alpha$). In fact, the dispersion relation is in general complex, but in many experimental [10], and numerical investigations [25, 26, 3] the length scale is considered real, which makes the linear disturbance growth temporal. In general, the transition process will be due to a spatio-temporal mechanism.

In practice, there is no need to impose any length scale as there are many 2D and 3D direct numerical simulations (DNS) where the initial interface is purely taken as planar [20, 14, 21, 27, 12] to replicate the experiments reported in [22, 28, 17] where the fluid interface is initially plane and the onset of RTI has been shown to originate from the edges and corners of the initial interface [12]. In contrast, many other simulations required the onset of the RTI to be at lower wavenumbers due to numerical reasons. The issue of numerical and experimental space-time resolution has been explained in [29], highlighting the improved experimental capabilities in recent times while discussing the parameter ranges of the Liepmann-Taylor scale also in the context of the design of capsules for inertial confinement fusion [30]. The fuel mixture is bombarded by a laser for the initiation process of fusion. RTI occurs at two instances: once during the initial implosion of the target and when the high-temperature gas mixture is decelerated in encountering the colder outer layer. The turbulent mixing during the RTI brings a colder outer layer of fuel to the core and, in the process, can suppress the ignition at the core [9]. The turbulent mixing following RTI is significantly influenced by the created vorticity field, not only during nuclear fusion but also during a volcanic eruption and in other geophysical events. Thus, the role of length scales in RTI is important, and it has been experimentally investigated in [25] by studying the distinction between experiments with and without imposed length scales. The authors in [12] have compared the length scale of the spikes (vortical structures originating from the heavier to lighter fluid) in the experiments reported in [22, 28, 16, 17, 25].

The above discussion indicates that the onset of RTI and associated length scales are still not unambiguously resolved. We have already noted that high accuracy compact schemes used for solving compressible NSE without the Stokes’ hypothesis enable one to track the onset of RTI to be from the edges and corners of the interface [20, 14, 21] for 2D RTI problem. The authors in [29] used a 2D multi-physics simulation package, HYDRA [31, 32] with a spatial resolution of the order of $3\mu m$ and time steps of the order of $100ps$. It was noted that for the nuclear fusion experiments using radiography, the ideal space-time resolution should be of the order of $2\mu m$ and $0ps$ to capture the evolving spikes! In contrast, the authors in [12] have used various grids with a representative one using $(320\times 640\times 320)$ in the computational box of dimension of $(150mm\times 150mm\times 150mm)$ i.e. about 65.536 million points with 160 cores. The results obtained matched the mixing layer and growth rate data of Read [22]. In the present research, we report results obtained using $(1280\times 2560\times 1280)$ i.e. 4.19 billion points for the same 3D setup of [22] using 19,200 cores. Thus, the spatial resolution has improved significantly, and the time step has been refined to $7.69\times 10^{-8}sec$ (with a non-dimensional time step of $6.25\times 10^{-7}$, as compared to a time step of $10^{-6}$ in [12]. More importantly, a novel parallel algorithm (non-overlapping high accuracy parallel (NOHAP) compact scheme) proposed in [33] is used. The parallelization scheme ensures the removal of any error originating at the sub-domain boundaries due to boundary closure schemes mandatory for compact schemes. Furthermore, the Schwartz domain decomposition method in [12], required overlap points at the subdomain boundaries following the parallelization method introduced in [34], to avoid errors originating due to parallelization. The NOHAP scheme used in the present research has been used earlier for different problems of fluid mechanics, as in [35, 36, 37].

The objective of the present research is to investigate the onset of RTI by the parallel algorithm by reducing computational errors down to the level of equivalent sequential computing of the same problem. The manuscript is formatted in the following manner. In the next section, we describe the formulation and various parameters used to solve the RTI problem. In section 3, a brief description of the parallel computing strategy and the associated numerical method is presented for the RTI, showing the linear scaling property of the parallelization. In section 4, the onset of RTI is shown with the appearance of acoustic waves and their propagation. The properties of such ultrasonic waves are explored. This is followed by section 5, where numerical solutions are shown to establish the relationship of the acoustic trigger at the onset with the vortical development that leads to the evolution of the mixing layer. The paper closes with a summary and conclusions in section 6.

The convective fluxes are $\hat{E}$, $\hat{F}$, $\hat{G}$, and given as

Similarly, the viscous flux vectors are denoted as $\hat{E}_{v}$, $\hat{F}_{v}$, $\hat{G}_{v}$, which are given as

In the above, the variables $\rho^{*}$, $p^{*}$, $u^{*}$, $v^{*}$, $w^{*}$, $T^{*}$ and $e_{t}^{*}$ are the dimensional density, pressure, Cartesian components of fluid velocity, absolute temperature and specific internal energy, respectively. The specific heat ratio ($\gamma$) is equal to 1.4 for air as assembly of diatomic molecules. The stress tensor are $\tau^{*}_{ij}$, for $i,j=1;to;3$, which are related to the strain rates given by the gradients of velocity as [12],

The specific heat capacity ($C_{v}$) and thermal conductivity ($\kappa$) are used as constants here. The heat conduction terms $q^{*}_{x}$, $q^{*}_{y}$ and $q^{*}_{z}$ are given as

We follow the procedure used in [12] to nondimensionalize the governing equations, with the variables related to appropriate scales as,

There are in-depth discussions on the concept of the bulk viscosity and on the validity of Stokes’ hypothesis in [12]. The authors have related the second coefficient of viscosity ($\lambda$) with the dynamic viscosity ($\mu$), by the expression $\lambda^{*}+\frac{2}{3}\mu^{*}=\mu_{b}$, where $\mu_{b}$ is the bulk viscosity of the fluid. Stokes [41] noted that in an analysis, with and without the bulk viscosity, if it produces different results, then it can be due to non-negligible values of divergence of velocity $\nabla\cdot\vec{V}$ contributing, along with $\mu_{b}$ being zero. For compressible flows, $\nabla\cdot\vec{V}$ is a significant non-negligible quantity. Past simulations of RTI [12, 14, 27] have shown the necessity to include the bulk viscosity to capture the early time behavior of RTI. Here, such an approach based on the acoustic dispersion and attenuation measurements in [15] for $\mu_{b}$ is adopted. Ash et al.’s experimental measurements provide values of $\mu_{b}$ for air as a function of temperature. Using a linear regression analysis of the experimental data, a relationship is drawn between $\mu_{b}$ and $T^{*}$ (in K) as [12],

Thus, the non-dimensional normal stress components are obtained as [12],

Finite difference compact schemes exhibit superior spectral resolution [42, 40] than the explicit schemes of same order, which motivated researchers to use them in many DNS and large eddy simulations [43, 20, 12, 35, 44, 45, 46, 47, 48]. The approximation of $m^{th}$ derivative of a function $f$ at $i^{th}$ point using a compact scheme is represented as,

The coefficients of compact schemes ($r$, $s$) derived for uniform grid spacing are independent of spatial location ($i$). The compact schemes can not be used at the first $max(n_{l},n_{r})$ points from the boundaries of computational domain and require boundary closure schemes as designed in [40, 49]. The compact schemes form a coupled system of equations, which can be represented as,

In many compact schemes [42, 40, 50], the $\left[A\right]$ and $\left[B\right]$ matrices are tri- and penta-diagonal matrices, respectively. The non-zero entries of $\left[A\right]$ and $\left[B\right]$ of OUCS3 scheme with $\eta=0$ [51, 40] for the $i^{th}$ node is given as,

The derivatives at a point depend on the derivatives at the neighboring points, which are obtained from the solution of Eq. (27). The tri-diagonal matrix algorithm [52] is used to solve the system of equations and get the derivatives at the grid points. The one-dimensional schematic of a non-overlapping domain decomposition used here is shown in Fig. 2. The $i^{th}$ and $(i+1)^{th}$ grid points represent interior grid points in the computational domain, which are now the last and the first grid points of the $p^{th}$ and $(p+1)^{th}$ processors, respectively in the parallel computing. The system of equations given in Eq. (27) have to be decoupled at $i^{th}$ and $(i+1)^{th}$ points so that, $p^{th}$ and $(p+1)^{th}$ processors can work in parallel, and solve the respective system of equations. An explicit scheme is to be used at the sub-domain boundaries in each processor to decouple the system of equations, which is the proposed sub-domain closure. The difference in the spectral resolution of the compact scheme used in the interior of the domain and the explicit scheme used at the sub-domain boundaries alters the numerical properties at a few grid points near the sub-domain boundaries [33]. This is the source of additional error in parallel computing of derivatives using compact schemes, a topic which has not been addressed except in [33].

There have been efforts in the past to reduce this error due to the parallelization. In [34], the domain is discretized with overlap points at the sub-domain boundaries, and one-sided implicit closure is used at the sub-domain boundaries. A large number of overlap points reduce the parallelization error, but which also affects the parallel performance of the simulation. There are also efforts in developing parallel compact schemes without requiring large overlap points. Such methods require one or two overlapping points. The use of an explicit central scheme at the extended sub-domain boundaries reduces bias at the sub-domain boundary point. However, instead of increased accuracy of the closure schemes, discontinuous resolution at the sub-domain boundaries has been noted as a source of error [33]. To improve the efficiency of large-scale simulations, one would like to decompose the domain with fewer overlap points without compromising the accuracy by filtering the solution [48, 47, 34]. Kim [45] used an optimized boundary closure scheme, which was extended in [46] by including halo nodes near sub-domain boundaries. Fang et al. [44], have proposed a parallel compact scheme without overlap points using a similar approach in [53], with an explicit scheme at the sub-domain boundary of the same order. However, these methods can not fully eliminate errors due to the parallelization of compact schemes. Moreover, the sub-domain closures present in [53, 46, 44] are specifically designed to work with a chosen compact scheme. In contrast, the NOHAP sub-domain closure [33] used in the present work does not require overlapping points at the sub-domain boundaries. Computationally this is equivalent to solving the sub-domain boundary point as equivalent to an interior point. Thus, this is equivalent to treating the full domain as part of a single domain as used in sequential computing, retaining the ability to control the error decided upon by machine precision. This is explained next.

The coefficients of the equivalent explicit OUCS3 scheme at an interior node are given in the Appendix. The magnitude of coefficients of $\left[C\right]$ matrix decay exponentially [33] with respect to diagonal. The semi-bandwidth $(n_{b})$ of a compact scheme for double-precision computations is the nodal location from the diagonal of $\left[C\right]$ matrix beyond which the coefficients are less than $10^{-16}$. The $n_{b}$ of OUCS3 scheme is found to be $48$ [33]. The solution of Eqs. (27) and  (29) with correct $n_{b}$ are identical up to machine precision level at which differences arise due to round-off error in computing. The solution via equivalent explicit scheme given in Eq. (29) requires $96$ floating-point operations. This is used only to obtain derivatives at the sub-domain boundaries and decouple the system of equations in Eq. (27). Then processors independently solve Eq. (27) to obtain derivatives from second to second-last points in each sub-domain. The NOHAP sub-domain closure is obtained from the compact scheme used in the interior domain. Thus, the spectral resolution is unchanged across the sub-domain boundaries. The NOHAP sub-domain closure eliminates the error due to parallelization completely up to machine precision level [33]. This methodology of sub-domain boundary closure can be used with any parallel schemes to remove the error at sub-domain boundaries for the chosen precision level. This is established with the help of the DNS of RTI next using the OUCS3 scheme.

The parallel performance of NOHAP closure is evaluated from the 3D RTI simulations using $1920$, $2064$, $2400$, $4800$, $9600$, and $19200$ processors. All the simulations are performed in the Cray XC40 machine with Intel Haswell $2.5~{}GHz$ CPUs that are installed at SERC, IISc, Bangalore, India. The computing node contains $24$ CPUs and $128~{}GB$ RAM. Floating-point computations are performed with double-precision accuracy. The simulations require $\approx 3.5~{}TB$ of RAM to solve the 3D compressible NSE in the chosen grid size with ($1280\times 2560\times 1280\approx 4.19$ billion points), and the grid size is kept the same for the tests. The code does not perform writing/reading data to/from the file system during the test. The time taken to perform time integration is shown as a function of the number of processors in Fig. 3. The dash-dotted line in the figure shows the linear scalability based on the results obtained with $1920$ and $2064$ processors. The 3D RTI simulations with NOHAP closure display linear scaling up to $19200$ processors. This is more than $99\%$ of parallel efficiency with respect to $1920$ processors result observed up to $4800$ processors, and $91.08\%$ of performance is retained while using $19200$ processors. The results demonstrate an excellent parallel efficiency and linear scaling of the solver with NOHAP closure. A careful design and implementation of NOHAP closure allow one to retain the high accuracy of compact schemes along with a good parallel performance, demonstrating the potential of NOHAP closure for even exascale computing.

The present results are computed using ($20\times 48\times 20$) = 19,200 processors and NSE solver takes $\approx 0.85s$ per time step. Onset of RTI is captured by solving NSE up to $t^{*}=0.7\sqrt{L/g}$, which requires $3,686,400$ core-hours of computation with Intel Haswell $2.5~{}GHz$ CPUs installed at SERC, IISc, Bangalore, India. The OUCS3 scheme [51, 40] with NOHAP sub-domain closure is used to discretize the convective fluxes, and the second-order explicit finite difference scheme is used to discretize the viscous flux and other gradients. The four-stage, fourth-order Runge-Kutta scheme [40] is used to perform time integration. The implicit 3D filter [54, 40] with 0.135 as filtering coefficient is applied at every $50$ timesteps is applied to eliminate high wavenumber oscillations in the simulation, and the results are discussed in this section.

At the initial unperturbed state, one notes that at the interface of the compartments, the density and pressure gradient act in the upward and downward directions, respectively. The impulsive removal of partition between the air masses propagates the information as pressure waves in the direction normal to the interface as shown in Fig. 4. The nondimensionalized instantaneous pressure field is subtracted from the initial hydrostatic pressure and plotted for $(z^{*}=0.5L)$-plane at indicated times in the figure. The pressure waves are reflected at the top and bottom walls, which travel towards the interface. The interface further partially reflects the traveling pressure waves, and a complex wave system is formed. The magnitude of the disturbance pressure is in the order of $10^{-6}$, which is many orders of magnitude smaller than the hydrostatic pressure. The error triggered by an improper parallelization of compact schemes is shown to be of a similar order of magnitude [33]. In contrast, the NOHAP closure removes parallelization errors at the sub-domain boundaries and accurately computes this multi-physics problem. The zoomed view in the bottom-right frame of the figure shows a visible vortical structure of the evolving RTI at the indicated time.

Figure 4 shows that the disturbance pressure field at the initial times is varying only in the direction normal to the interface. Thus, the variation of disturbance pressure at $x^{*}=z^{*}=0.5L$ is plotted in Fig. 5 at the indicated time instants to understand the evolution of disturbances. The discontinuity at the interface excites wavenumbers with an amplitude proportional to $1/\alpha$ [24]. The pressure waves’ genesis is from the interface, which travel towards the top and bottom walls. The high wavenumber components marked as ($HW$) in the frames diffuse quickly over time. The components identified as $A$ and $B$ in the figure propagate towards the top and bottom walls, respectively, and are reflected by the walls. The component $A$ is partially reflected by the interface at $t=0.0065$, and component $C$ is created, which propagates in the opposite direction, as shown in the last frame of the figure.

The corresponding Fourier transform of disturbance pressure field is shown in Fig. 6 for the indicated time instances. The peaks at low wavenumbers ($\alpha^{*}<600/L$) in Fig 6 correspond to the variation of disturbance pressure at the extreme $y$ locations [24] in Fig 5. The high wavenumber components ($\alpha^{*}>7600/L$) created by the discontinuity in the initial condition diffuse over time. The central wavenumber and corresponding magnitude of the pressure wave are marked with circles in Fig. 6. The strength and wavenumber of the pressure wave decay over time.

The disturbance pressure at $x^{*}=z^{*}=0.5L$ is plotted in the $(y,t)$-plane, in Fig 7. The peaks of the pressure waves at early times are connected with dash-dotted lines in the figure. The slope of the dash-dotted line provides the speed of propagation of the pressure waves in the fluid domain. The slope of the lines are equal to the local speed of sound, which is given as,

on the top and bottom side of the interface, respectively. The dispersion relation of the linear convection equation ($\omega_{0}=a\alpha$) [40] provides the range of frequency of acoustic waves as, $0.8~{}MHz$ - $3.73~{}MHz$. These ultrasonic sound waves serve as a deterministic mechanism that triggers the RTI of air masses separated by an unperturbed interface and ensures the reproducibility of the simulations.

In the 2D RTI simulation [14], it has been noted that the initial 1D pressure pulses become 2D complex wave systems, and one observes the visual signature of vortical structures that leads to the evolution of mixing at the interface. In the present 3D RTI simulation, it is important to ascertain how truly the onset of RTI takes place and its relation to initial ultrasounds reported so far in the above. It has been noted by Chandrasekhar [23] and in visualization experiments [38, 22] that the visible onset of RTI is due to vorticity evolution at the interface, and this has been attributed to the baroclinic instability. Thus, in the following section, one is interested in pinpointing the onset of RTI with the evolution of vorticity at the interface plane and whether the other terms of vorticity transport equation also contribute to it.

Video-1 (click caption/figure to play the video): The evolution of disturbance pressure at $z=0.5$ plane.

The evolution of temporal growth of vorticity in the interface plane is investigated using Eq. (1). This equation provides the clue of the temporal growth rate of vorticity that is typical of the present simulation based on the specific formulation. The various terms present on the right-hand side of Eq. (1) contribute to the growth rate of vorticity at the ($y^{*}=0.5L$)-plane, which are due to (i) vortex-stretching (typical of 3D formulation), (ii) contribution due to dilation and compression locally by the divergence of the velocity term (due to bulk viscosity contribution) and also (iii) due to misalignment of density gradient from the pressure gradient (baroclinic term). These are present, apart from viscous term contribution an increase of enstrophy [55]. The above three terms on the right-hand side of Eq. (1) are plotted for an early time ($100\Delta t$) in Fig. 8. These stretching, divergence and baroclinic terms are plotted on the top, middle, and bottom frames, respectively. While the left column shows the $x$-component (which is identical to the $z$-component), the right column shows the $y$ component of each term. The results show that the growth of RTI begins from the corner and edges of the interface, along with the temporal growth of vorticity triggered initially due to the acoustic pulses. The maximum amplitude of each term is shown in each frame, which clearly indicates that the baroclinic term contributes more than $97\%$ to the total growth of vorticity at the interface. The contribution by the bulk viscosity and vortex stretching terms contribute to the rest by almost the same order of magnitude. Thus, the acoustic pulses and their reflections (whose origin can be traced back to the ultrasound created at the onset of the experiment) from the walls cause the misalignment of the pressure and density gradient vectors, leading to baroclinic vorticity production at the interface. The convection created by the generated vorticity enhances the mixing of the heavy and lighter fluids as investigated in [22, 14, 12]. The temporal evolution of disturbance pressure at $(z^{*}=0.5L)$-plane is shown in Video-1. The ultrasonic acoustic waves predominately travel in the direction orthogonal to the initial interface and suffer multiple reflections from the walls and the interface. The frequency and wavenumber of the acoustic waves decay over time. The bottom-right frame in the animation shows the onset of visible structures of RTI due to background acoustic disturbance.

The onset mechanisms involved in the buoyancy-driven Rayleigh-Taylor instability (RTI) are numerically investigated in the present work by performing 3D DNS of a flow field consisting of air at two different temperatures (and hence densities), initially separated by a non-conducting partition. The onset of the numerical experiment involves the removal of this partition at $t=0$. The computations involved solving the 3D NSE for compressible flows by adopting a parallel algorithm which allows one to operate with computing error akin to that of an equivalent sequential computation. This NOHAP strategy for parallelization allowed one to capture very high-frequency acoustic signals in the RTI, which are the harbingers of the onset mechanism, by computing using $\approx$ 4.2 billion mesh points and a time-step of $7.68\times 10^{-8}$.

The parallel performance of the NOHAP scheme has been demonstrated by a series of 3D RTI simulations using 1920, 2064, 2400, 4800, 9600, and 19,200 processors. A $91.08\%$ retention of parallel performance has been demonstrated using 19,200 processors, suggesting a linear scaling of the present solver. Pressure waves reflected back and forth from the interface to the top/bottom walls were captured using the error-free parallelization tool of NOHAP. These pressure waves, in turn, form a complex wave system at the interface, which aids the formation of baroclinic vorticity. This vorticity is the seed for the onset of the RTI, as has been reported in prior works [23, 11, 12]. The high wavenumber components of these acoustic waves diffused quickly with time, and the typical range of frequencies associated with these are between 0.8 to 3.73 MHz. These ultrasonic waves serve as a deterministic trigger for the generation of baroclinic vorticity, and hence the onset of the RTI.

The onset and propagation of RTI are from the junction of the corner and edges of the 3D box simulated with the interface. This is followed by a spatio-temporal growth of the vorticity triggered at these junction points due to the acoustic waves. The further evolution of the instability from this onset to subsequent development of the mixing layer is aided by the convection of the generated vorticity. An analysis of the terms of the vorticity transport equation revealed that baroclinic term contributes to more than $97\%$ of vorticity generation during onset, confirming the commonly adopted viewpoint regarding RTI inception [23, 22].

The authors are grateful to National Supercomputing Mission, Government of India, Grant No. DST/NSM/R$\&$D_Exascale/2021/17 for sponsoring the project and providing full financial support for the problem solved here as part of an NSM initiative. The authors thank SERC, IISc, Bangalore, India, for providing the computing facility for the simulations.