A Relaxation Filtering Approach for Two-Dimensional Rayleigh–Taylor Instability-Induced Flows

In this paper, we investigate the performance of a relaxation filtering approach for the Euler turbulence using a central seven-point stencil reconstruction scheme. High-resolution numerical experiments are performed for both multi-mode and single-mode


Introduction
Rayleigh-Taylor instability (RTI) is an interfacial hydrodynamic instability that occurs at the interface separating two fluids of different densities in the presence of relative acceleration [1]. Understanding the behavior of the RTI-induced flows are of great importance because of the prevalence of such instability phenomena in many natural, industrial, and astrophysical systems with unstably stratified interfaces, such as coastal upwelling near the surface of the oceans [2], atmosphere and clouds [3], plasma physics such as magnetic or inertial confinement fusion implosions [4], the ignition of supernova [5,6], air bubble formation in the blood of deep sea divers [7], premixed combustion [8,9] and many more. In general, RTI phenomenon is one of the easiest hydrodynamics instabilities to observe, for example, if we invert a glass filled with water, the RTI occurs which makes the water falling [10,11]. Lord Rayleigh first described theoretically an instability that occurs in Cirrus cloud formation when a dense fluid is supported by a lighter one in a gravitational field [12]. Later Sir G. Taylor demonstrated the same instability experimentally for accelerated fluids [1], and honouring to their contributions, this instability is named after Lord Rayleigh and Sir G. Taylor which is the Rayleigh-Taylor instability. A detailed overview of the application of RTI phenomena with definitions, physical interpretations and terminologies can be found in [13][14][15][16][17]. Although RTI is a part of many diverse areas of scientific research, and there have been a substantial body of works conducted on this of a turbulent flow governed by Eulerian hyperbolic conservation laws, ILES methodology can be a good choice to consider which is proven to show a good performance on resolving turbulent flows with shock and discontinuities [45][46][47][48]. One of the popular ILES framework is to use an upwind scheme (e.g., Weighted Essentially Non-Oscillatory (WENO) scheme) along with a Riemann solver to incorporate the artificial through the use of numerical truncation errors [49]. The upwind-biased and nonlinearly weighted WENO schemes are widely used in resolving highly compressible turbulent flows because of their robustness to capture discontinuities in shock dominated flows and high order of accuracy in preserving turbulence features [50][51][52]. It should be mentioned that the development of improved WENO scheme is an active research field and still, there are a lot of works going on in this research direction [53][54][55][56]. Another candidate modeling approach is the explicit filtering approach using relaxation filtering which add dissipation on the truncated scales in LES through a low-pass filter [57][58][59][60]. In this approach, an additional low-pass spatial filter is used to estimate the effect of unresolved scales. The selection of relaxation filter also affects the solution field for explicit filtering approach and there are a significant number of literature available on the formulation of suitable and efficient LES filters [61][62][63]. In this work, we use the sixth-order symmetric central scheme with a 7-point stencil Simpson's filter (SF7) as a relaxation filter for RTI test case. We also implement ILES scheme combined with Roe and Rusanov Riemann solver to compare the results obtained by our relaxation filtering solver. The main purposes of this paper are to simulate RTI-induced flow (for both single and multi-mode perturbation) using a relaxation filtering approach to observe the resolution capability of this scheme, analyze the flow behavior by observing the density field contours, and compare the results obtained by relaxation filtering scheme and ILES scheme through kinetic energy (with and without density-weighted velocity) and power density spectra plots. The results show that the relaxation filtering scheme captures more scale in inertial subregion whereas the ILES scheme resolves more scales in high wavenumber regime. For both multi-and single-mode RTI, the kinetic energy spectra plots tend to follow the k −11/5 scaling law. On the other hand, the power density spectra plots are observed to be aligned to k −7/5 at high resolution. More rigorous derivations and mathematical analyses of scaling laws can be found elsewhere [64][65][66][67][68][69]. Also, the density contour plots for single-mode RTI reveal that the symmetry of the falling spike break at high resolution because of the formation of secondary instabilities from smaller scales resolved at high resolution. However, the lower resolution simulations resolve the symmetry for all schemes since the numerical dissipation surpasses the formation of the secondary instabilities. Also, it has been seen that the filter strength of the relaxation filter allows us to add more or less dissipation to the solver which eventually affects the flow behavior.
The rest of the paper is organized as follows. Section 2 gives a brief description of the governing equations. Section 3 illustrates the numerical methodology implemented in this study. In Section 4, we detail the results obtained by the numerical schemes used in our investigation along with the problem definitions of the RTI test problem. We demonstrate our findings through high-and coarse-resolution density field contour and density-weighted energy spectra plots. Section 5 gives the summary of our findings and conclusions.

