CFD Instability Examples

These instability browser simulations are meant to showcase how different conditions can affect the development of an instability, both in real life and in simulations. They are not meant to be completely physically correct (do NOT use these simulations as results for anything!). Detailed three-dimensional simulations of these instabilities can require substantial computing resources. These simplified examples make some of the same questions accessible in a browser. Having access to this as a pedagogical tool can be really helpful in explaining what researchers care about and what they need to look for. The rest of this accordion section will be a slightly more technical explanation of what is happening. The code can be found in the website’s source repository. Please let me know if you encounter unexpected behavior or crashes. Enjoy!

How the calculations work

The simulations follow the motion of two materials in a two-dimensional domain using an Eulerian description. Both materials share a velocity field, while the relative amount of each material determines the density throughout the domain. The calculations use dimensionless quantities, allowing changes in density contrast, acceleration, and initial conditions to be compared without choosing a particular physical scale. The domain is divided into a fixed grid of cells or zones. Each cell stores the proportion of each material it contains, while velocities are stored along the cell edges. This arrangement makes it possible to calculate how fluid moves between neighboring cells. (Staggered MAC grid)

Material description

The materials begin on opposite sides of a wavy interface, or perturbation. Its shape is set by the number of waves and the amplitude relative to their wavelength. Across the interface, the material fraction changes smoothly over a small distance, which we refer to as a diffuse interface. The width of this transition can affect how the instability develops and how the materials become interleaved. Intermediate colors represent material fractions within a grid cell. The Atwood number describes the density difference between the materials relative to the sum of their densities. The Rayleigh–Taylor simulation begins at rest. With the heavier material above the lighter material, downward acceleration drives the instability. The Kelvin–Helmholtz simulation begins with the upper and lower layers moving in opposite directions. A smooth velocity transition connects the layers, and a small velocity disturbance helps initiate the rolling motion.

Eulerian flow and methods

The solver, the numerical code that advances the simulation, updates the material distribution by calculating how much material enters and leaves each cell. These exchanges depend on the velocity at each cell edge and the nearby material fractions. A monotonized central slope limiter reduces artificial oscillations around steep changes in the material distribution. (An illustrated explanation of slope limiting) In order to move the material distribution forward in time, an explicit two-stage solver is used. The first stage predicts an updated distribution. The second stage advances that prediction again and averages the result with the original distribution. (SSP-RK2 time integration, also described as Heun’s method)

Velocity is transported by tracing backward through the flow to estimate where the fluid came from during the time step. A midpoint estimate is used when tracing this path, and the velocity is interpolated from nearby grid values. (Semi-Lagrangian advection) Acceleration then changes the vertical motion according to the local density. When viscosity is enabled, an additional update smooths velocity differences between neighboring cells. Larger viscosity values increase this smoothing. (Body forces and viscous diffusion)

The fluids are treated as incompressible, so the model does not include compression or sound waves. This approximation is useful for studying flows where compressibility effects are small. After updating the velocity, the solver calculates a pressure correction that brings the flow entering and leaving each cell into balance. The local material density is included when calculating this correction. (Pressure projection) The pressure calculation uses a sequence of coarse and fine grids. Coarser grids help correct errors that extend across the domain, while finer grids resolve more local differences. The solver repeats this process and checks the remaining imbalance. (Multigrid methods)

The calculation time step changes as the simulation develops. Flow speed, grid spacing, acceleration, and viscosity are used to determine how far the solver can safely advance. Faster motion and finer grids generally require smaller time steps. (CFL-based time stepping) The save frequency sets the spacing between recorded snapshots. The solver takes calculation steps as needed to reach each scheduled output time.

For both simulations, the left and right edges of the domain are connected, so material leaving one side enters the other. The top and bottom act as stationary walls. Fluid can move along these walls but cannot pass through them. (Periodic and free-slip boundary conditions)

Why this runs well in a browser

The simplified two-dimensional model keeps the computational work manageable, and WebGL 2 allows the calculations to run on the graphics processor. The fluid fields are stored in textures, allowing many grid cells to be updated in parallel. It is honestly wild to me how well this works. (GPU fluid simulation) The displayed colors come from the calculated material fractions. Each endpoint color represents one material, while intermediate colors represent cells containing portions of both. Changing the color palette updates the display without changing the simulation.

Major inspiration for being able to do is credited to this website. While I have wanted to this for a couple years prior to the initial posting of these apps, most of my attempts resulted in bad results. GPU acceleration was a pretty big epiphany for getting this done.