A shear flow
Channel: Nils Berglund
Date: 2022-10-26

In a shear flow, several layers of fluid move at different speeds, causing vortices (or eddies) to form where the layers touch each other. This simulation starts with horizontal layers, whose velocity has a small dependence on the x-coordinate. Colors depend on the vorticity of the velocity field, which measures how much it rotates. In the second part, a logarithmic scale is used. This simulation is strongly inspired by the description on the page http://ammar-hakim.org/sj/je/je13/je13-incomp-euler-2d.html
The incompressible Euler equations describe the velocity field of a non-viscous, incompressible fluid. They are strongly non-linear, because of the presence of a so-called advection term in what is known as the material derivative. Another difficulty in solving these equations numerically is the incompressibility condition, which requires the velocity to remain divergence-free. Since most numerical schemes will not conserve the divergence, a common approach is to "project" the solution on the space of divergence-free flows after each integration step. 
I used here a different approach, based on the so-called stream function. By writing the velocity as the curl (rotational) of a stream function, directed along the z-axis perpendicular to the plane of the two-dimensional fluid, one ensures that the velocity is always divergence-free. The incompressible Euler equations are then equivalent to a partial differential equation for the stream function, and a Poisson equation linking stream function and vorticity. Since I had no solver for the Poisson equation at hand, I used a forced heat equation for the evolution of the vorticity, which admits the solution of the Poisson equation as a stationary solution. This induces additional errors, which are hopefully not too important. There is also a weak dissipation term, proportional to the Laplacian of the stream function, to prevent numerical blow-up. 
The initial state for the stream function has the form delta*cos(ax) + g(y), where g(y) = 1/(rho*cosh((y+b)/rho)^2) for positive y, and g(y) = -1/(rho*cosh((y-b)/rho)^2) for negative y. There are periodic boundary conditions on the shown rectangular simulation region. 

Render time: 28 minutes 44 seconds
Color scheme: Turbo, by Anton Mikhailov
https://gist.github.com/mikhailov-work/6a308c20e494d9e0ccc29036b28faa7a

Music: At The Shore - The Dark Contenent by Kevin MacLeod is licensed under a Creative Commons Attribution 4.0 licence. https://creativecommons.org/licenses/by/4.0/
Source: http://incompetech.com/music/royalty-free/index.html?isrc=USUAN1100770
Artist: http://incompetech.com/

See also https://images.math.cnrs.fr/Des-ondes-dans-mon-billard-partie-I.html for more explanations (in French) on  a few previous simulations of wave equations. 

The simulation solves the 2D incompressible Euler equation for the stream function and vorticity. The solution of the Poisson equation for the vorticity is approximated by the solution of a forced heat equation, having the solution of the Poisson equation as stationary solution.
C code: https://github.com/nilsberglund-orleans/YouTube-simulations
https://www.idpoisson.fr/berglund/software.html
Many thanks to Marco Mancini and Julian Kauth for helping me to accelerate my code!