8
$\begingroup$

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

enter image description here

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:

$\endgroup$
5
13
+100
$\begingroup$

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]

Figure 1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.