Governing Equations
In our study, we consider two-dimensional Euler equations in their conservative dimensionless form as underlying governing equation for the Rayleigh-Taylor instability-induced flow evolution and can be expressed as: where F, G account for the inviscid flux contributions to the governing equation and S represents the gravitational term acting on the vertically downward direction (i.e., g = −1). The quantities included in q, F, G and S are: Here, ρ, p, e, u, and v are the density, pressure, total energy per unit mass, and the horizontal and vertical velocity components, respectively. The total enthalpy, h and pressure, p can be obtained by: where γ = 7/5 is chosen as the ratio of specific heats in our study. We refer the reader to [50,70] for details on the development of the eigensystem of the equations to devise hyperbolic conservation laws.

Numerical Methods
To develop the computational algorithm for our test problem governed by hyperbolic conservation laws, we formulate a finite volume framework by using different numerical strategies and schemes. In this section, we briefly introduce the numerical methods considered in the present study. We use the method of lines to cast our system of partial differential equations given in Equation (1) in the following form of ordinary differential equation through time: where q i,j is the cell-averaged vector of dependent variables, and £(q i,j ) represents the convective flux terms in the governing equation which can be expressed in the following discretized form: Here, F i±1/2,j are the cell face flux reconstructions in x-direction and G i,j±1/2 are the cell face flux reconstructions in y-direction. We use the optimal third-order accurate total variation diminishing Runge-Kutta (TVDRK3) scheme [71] for the time integration: where the time step, ∆t should be obtained by (satisfying the Courant-Friedrichs-Lewy (CFL) criterion): , η ∆y max(|v|, |v + a|, |v − a|) , where a is the speed of the sound that can be computed from the primitive flow variables (i.e., a = γp/ρ). In our current investigation, we use η = 0.5 for all the simulations (η ≤ 1 for numerical stability). For the cell face flux reconstructions, we have implemented the ILES and relaxation filtering modeling approaches on our test problem which will be discussed briefly in the subsequent sections.

ILES Approach
To develop our ILES framework, we first use the WENO interpolation scheme to reconstruct the left and right state of the cell boundaries. Later, we calculate the fluxes at cell edges from the reconstructed left and right states using a Riemann solver. The finite volume framework of a system of Euler conservation equations usually requires a Riemann solver to avoid the Riemann problem [72]. The damping characteristics of nonlinear WENO schemes acts as an implicit filter to prevent the energy accumulation near the grid cut-off [73,74]. In this work, we use the 5th order accurate WENO scheme followed by two widely used Riemann solver, Roe and Rusanov Riemann solver, to determine the flux at cell boundaries.

