2-D test problem 3: The Kelvin-Helmholtz Instability
Linear regime:  Two layers of fluid sliding past each other with relative speed 2v0 form a shear layer (Fig. 1) and is subject to the Kelvin-Helmholtz instability for certain ranges of the relative Mach number, 2M, where M = v0/cs, and where cs =
(gp0/r0)1/2 is the sound speed.  To approach this problem analytically, one examines perturbations to the linearised fluid equations, from which the following conclusions may be drawn:

  

1.  For 2M > 23/2 ~ 2.83, the shear layer is unconditionally stable to perturbations;

      

2.  The shear layer is most unstable for 2M = 31/2 ~ 1.73;
   
 
 
Figure 1: A shear layer
 

3.  Where the shear layer is unstable, the growth rate of the perturbation is inversely proportional to its wavelength.

       

The latter conclusion cannot be true indefinitely; at some point, the effects of viscosity and/or surface tension damp out the perturbations before they can grow.

        

ZEUS-3D simulations were performed on a 500 x 500 uniform grid resolving the domain (x,y) = (-0.5:0.5,-0.5:0.5), with an additional 500 (uniform) x 150 (geometrically ratioed) zones in each of two "buffer regions" (x,y) = (-0.5:0.5,-10:-0.5) and (x,y) = (-0.5:0.5,0.5:10) to insulate the region of interest from unwanted reflections from the upper and lower boundaries.  Others have used reflecting boundary conditions at y = +/-0.5 with success (e.g., Ryu et al., 2000, ApJ, 545, 474), but this works only for relatively low Mach numbers (2M < 1.2).  At higher speeds, shocks are excited and much more of the region above and below the shear layer become dynamically important.  Periodic boundary conditions were imposed on the left and right boundaries.

For a variety of relative Mach numbers, a shear layer of thickness 0.04L (where L = 1 is the width of the box) was initialised with r0 = 1, cs = 1, and perturbed with a sinusoidal transverse velocity with amplitude 0.001 and a specified wavelength.  The non-zero thickness of the shear layer is to prevent minute numerical perturbations from breaking up a too-thin shear layer into small vortex tubes.  Figures 2 and 3 below show contour images of the vorticity component normal to the image plane (particularly effective at visualising the shear layer) for many of these simulations before the onset of non-linearity so that results may be compared directly with conclusions from linear theory.
 
 
Figure 2.  ZEUS-3D images showing the development of the K-H instability after t = 3 sound-crossing times for the relative Mach numbers indicated.  Highest growth rate is for 1.25 < 2M < 1.5, whereas linear theory predicts highest growth rate at 2M = 1.73.  This discrepency may result from the non-zero thickness of the numerical shear layer, as opposed to infinitesimally thin shear layer assumed by linear theory. Not shown are 2M = 2.5 and 2.75 which grow too slowly to be seen at t = 3 (but develop by t = 12), and 2M = 3.0 which remains stable through t = 12. Linear theory predicts stability at 2M = 2.83.
 
 
 
Figure 3.  ZEUS-3D images showing the development of the K-H instability for 2M = 1.5 at t = 0.75 sound-crossing times for the wavelengths shown.  L = 1 is the length of the box.  These results agree qualitatively with the linear result that the growth rate is inversely proportional to the wavelength. Note that non-linear effects are already apparent in the l = L/4 case.
Non-linear regime:  The thumbnails below link to animations of the normal vorticity for the indicated relative Mach numbers of the shear layer, and carry the simulations to t = 12, well beyond when linearity ceases to be valid.  The sense of the shear is reversed in the animations compared with the contour images above.  Thus, the top layer is moving to the left in the animations, and the bottom layer is moving to the right.
 
2M = 0.25
2M = 0.50
2M = 0.75
2M = 1.00
2M = 1.25
2M = 1.50
                
                
2M = 1.75
2M = 2.00
2M = 2.25
2M = 2.50
2M = 2.75
2M = 3.00
 
