# How to create a fluid animation like this?

When two fluids move past each other at different speeds, complex instabilities arise that look really cool:

The Navier-Stokes equations are numerically intense but can usually predict what will happen, is there any good way to plot these bi-fluid fractals with Dynamics?

These are related questions, but don't show how to generate an efficient Animate[] or GIF:

We cannot accurately reproduce what is shown in the animation, since the initial data are unknown. We can show a similar mixing process at the boundary of two oncoming flows. We use FEM and the NSE integration method described in Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM.

Calculation of the velocity:

Needs["NDSolveFEM"]
k0 = Pi; L = 5; eps = 10^(-3);
reg = Rectangle[{0, 0}, {2*L, L}];
k = 70; ν = 5/10^3; t0 = 1/400; f = TranslationTransform[{2*L, 0}];
bcP = {
DirichletCondition[{u[x, y] == -1, v[x, y] == 0}, y == 0 && eps <= x <= 2*L - eps],
DirichletCondition[{u[x, y] == 1, v[x, y] == 0}, y == L && eps <= x <= 2*L - eps],
PeriodicBoundaryCondition[u[x, y], x == 0, f],
PeriodicBoundaryCondition[v[x, y], x == 0, f]};
UX[0][x_, y_] := Tanh[k0*(y - L/2)];
VY[0][x_, y_] := 0.3*Sin[k0*x]*Sin[Pi*(y/L)]^2;
P0[0][x_, y_] := 0;
Do[
{UX[i], VY[i], P0[i]} =
NDSolveValue[
{
{
Inactive[Div][{{-μ, 0}, {0, -μ}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + Derivative[1, 0][p][x, y] + (u[x, y] - UX[i - 1][x, y])/t0 + UX[i - 1][x, y] D[u[x, y], x] + VY[i - 1][x, y] D[u[x, y], y],
Inactive[Div][{{-μ, 0}, {0, -μ}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + Derivative[0, 1][p][x, y] + (v[x, y] - VY[i - 1][x, y])/t0 + UX[i - 1][x, y] D[v[x, y], x] + VY[i - 1][x, y]*D[v[x, y], y],
Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]
} == {0, 0, 0} /. μ -> ν,
bcP, DirichletCondition[p[x, y] == P0[i - 1][x, y], x == 2*L || x == 0]
},
{u, v, p},
Element[{x, y}, reg],
Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.02}}
],
{i, 1, k}
];


Interface function diffusion calculation:

cn[0][x_, y_] := Tanh[k0 (y - L/2)];
Do[cn[i] =
NDSolveValue[{.3 ν (D[c[x, y], x, x] +
D[c[x, y], y, y]) == (c[x, y] - cn[i - 1][x, y]) +
UX[i - 1][x, y] D[c[x, y], x] + VY[i - 1][x, y] D[c[x, y], y],
DirichletCondition[c[x, y] == -1, y == 0 && eps <= x <= 2 L - eps],
DirichletCondition[c[x, y] == 1, y == L && eps <= x <= 2 L - eps],
PeriodicBoundaryCondition[c[x, y], x == 0, f]},
c, {x, y} ∈ reg,
Method -> {"FiniteElement", "InterpolationOrder" -> {c -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.02}}];, {i, 1, k}]


Create frames and animations:

frame = Table[
ContourPlot[cn[i][x, y], {x, 3, 6}, {y, 1.75, 3.25} ,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
Frame -> False, PlotRange -> All,
Contours -> Table[j^2 Sign[j], {j, -1, 1, .1}],
ContourShading -> Automatic, ContourStyle -> Green], {i, 0, k, 1}];

Export["C:\\Users\\name\\Desktop\\cnu03_9.gif", frame, AnimationRepetitions -> Infinity]