The Rayleigh-Taylor instability is a well-known benchmarck for CFD codes. The idea is to start with two phases, on on top of the other, the lighter one being underneath. The interface is slightly perturbed, and this plume appears. I describe a quick and dirty way of getting this instability.

The usual computation is as follows:

- The simulation cell is a rectangle, dimensions . I.e. an aspect ratio of 4, and a 2D case.
- The density ratio is (I think) set to 3. I.e. the lighter phase is three times as light as the heavy one. I think the actual values are “1” and “3” in whichever reduced units that are used. Atwood number is therefore 1/2.
- Fluids are initially at rest, so there is no input typical velocity. A natural one would be , similar to the velocity that a fluid would have if it falls a distance
*d*under gravity. - The Reynolds number is therefore fixed as , where (I think), the lighter density and viscosity are used.
- Boundary conditions for the velocity are: no slip at the top and at the bottom, and slip (this is quite important) at the left and right walls. Pressure: zero normal gradient at all walls. There may be a symmetry plane at the
*x*=0 line, in order to avoid half of the domain. - The initial interface is perturbed by setting it as a cosine shape:
- Lots of details can be found on Guermond & Quartapelle (2000). A projection FEM for variable density incompressible flows. Journal of Computational Physics, 165(1), 167-188. It’s easy to get its pdf.

OK, this is all quite doable in OpenFOAM. As explained in this Eric Paterson at a workshop in Chalmers, the only trick seems to be to use funkySetFields. This utility is now part of swak4Foam, the Swiss Army Knife for FOAM. Now, installation for new releases of OpenFOAM may be a bit tricky, but the procedure is well documented.

Here’s the thing: with my particle simulations I still don’t know how to perform multiphase simulations: different densities, and viscosities, etc. I do know how to carry a color field around with particles: you just set it at time zero and never change it. O top of that, I can only do periodic boundary conditions, and on a square. So, this is what I did

- The two fluids have the same physical parameters
- But, some funny gravity acts upwards for one fluid and downwards for the other. If the color function is , either 0 or 1 for the two phases, this force per unit mass would be
- Since I can only do a square, I perturb my interface as . That should give us four plumes, instead of one

Now, that I can do! Also, it is not hard to hack interFoam to do the same. Just take the geometry from, say “cavity”, set cyclic patches at the boundary. Then, on the code, make sure “p” and not “p_rgh” is the field that enters the velocity equation, then add an extra term on the right-hand side, **+ (2 * alpha1 – 1) * rho * g** (just like that!). Get the control files from, e.g. dam break, and you are all set to run.

Btw, my method is not yet published, but it’s similar to SPH, MPS, or pFEM.