The formation of a "cat's eye" (a strong and sustained vortex tube prominant near the end of most of the animations) is a well-known effect and has been observed for years in simulations peformed by finite-difference codes (it is worth noting that most SPH codes cannot reproduce the K-H instability in general, and the cat's eye in particular).  It is the end state for all unstable Mach numbers, with the most rapid growth at 2M ~ 1.5.  For 2M > 1 (supersonic relative speeds), shocks are excited above and below the shear plane which cause much more of the fluid to become involved.  These animations have been "bracketed" so that crimson means zero vorticity and the deepest blue means the most negative value attained (always along the original shear layer at t = 0).  "Blood-red" indicates positive vorticity only found in the shocks.
Periodic boundary conditions: Periodic boundary conditions in dzeus35/36 are accurate to machine round-off error.  Unlike an unsplit scheme where boundary conditions must be applied only at the end of each MHD cycle, boundary conditions must be applied frequently during each MHD cycle in ZEUS-3D owing to its operator and directionally split methods.  A significant effort has been made to ensure that boundary conditions are applied as often as needed (but only as often as
needed) so that periodicity is maintained to machine accuracy.  The same is true for all flavours of the reflecting boundary conditions.

The thumbnails to the right are linked to mpegs of half-resolution simulations of the K-H instability at 2M = 1.0 with two periods.  To machine accuracy, each period is identical to the other in all variables, indicating the periodic boundary condition behaves exactly as the central column of zones.  Shown are both vorticity (with a reversed colour palette compared to the animations above) and density using a Courant number of 0.50.
     
vorticity
density
 
An unexplained numerical instability: Animations linked to the thumbnails to the right are from simulations identical to those immediately above, except with a Courant number of 0.75.  Oscillations superposed over what appear to be the correct solutions are apparent in all variables, particularly density, and as yet I have no explanation for this.  The instability is "cured" by taking a lower Courant number (in this case, 0.5), though I do not think this is a normal CFL violation since the amplitude of the oscillations stabilises at several percent of the underlying solution.  CFL violations typically grow without bound.  In these simulations, the wavelength is nine zones, regardless of the total number of zones across the box.
  
vorticity
density
The higher the relative shear velocity, the lower the Courant number needed to stabilise the flow.  At 2M = 1.5, a Courant number of 0.35 is needed for stability and at 2M = 2.5, 0.25.  Increasing either form of the artificial viscosity can also cure the problem, but I am convinced it is only because this also decreases the time step.  Thus, I do not believe this is a problem in upwindedness (which the artificial viscosity is supposed to stabilise), though I have not ruled this out.  It is clearly hydro in nature, since there is no magnetic field in these simulations, though curiously it does resemble the striping observed in advecting a magnetic flux loop, which is also cured by lowering the Courant number (problem 4 in the 2-D Gallery).

Jim Stone agrees that this doesn't seem like a CFL violation nor an upwindedness problem, and has suggested this instability may be related to the operator split nature of the momentum equation, as he does not see this behaviour in ATHENA.
Kelvin-Helmholtz instabilities in nature: The atmosphere surrounding a rapidly rotating planet is a particularly good natural laboratory for K-H instabilities since the coriolis effect is forever driving shear layers with relative speeds comparable to the rotation speed of the planet.  On the Earth, K-H instabilities are occasionally observed in high clouds (Fig. 4).  However, with its higher rotation speed, K-H instabilities are common-place in the Jovian atmosphere with the Great Red Spot (observed since Galileo's time) being the most spectacular example of a sustained cat's eye in the solar system (Fig. 5).
 
     
     
Figure 4.  Clouds often form along shear layers (separating air masses of differing moisture and temperature), and are thus subject to the Kelvin-Hlemholtz instability, particularly for sufficiently high shear velocities (photo credit: Brooks Martner, NOAA/ETL).
  
Figure 5.  Jupiter's Great Red Spot is a vortex tube (cat's eye) suistained by transonic shear layers in the planet's atmosphere.  Note the smaller-amplitude K-H instabilities in the surrounding gas, as well as within the Red Spot itself (photo credit: NASA).