This 2D L&W Rayleigh Taylor (RT) test consists of a domain with a constant gravity field with uniform acceleration.

• t = 8.5
• Adiabatic index: Gamma = 5/3
• Acceleration: a = -0.1 (downwards)

Hyrdrostatic initial conditions and compressibility

The grid is set up in hydrostatic equilibrium. With a constant gravity field this implies an exponential vertical distribution, with a scale height for each component depending on the isothermal sound speed a^2 = P/d.

Although L&W do not specify the isothermal sound speed, or alternatively the ratio of pressure on density at the interface, it appears that the ratio P/d was approximately 1.0 at the interface. This results in a dynamically cool system, with some compressibility. The density in the upper, colder, layer is distinctly exponential when viewed with the contoured colour table used here, although it is less obvious in the L&W figures. This affects the effective Attwood number : that is commonly used to describe the RT instability, in that the density contrast effectively changes as the bubble/stalk structures grow, whereas in the incompressible limit the upper and lower densities remain constant.

The appearance of the L&W results compared with the front-tracking results they presented from the Los Alamos code seems to indicate that the Los Alamos curve may be from an incompressible simulation, as the extent of the bubble vortices and the x-height of the lowest part of the dense stalks differ somewhat.

The initial tests here use P/d = 1.0, and give the best comparison to the L&W codes. Additional hotter simulations. with P/d = 10.0 and P/d = 100.0 giving successively less compressible fluids, approach the Los Alamos solution.

Critical setup of the initial interface

The well known property of the RT instability is the exponential initial growth rates, with the fastest rates corresponding to the highest wavenumbers. Consequently, any minute deviations from a smooth initial conditions give rise to high wavenumbers. perturbation which will grow rapidly. High resolution methods will model this property best. As the sinusoidal wave is averaged over a finite grid, even with some smoothing of the interface, the high resolution methods will break up into a complex multi-mode interaction between the dense and low density regions.

Here the sinusoidal interface between the high and low density regions is convolved with a Gaussian FWHM of 1.0 cell, sampled over ±3.16 sigma with 1024x1024 subcells. With any less smoothing the initial pixel aliasing causes unwanted structure to evolve uncontrollably.

Suppression of high frequency modes

With sufficient numerical diffusion, the higher modes can be suppressed, leaving just the principle mode to grow. This is reflected in the results of L&W, where the diffusive central methods produced single RT bubble/stalk structures, while PPM, WENO, and VH1 broke up into complex structures.

Here Fyris performs like the other high resolution methods, with small initial high wavenumber deviations from a continuous sin wave interface that grow and eventually distort the final state.

However, Fyris has the option to increase the numerical diffusion incrementally to deliberately suppress the unwanted high modes, while leaving the principle mode to grow. The process used to flatten the interpolation parabola, to reduce the order of the method in the presence of shock waves in order to satisfy the Godunov theorem, is used to increase the diffusion.

The minimum degree of flattening is normally 0.0, for no flattening, giving the maximum available order to the interpolation. A flattening value of 1.0 represents a 100% flattening of the mean slope, so the left, central and right values are equal, and the method reduces to first order Godunov fluxes.

Normally the flattening is assigned dynamically, as described in the PPM method. To add global diffusion in Fyris, the minimum flattening can be raised from 0.0, normally in the range of 10-30% flattening. This means that the slope from the conservative interpolation of the left and right values is reduced by 10-30%, while retaining parabolic interpolation between the less extreme endpoints. As flattening changes from 0.0 to 1.0 the order is continuously reduced from 3rd to 1st order spatial interpolation.

In the RT test, minimum flattening of 0.1, 0.2 and 0.3 are trialled, and the results compared with the high resolution runs.

Square and rectangular cells

The L&W test was defined on a 100x400 grid, with coordinates which range 0-1/6, 0-1. This results in cells with a rectangular domain, with 50% higher resolution in the x direction. This raises the possibility of anisotropic smoothing or numerical diffusion. Here a set of simulations. with square cells, 100x600 cells are tested, giving the same resolution in the y direction as the original L&W test in the x-direction.

The code and configuration files plus comparison data for Fyris to run the 2D Rayleigh Taylor problem will be available soon.