Weno Reconstruction
The WENO scheme is first introduced in [75] for problems with shocks and discontinuity to get an improvement over the essentially non-oscillatory (ENO) method [76,77]. In this work, we use an implementation of the WENO reconstruction using 7-point stencils (i.e., updating any quantity located at index i depends on the information coming from i − 3, i − 2, ..., i + 3) which can be written as: . (11) Here, q L i+1/2 and q R i−1/2 are the left state (positive) and right state (negative) fluxes, respectively, approximated at midpoints between cell nodes. The left (L) and right (R) states correspond to the possibility of advection from both directions. Since the procedures are similar in the y-direction, we shall present stencil expressions only in the x-direction for the rest of this document. w k are the nonlinear WENO weights of the kth stencil where k = 0, 1, ..., r and r is the number of stencils (r = 2 for the WENO5 scheme). The nonlinear weights are proposed by Jiang and Shu [78] in their classical WENO-JS scheme as: but the nonlinear weights defined by the WENO-JS scheme are found to be more dissipative than many low-dissipation linear schemes in both smooth region and regions around discontinuities or shock waves [79]. In our study, we have used an improved version of WENO approach proposed by [80], often referred to as WENO-Z scheme. One of the main reasons behind selecting WENO-Z can be less dissipative behavior than classical WENO-JS to capture shock waves. Also, there is a smaller loss in accuracy at critical points for improved nonlinear weights. The new nonlinear weights for the WENO-Z scheme are defined by: where β k and p are the smoothness indicator of the kth stencil and a positive integer, respectively. Here, = 1.0 × 10 −20 , a small constant preventing zero division, and p = 2 is set in the present study to get the optimal fifth-order accuracy at critical points. The expressions for β k in terms of cell values of q are given by: d k are the optimal weights for the linear high-order scheme which are given by:

Roe Riemann Solver
Based on the Godunov theorem [72], Roe developed an approximate Riemann solver, known as the Roe Riemann solver [81]. In our computational algorithm, we use the flux difference splitting (FDS) scheme of Roe [81] where the exact values of the fluxes at the interface can be computed in the x-direction by: Here, ∆F is the flux difference, calculated as: where Here, ∆ denotes the difference between right and left state fluxes for the variables ρ, p, u, v (e.g., ∆u = u R − u L ), and eigenvalues are defined as λ 1 = |ũ|, λ 2 = |ũ +ã| and λ 3 = |ũ −ã|, whereã is the speed of the sound at averaged state. In the equations, the tilde represents the density-weighted average, or the Roe average, between the left and right states. The Roe average values can be found by: where the left and right states of the un-averaged conserved variables are available from the WENO5 reconstruction described earlier. However, it is later realized that the stationary expansion shocks are not dissipated appropriately by this method. To fix the entropy in the expansion shocks, Harten proposed the following approach [82] replacing Roe averaged eigenvalues by: Here, = 2κã where κ is a small positive number, is set 0.1 in our computations. Similarly, in y-direction, λ 1 = |ṽ|, λ 2 = |ṽ +ã| and λ 3 = |ṽ −ã|. The interfacial fluxes in y-direction can be estimated by: where with

Rusanov Riemann Solver
Rusanov proposes a Riemann solver based on the information obtained from maximum local wave propagation speed [83], sometimes referred to as local Lax-Friedrichs flux [84,85]. The expression for Rusanov solver in the x-direction is as follows: where the right constructed state flux component, i+1/2,j ) and the characteristic speed, c i+1/2 =ã + |ũ|. The density-weighted average of the conserved variables can be calculated by Equation (18). Similarly, the expression for Rusanov solver in y-direction is: where c j+1/2 =ã + |ṽ|.

