How to create a fluid animation like this?

Question

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:

• Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM
• Nonlinear FEM Solver for Navier-Stokes equations in 2D

Answer 1

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["NDSolve`FEM`"]
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]