Central Scheme with Relaxation Filtering (Cs+Rf) Approach
In our relaxation filtering approach, we consider a symmetric flux reconstruction using a purely central scheme (CS) combined with a low-pass spatial filter, 7-point stencil Simpson's filter (SF7) in our case, as a relaxation filter (RF). We denoted this solver as CS+RF. For the cell interfacial reconstruction of the conserved quantity, the following symmetric non-dissipative scheme is used [86]: for the interpolation in x-direction, and similarly in y-direction, the conservative interpolation formula reads: where the stencil coefficients are given by: Here, q i,j represents the flow variables (at cell centers) given in Equation (4). The calculated fluxes from the relevant face quantities determined from the nodal values can be used in discretized finite volume equation. In this approach, we assume that the explicit filtering removes the frequencies higher than a selected cut-off threshold through the use of the low-pass spatial filter. A low-pass filter is commonly used in explicit filtering approaches which can be considered to be a free modeling parameter with a specific order of accuracy and a fixed filtering strength [87]. The filtering operation is done at the end of every timestep to remove high frequency content from the solution which eventually prevents the oscillations [58,88,89]. A discussion and analysis of the characteristics on different low-pass filter can be found in [90]. In our investigation, the expression for the sixth-order sequential RF for any quantity f is: where Here, the discrete quantity f i,j yields the filtered valuef i,j and the filtering coefficients are: and σ is a parameter that controls filter dissipation strength in a range of [0, 1] where σ = 0 indicates no filtering effect at all, i.e., completely non-dissipative and σ = 1 indicates the highest filtering effect, i.e., most dissipative with a complete attenuation at the grid cut-off wavenumber. The transfer function of the SF7 filter displays a trend of more dissipation with the increase of the parameter σ [49].

Results
In this section, we present our numerical assessment of the modeling approaches outlined in the previous section for both multi-mode and single-mode two-dimensional RTI test problem. We first illustrate the problem definitions of our test case which is followed by the results obtained by different numerical solvers. We perform our quantitative comparisons between the ILES and CS+RF models using the density contours, density-weighted kinetic energy spectra, and compensated density-weighted kinetic energy spectra plots. For comparative analysis, we obtain the high-resolution ILES and CS+RF solutions by using a parallel computing approach using the Open Message Passing Interface (MPI) framework [91,92]. A detailed discussion on the MPI methodology implemented in our study can be found in [49]. Using both high-and coarse-resolution simulation results, the scaling behaviors of the kinetic energy spectra plots are also investigated in this section.

Two-Dimensional RTI Test Problem: Case Setup
In our numerical experiments, we use a two-dimensional implementation of RTI using the aforementioned numerical schemes. In general, RTI arises at the interface of two fluids when a dense fluid is supported above a comparatively lower density fluid in a gravitational field or stay in the presence of relative acceleration. Since it is found in numerous literature that many properties related to the RTI-induced flows such as the overall growth rate of RTI mixing, dissipation scales, velocity field, and so on, more or less depend on the initial conditions of the flow domain [14,[93][94][95], we consider RTI with multi-mode or randomized perturbation and RTI with single-mode perturbation in our study. We first focus on the case of randomized initial perturbation where our computational domain is set (x, y) ∈ [0, 0.5] × [−0.375, 0.375] with the following initial conditions: p(x, y) = 2.5 − ρy.
Here, L y is set 0.75 and the amplitude of the perturbation is set at λ = 0.01 and α is a random number with a value in between 0 and 1. Since λ is updating itself at each grid point, we note that it is an implicit function of both x and y. On the other hand, for the single-mode RTI, the computational with the following initial conditions: where L x and L y is set 0.5 and 1.5 respectively with the similar amplitude of the perturbation as the multi-mode case, λ = 0.01. The similar two-dimensional RTI test problem set up has been used in various studies related to RTI [96,97]. Figure 1 shows the schematic of the computational domain for both cases of the RTI test problem where it can be seen that we consider the normalized gravity is acting in vertically downward direction in our problem definitions. We must note here that we apply the periodic boundary condition on the left and right boundaries, and the reflective boundary condition on the top and bottom boundaries of our computational domain for both test setup. To get a better understanding on the boundary conditions used in our test setup domain, we can consider an arbitrary two-dimensional domain illustrated in Figure 2. To apply the periodic and reflective boundary condition for our 7-point stencil scheme, we take three ghost points in each direction of the four boundaries of our computational domain. For periodic boundary condition on the left and right boundaries, the ghost point values of the time-dependent variables in vector q (from Equation (2)) can be computed by: where j = −2, −1, ...., N y + 3. On the other hand, our approximation for the reflective boundary condition is as follows: where i = 1, 2, ...., N x . For the parallelization, we do the domain decomposition in the y-direction and update the ghost points of the local domain by transferring the information from the adjacent domain. Although we implement our reflective boundary conditions as defined by Equation (39), we note that the boundary conditions are often applied on the velocity rather than momentum. We stress that simulations of unsteady compressible flows require an accurate control of wave reflections from the boundaries of the computational domain since such waves may propagate from the boundary and interact with the flow [98]. We plan to implement more accurate characteristics-based boundary conditions in our future studies.

RTI with Random (Multi-Mode) Perturbation
The nonlinear evolution of the Rayleigh-Taylor instability from multi-mode initial perturbations is studied based on the density field contour and density-weighted kinetic energy spectra to assess the performance of the underlying modeling schemes. Figure 3 shows the time evolution of the density field at high resolution using the ILES-Roe scheme. As it is shown in [99] that the DNS and ILES give similar results for global properties of RTI-induced mixing, we use the high-resolution results of ILES schemes to avoid the higher computational cost of DNS. Also, it is apparent in Figure 3 that several fine scale structures are captured using the ILES-Roe scheme because of its capability to resolve the smaller scales in high wavenumber region. It can be observed that the mixing growth rate is uniform along the interface at t = 1.6 with multiple modes. It has been seen before in [21] where authors found a uniform growth of the mixing region initially for an idealized initial condition whereas the experimental results of same test condition show the presence of dominant scale at the same time. In our study, even though there are some dominant scales or modes present at t = 4.0, there are a considerable amount of unmixed region can be seen at the same time which indicates the slow mixing rate for this initial condition. In Figure 4), we present the density contour plots for coarse-resolutions obtained by different ILES-Riemann solver combinations at final time of our simulation, i.e., at t = 4.0. As we can see, a clear difference in the growth of scales as well as mixing for both solvers. Since the Rusanov solver is more dissipative than the Roe solver [100], it is expected to have different evolution of the scales in the flow field. Similarly, we can observe different flow field evolution of CS+RF scheme for different filtering strength, σ in Figure 5. Since the higher value of σ adds more dissipation, this solver induces different amount of perturbation in the flow field with the evolution of time than the solver with lower value of σ. Hence, we plot the density-weighted kinetic energy spectra to get a better view in the performance of different solvers [27,[101][102][103][104]. To include these density effects, we define the energy spectrum built on density-weighted velocity vector which can be expressed as: where the density-weighted velocity components are u ρ (x, y, t) = ρ(x, y, t)u(x, y, t), We then can calculate the density-weighted kinetic spectra by following expressions: where k refers to the wavenumber along x-direction. We obtain the Fourier coefficients using a standard FFT algorithm [105]û where i refers to unit imaginary number, and (x i , y j ) determines the Cartesian grid. Since our domain is periodic only in x-direction, we note that our spectra calculations are averaged in y-direction as illustrated in Equation (43). The other statistical measures investigated in our study are the classical kinetic energy spectra and power density spectra. The energy spectra can be calculated using the following definition in wavenumber space [106]: where the velocity componentsû andv can be computed using a similar fast Fourier transform algorithm presented in Equation (44). To quantify the effect of the scale content of density field, we use the power spectrum that reflects the average packaging of density over different scales at any given time in the simulation. This may be given by the following expression: whereρ is the Fourier coefficients of the density field. For the validation of our spectra plots, we follow the well-established theory for two-dimensional RTI systems [14,43,44]. In his seminal paper, Chertkov [43] proposed a phenomenological theory corresponding to the Bolgiano scaling [107] which can be abstracted to k −7/5 scaling law for density or temperature and k −11/5 scaling law for velocity. In Figure 6, we can see that the spectra plots (on the left) obtained by the ILES and CS+RF schemes are showing a clear inertial subrange with the k −11/5 scaling. However, it is apparent that the CS+RF scheme results are the most aligned with the k −11/5 reference line. Moreover, we present the kinetic energy spectra plots without density weighting in Figure 7 which supports the conclusions of the density-weighted energy spectra plots. To validate further, we plot the regular and compensated power density spectra in Figure 8 where it can be seen that the density spectra for CS+RF scheme are following the k −7/5 scaling law.
The time evolution of the spectra shows similar statistical trends for all schemes. Therefore, we will only focus on the results at final time in our subsequent analyses. To compare the dissipation characteristics of the schemes, we place the density-weighted spectra of both ILES schemes in a single plot as well as for both CS+RF scheme with different filtering strength σ in Figure 9. It can be observed that the Rusanov solver is more dissipative than the Roe solver and the higher σ value adds more dissipation. These findings are consistent with the previously found results in the literature as well. Also, the spectra are following the reference k −11/5 scaling. The density-weighted spectra plots compensated by k 11/5 for the CS+RF scheme with different filtering strength show that all the lines are flat above the axis line. However, the kinetic energy spectra plots without the density weighting in Figure 10 exhibit a similar trend as the density-weighted ones. On the other hand, the power density spectra plots in Figure 11 show that the k −7/5 scaling law is maintained for both set of schemes. However, the compensated spectra plots indicate that the CS+RF scheme is more consistent with the scaling law than the ILES schemes. We present another set of density-weighted spectra plot varying grid resolution to show a comparison between the ILES and CS+RF schemes in Figure 12. It is apparent that the CS+RF captures more scales in the inertial subrange than the ILES schemes. However, the CS+RF scheme reaches the effective grid cut-off scales earlier than the ILES schemes. It is because the CS+RF solvers do the filtering once at the end of the simulation whereas the ILES solvers implicitly adds dissipation throughout the simulation. As a result, ILES schemes capture wide range of scale at high wavenumber even though they resolve comparatively less scales in the inertial subrange. Some key points can be seen from Figure 12 that the σ = 1.0 solver is the most dissipative among all solvers considered in this study, and σ = 0.4 solver captures more scales in the inertial subrange than the other solvers. Also, ILES-Roe solver resolves the highest range of scales in high wavenumber for both coarse and high resolution which explains the appearance of very fine small-scale structures in the density field contour plot obtained by ILES-Roe solver.  (c) (d) Figure 6. Time evolution of density-weighted kinetic energy spectra and compensated densityweighted kinetic energy spectra for the RTI problem with multi-mode perturbation obtained using different modeling approaches at a resolution of 16384 × 24576; (a) density-weighted spectra using ILES-Roe solver; (b) compensated density-weighted spectra using ILES-Roe solver; (c) density-weighted spectra using CS+RF (σ = 1.0) solver; (d) compensated density-weighted spectra using CS+RF (σ = 1.0) solver. (c) (d) Figure 11. Comparison of ILES (ILES-Roe and ILES-Rusanov) models and CS+RF (σ = 1.0 and σ = 0.4) models for the RTI problem with multi-mode perturbation showing the power density spectra and compensated power density spectra at different resolutions; (a) power density spectra using ILES solvers; (b) compensated power density spectra using ILES solvers; (c) power density spectra using CS+RF solvers; (d) compensated power density spectra using CS+RF solvers.

RTI with Single-Mode Perturbation
Numerical simulation of flows with RTI is comparatively challenging because the instability grows from the small scales of the flow field [19,108]. Since the analytical modeling can be done for single-mode RTI, the numerical study of the RTI with the single-mode perturbation setup has been started very early [33,109] and still being studied extensively to understand and explain the nature of RTI-induced flows [20,31,32,[110][111][112][113]. For our analyses of the single-mode perturbation case, we first present the time evolution of the density field results obtained by ILES-Roe solver at a high resolution of 8192 × 24576 in Figure 13. It is observed in many studies that the tips of the spikes of a single-mode RTI-induced flow always maintains the symmetry [14]. Yet in our simulation, the line of symmetry within the spike of the single-mode RTI in Figure 13 is seen broken at time t = 4.5. This phenomenon was observed and well explained by Ramaprabhu et al. [31] at late-time of the RTI flow simulation which the authors referred to as "chaotic mixing" at the late-time regime. With simulations of the Euler equations, it can be also seen in [96] that the less dissipative schemes show this interface breaking up, while the more dissipative schemes suppress the instability. In our high-resolution simulation, the presence of small-scale structures can be seen in very early stage which lead to secondary instability, i.e., the KH vortex formation as well as chaotic mixing. However, for small or modest grid resolutions, the numerical viscosity suppresses the small-scale structures and preserves the symmetry which can be seen in Figures 14-17. In Figure 14, we present the state of the density field at t = 2.7 (top row) and t = 4.5 (bottom row) for the simulations by using the ILES-Rusanov solver at different grid resolutions. It is apparent that the 256 × 768 and 1024 × 3072 resolution results hold the symmetry. But the 4096 × 12288 result shows the development of smaller scales at t = 2.7 which leads to the loss of symmetry at final time, t = 4.5. Similar conclusions can be made for the ILES-Roe solver results in Figure 15. However, the loss of symmetry can be observed even in 1024 × 3072 resolution simulation for the ILES-Roe scheme since the ILES-Roe solver is less dissipative compared to the ILES-Rusanov solver. For both CS+RF schemes in Figures 16 and 17, the symmetry holds for lower resolutions and breaks for higher resolution. Since there is no physical viscosity in Euler simulations, we note that the breakup of the interface and the loss of symmetry can be due to the numerics. When increasing the resolution, the loss of symmetry in RTI problems have been also demonstrated in the literature (e.g., see [55,114,115]). Similar observations can be seen when we use the higher-order numerical schemes. We also refer to [116,117] for an illustration of symmetry breaking and increasing mixing in Richtmyer-Meshkov instability problems for solving Euler equations. Figure 18 presents the density field plots at the final time, t = 4.5, to get a comparative idea between the performance of the ILES schemes. We can observe in Figure 18 that the symmetry is maintained in lower resolutions, but starts to break with the increase of the resolution for both ILES solvers. If we look at the 4096 × 12288 resolution results for both ILES solvers, we can see that the ILES-Roe scheme result is more deviated from the symmetry than the ILES-Rusanov scheme because of the dissipative behavior of the ILES-Rusanov scheme. Based on these findings, we can say our two-dimensional simulation results are consistent with the findings in [31] for three-dimensional RTI case. Additionally, the dimensionless Atwood number defined as: is set as A ∼ 0.33 in our case, and it is lower than 0.6, which indicates the formation of reacceleration phase in the flow field due to the secondary KH instabilities. As suggested in the literature [31], these secondary instabilities can be responsible for the change in the usual behavior of the spikes in single-mode RTI flows. These findings are also supported by the works of Liska and Wendroff [96] where they showed that the less dissipative schemes result in an interface break up while the instability might be suppressed by high dissipative schemes. The same observations can be found at final time in Figure 19 that the higher resolution results start to break the symmetry of the spike for both CS+RF schemes. However, the density fields obtained by the CS+RF and ILES schemes seem different due to different amount of dissipation added to the system by different solvers which would eventually lead to different evolution of the flow fields. To get a more precise understanding on the simulation results, we next focus on the density-weighted kinetic energy spectra plots. We can say from the time evolution of the kinetic energy spectra plot in Figure 20 that the trends of the spectra are similar at late-stage of the simulation. We can observe that the density-weighted spectra analysis for ILES-Roe scheme shows an inertial subrange following k −11/5 scaling law. On the other hand, the kinetic energy spectra for CS+RF scheme follow the k −11/5 scaling in Figure 21. To validate our findings further, we present the power density spectra plots for both ILES-Roe and CS+RF (σ = 1.0) schemes in Figure 22. The power density spectra display a good alignment with the k −7/5 reference line. Since the time evolution of the field for both schemes follow a similar trend, we can consider the solutions at final time for rest of our analysis.  Figure 23 shows that the ILES-Rusanov solver is more dissipative than the ILES-Roe solver as expected and the CS+RF scheme with σ = 1.0 is more dissipative than the σ = 0.4 solver. One interesting point can be noticed that the density-weighted spectra for the CS+RF scheme tend to deviate from the reference k −11/5 line at high wavenumber; however, the kinetic energy spectra in Figure 24 and the density-weighted spectra in Figure 25 clearly show that the spectra for the CS+RF scheme follow the reference scaling laws. On the other hand, the ILES spectra also maintain the inertial subrange following the k −11/5 and k −7/5 laws. Finally, similar to the previous section of multi-mode RTI case, we find the CS+RF solver captures more scales in the inertial range than the ILES solvers as shown in Figure 26. However, the ILES solvers resolve more scales in the high wavenumber region. This explains the reason we have seen different density field evolution for different solvers and appearance of smaller scales in ILES solvers than the CS+RF solvers. Since the ILES-Roe solver is least dissipative among the other solvers, we can see in the density contour plots that the ILES-Roe solution deviates most from the symmetry.   (c) (d) Figure 20. Time evolution of density-weighted kinetic energy spectra and compensated densityweighted kinetic energy spectra for the RTI problem with single-mode perturbation obtained using different modeling approaches at a resolution of 8192 × 24576; (a) density-weighted spectra using ILES-Roe solver; (b) compensated density-weighted spectra using ILES-Roe solver; (c) density-weighted spectra using CS+RF (σ = 1.0) solver; (d) compensated density-weighted spectra using CS+RF (σ = 1.0) solver.
(a) (b) (c) (d) Figure 21. Time evolution of kinetic energy spectra and compensated kinetic energy spectra for the RTI problem with single-mode perturbation obtained using different modeling approaches at a resolution of 8192 × 24576; (a) kinetic energy spectra using ILES-Roe solver; (b) compensated kinetic energy spectra using ILES-Roe solver; (c) kinetic energy spectra using CS+RF (σ = 1.0) solver; (d) compensated kinetic energy spectra using CS+RF (σ = 1.0) solver.

Summary and Conclusions
In this paper, we put an effort to show the performance of a relaxation filtering approach using central scheme (CS+RF) on resolving the flows resulting from Rayleigh-Taylor hydrodynamic instability, and compare the simulation results with the results obtained by two common ILES-Riemann solver schemes. To assess the performance of the solvers, we use the density field contours and different spectra plots. We further analyze the resolution capacity of both CS+RF and ILES schemes as well as the flow nature at high-resolution simulation using CS+RF and ILES schemes. To validate the observations from the field plots, we use the statistical tools, i.e., kinetic energy and power density spectra plots for both high and coarse resolutions which show consistency with the existing results in the literature. In our investigation, we consider the two-dimensional RTI test problem with two different initial conditions. From the simulation results of both cases, we come to this conclusion that the CS+RF schemes capture more scales in the inertial subregion whereas the ILES schemes resolve a wide range of scales in high wavenumber region. The ILES-Rusanov scheme is more dissipative than the ILES-Roe scheme, and hence, the ILES-Roe scheme tends to deviate more from symmetry in the spike of single-mode RTI case. On the other hand, it is also observed that the dissipation can be controlled by σ parameter for CS+RF scheme which also affect the perturbation as well as the evolution of the flow field. Furthermore, we observe that the kinetic energy spectra follow k −11/5 scaling law for both multi-mode and single-mode RTI case whereas the power density spectra plots are seen to be more align to k −7/5 line at different resolutions. Also, we observe a chaotic mixing at late-time stage for single-mode RTI case at high resolution. It is because of the formation of secondary KH instabilities at high-resolution simulation of single-mode RTI case. The higher numerical dissipation due to the coarser resolution suppresses the formation of secondary instabilities which is the reason behind the preservation of symmetry for our coarse-resolution simulation results. Overall, we believe the study of relaxation filtering approach using CS will be a good contribution to the numerical study of RTI-induced flows as well as for understanding the nature of the flow field with the evolution of the instability.