Message Boards Message Boards

Rogue wave: stable fluids algorithm for air/water interface simulation

SUPPLEMENTARY WOLFRAM MATERIALS for the ARTICLE:

Sergio Manzetti, Alexander Trounev (2021).

A Navier-Stokes model for Rogue wave simulation.

ResearchGate, Technical Report. https://www.researchgate.net/publication/354527324


enter image description here

ABSTRACT

Rogue waves are anomalous phenomena occurring over large sea bodies, where they are said appear out of nowhere and disappear without a trace. Simulations of rogue waves have been carried over the last 40 years, with several models published. In this paper we investigate the formation of rogue waves by two models: I. the Navier-Stokes equation where we take into account density and viscosity gradient on the air-water interface due to diffusion and II. a Korteveg de Vries-type model for quantum jumps which we developed earlier. We derive also a stable fluids algorithm which we use to compute nonlinear waves interaction produced by the wind. The results are discussed and compared.

This code solves problem of viscous incompressible flow with gravitational force in a rectangle with periodic boundary condition on the left and right side and with Dirichlet condition on the top and bottom side. In the initial condition fluid velocity is periodic wave, and density has unit step like distribution on the air/water interface. Some other application of this code has been discussed on https://mathematica.stackexchange.com/questions/246091/stable-fluids-code-for-electromagnetic-mixture-application

This model has been discussed in our report A Navier-Stokes model for Rogue wave simulation

Two phase model of air-water interface

Let us consider the system of equations describing nonlinear waves on the air-water interface. As it well known the air and the water can be considered as viscous incompressible fluids with density $\rho_a, \rho_w$, and dynamical viscosity $\mu_a, \mu_w$ consequently. Taking into account the gravity force and diffusion on the air-water interface, we have $$\nabla.\bf{u}=0\\ \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}.\nabla)\mathbf{u}+\frac{\nabla P}{\rho_i}=\nu_i\nabla ^2\mathbf{u}+\mathbf{f}\\\ \frac{\partial \phi}{\partial t} +(\mathbf{u}.\nabla)\phi=\frac{\nu_i}{Sc_i}\nabla ^2\phi\\ $$ Here it is indicated $\rho_i, i=a,w $ - air and water density, $\mathbf{u}=(u_x,u_y)$ - flow velocity, $P$ - pressure; $\bf{f}$ - force acting on the volume of the air-water mixture; $\phi$ - the interface function describing the averaged density of the air-water mixture; $Sc_i$ - analog of the Schmidt number characterized water diffusion in the air and air diffusion in the water. Note that we considered averaged effect of mass transfer on the air-water interface including drops and bubbles. Let us define the Cartesian coordinate system so that the $y$ axis is directed against the direction of the gravitational acceleration vector and the $x$ axis is parallel to horizon. Let suppose that the water surface relief is described by the equation $y=r(t,x,y)$ - Figure 1. Figure 1. Flow geometry on the air-water interface with waves of small (left) and large (right) amplitude

We set the boundary conditions for the flow parameters on the top and bottom part of the boundary layer and periodic boundary conditions by coordinate $x$ as follows:

$$ y=0: \mathbf{u}=0,\phi=\rho_w\\ y=H: \mathbf{u}=(U_0,0), \phi=\rho_a\\ \mathbf{u}(t,0,y)=\mathbf{u}(t,L,y)\\ \phi(t,0,y)=\phi(t,L,y). $$ Here $H$ is the height of the boundary layer,$U_0$ is the wind velocity, $L$ is the period of the wave tray. We assume that at the initial time the flow velocity and interface function are given by

$$ t=0, y<H/2:u_x = 0, u_y=V_0 \sin(2 \pi n x/L), \phi =\rho_w\\ t=0, y\ge H/2: u_x=U_0, u_y=V_0\sin (2 \pi n x/L), \phi =\rho_a\\ $$ This problem can be solved with using numerical methods and an appropriate turbulence model. We have used stable fluids algorithm (Stam1999, Stam 2000) to solve 2D Navier-Stokes equations and to simulate nonlinear waves interaction on the air-water interface with a given wind velocity on the top of the boundary layer. While in the standard gravity wave theory the air-water interface is considered mainly as a potential flow, in our approach we take into account velocity, viscosity and density gradients on the interface. Since air/water density ratio is about $10^{-3}$ we have some challenging numerical problem. To solve this problem we made some simplifications in the basic Navier-Stokes equations. First, we suppose that velocity is a continues field at $t>0$ so that flow in the water and flow in the air is an united flow continuously distributed from the bottom to the top and the periodic one in the $x$ direction. We neglect by the surface tension due to large scale of the waves. Therefore we don't need any boundary condition on the air-water interface. Second, at $t>0$ we introduce continues density as $\rho=\phi$ to compute velocity field on every step. Third, we use stable fluid algorithm in time in the very specific order described below.

Stable fluids algorithm

1) Solve advection equation with using boundary conditions, initial data from the previous step, and an implicit algorithm

$$ \frac{\partial \mathbf{u}_1}{\partial t}+(\mathbf{u}_1.\nabla) \bf{u}_1=0 $$

2) Solve diffusion equation with $\nu_i=\nu_i(\phi)$, initial data $\mathbf{u}_1$, boundary conditions and with using, for example, Gauss-Seidel relaxation algorithm

$$ \frac{\partial \mathbf {u}_2}{\partial t}-\nu_i \nabla^2 \mathbf{u}_2=0$$

here $$\nu_i(\phi)=\nu_w \phi ^k, k=-0.4029$$ for the air temperature of 20C.

3) Add force to the velocity field from the previous step as follows

$$ \mathbf{u}3=\mathbf{u}2+\mathbf{f} dt\ $$

To simulate force acting on the volume of air-water mixture we used approximation $$ f_y=-\frac{(\rho_w-\phi)(\phi-\rho_w/2)(\phi-\rho_a)}{\rho_w^2(\rho_w-\rho_a)Fr^2} $$ 4) Make projection step. Here we can use two models. First model is standard projection (Stam 1999) by solving Poison equation $$ t>0, \nabla^2 q=\nabla . \bf{u}_3\\ \bf{u}_4=\bf{u}_3-\nabla q $$ Note, that this step allows us to define divergent free velocity field.

5) Update velocity field $\bf{u}_4\longrightarrow \bf{u}_5$ with using boundary conditions.

6) Make diffusion step with interface function by solving diffusion equation with initial condition from the previous step, with boundary conditions, and with $\nu=\nu_i(\phi)/Sc$ as follows
$$ \frac{\partial \phi_1}{\partial t}-\nu \nabla^2 \phi_1=0 $$ 7) Make advection step by solving advection equation with initial data from step 6 with boundary conditions using an implicit algorithm (from step 1, for example),
$$ \frac{\partial \phi_2}{\partial t}+(\bf{u}_5.\nabla) \phi_2=0 $$ 8) Update interface function with boundary conditions.

9) Return to step 1 with updated velocity and interface function.

This algorithm can be also used in some different order, for example, in the beginning we can compute step 3 (Stam1999). Also we can make steps 6-8 first, and then compute velocity field with steps 1-5. The question about stability of the stable fluids algorithm not really solved yet. From our experience, we can't do arbitrary time step, but $dt$ is limiting by the grid size as usual for more precise computation. In our computations we have used $dt=2/N$ for $N\times N$ grid. In general case the numerical solution depends on the Froude number $Fr=\frac{U_0}{\sqrt{g L}}$, Reynolds number $Re=\frac{U_0 L}{\nu_w}$, Schmidt number and number of waves in the initial data. In our computations we fixed the Froude number, so that we can put $H=L=1$. We also fixed the Reynolds number, therefore we can put $U_0=1.5, V_0=0.5$. The algorithm 1-9 has been implemented with FDM and compiled to C with using Mathematica 12.3.

Code to simulate air/water interface with $Re=10^4, Fr=1$

rhoWater20C = 1.027; nuW20C = 0.01007; rhoAir20C = 0.001204; nuA20C = 0.151;dif = 1/10000; pec = .1; U0 = 1.5; V0 = .5; dn0 = 0.997658; dn1 = 0.514102; kap = 1; n = 81; Fr = 1; F0 = 1; n1 = n + 1; sm = 600; r = 20; den = 
 ConstantArray[dn1 (1 + dn0 Tanh[-kap Range[-n1/2, n1/2]]), n1];u0 = ConstantArray[0, {n1, n1}]; Do[
  u0[[i, j]] = U0 (1 + Tanh[kap (j - n1/2)])/2;, {i, n1}, {j, n1}];
 v0 = ConstantArray[0., {n1, n1}]; Do[
  v0[[i, j]] = V0 Sin[10 Pi (i - 1)/n];, {i, n1}, {j, n1}];periodic[n_, up_, ud_, ub_] := 
    Module[{bd = ub}, Do[bd[[1, i]] = .5 (bd[[n, i]] + bd[[2, i]]);
      bd[[n + 1, i]] = bd[[1, i]]; bd[[i, 1]] = ud;
      bd[[i, n + 1]] = up;, {i, 2, n}];
     bd[[1, 1]] = .5 (bd[[2, 1]] + bd[[1, 2]]);
     bd[[n + 1, n + 1]] = .5 (bd[[n, n + 1]] + bd[[n + 1, n]]);
     bd[[n + 1, 1]] = .5 (bd[[n, 1]] + bd[[n + 1, 2]]);
     bd[[1, n + 1]] = .5 (bd[[1, n]] + bd[[2, n + 1]]); bd];

  diffuse[n_, r_, a_, c_, c0_] := 
    Module[{c1 = c}, 
     Do[Do[Do[
         c1[[i, j]] = (c0[[i, 
               j]] + (a/den[[i, j]]^.4029) (c1[[i - 1, j]] + 
                 c1[[i + 1, j]] + c1[[i, j - 1]] + 
                 c1[[i, j + 1]]))/(1 + 4 a/den[[i, j]]^.4029);, {j, 2, 
          n}];, {i, 2, n}];
      Do[c1[[1, i]] = c1[[n, i]]; c1[[n + 1, i]] = c1[[2, i]];
       c1[[i, 1]] = c0[[i, 1]];
       c1[[i, n + 1]] = c0[[i, n + 1]];, {i, 2, n}];
      c1[[1, 1]] = .5 (c1[[2, 1]] + c1[[1, 2]]);
      c1[[n + 1, n + 1]] = .5 (c1[[n, n + 1]] + c1[[n + 1, n]]);
      c1[[n + 1, 1]] = .5 (c1[[n, 1]] + c1[[n + 1, 2]]);
      c1[[1, n + 1]] = .5 (c1[[1, n]] + c1[[2, n + 1]]);, {k, 0, r}];
     c1];

  advect[n_, d0_, u_, v_, dt_] := 
    Module[{x, y, d1, dt0, i0, i1, j0, j1, s0, s1, t0, t1}, 
     d1 = ConstantArray[0, {n + 1, n + 1}]; dt0 = dt n;
     Do[Do[x = i - dt0 u[[i, j]]; y = j - dt0 v[[i, j]];
        i0 = Which[x <= 1, 1, 1 < x < n, Floor[x], True, n];
        i1 = i0 + 1;
        j0 = Which[y <= 1, 1, 1 < y < n, Floor[y], True, n];
        j1 = j0 + 1; s1 = x - i0; s0 = 1 - s1; t1 = y - j0; t0 = 1 - t1;
        d1[[i, j]] = 
         s0 (t0 d0[[i0, j0]] + t1 d0[[i0, j1]]) + 
          s1 (t0 d0[[i1, j0]] + t1 d0[[i1, j1]]);, {j, 1, n + 1}];, {i, 
       1, n + 1}]; d1];

  project[n_, r_, u0_, v0_, u_, v_] := 
    Module[{ux = u, vy = v, div, p}, 
     p = ConstantArray[0, {n + 1, n + 1}];
     div = ConstantArray[0, {n + 1, n + 1}];
     ux = ConstantArray[0, {n + 1, n + 1}];
     vy = ConstantArray[0, {n + 1, n + 1}];
     Do[div[[i, 
         j]] = -.5/
          n (u0[[i + 1, j]] - u0[[i - 1, j]] + v0[[i, 1 + j]] - 
           v0[[i, j - 1]]);, {i, 2, n}, {j, 2, n}];
     Do[Do[Do[
         p[[i, j]] = (div[[i, 
               j]] + (p[[i - 1, j]] + p[[i + 1, j]] + p[[i, j - 1]] + 
                p[[i, j + 1]]))/4;, {j, 2, n}], {i, 2, n}];, {k, 0, r}];
     Do[ux[[i, j]] = u0[[i, j]] - .5 n (p[[i + 1, j]] - p[[i - 1, j]]);
      vy[[i, j]] = 
       v0[[i, j]] - .5 n (p[[i, j + 1]] - p[[i, j - 1]]);, {i, 2, 
       n}, {j, 2, n}]; {ux, vy}];

  Fx[t_, x_, y_] := 0;
  Fy[t_, x_, y_] := -1/Fr^2;

  onestep[n_, step_, r_, a_, uin_, vin_, dt_, c_] := 
   Module[{u1, v1, f1, f2, u, v, u0, v0}, 
    f1 = ConstantArray[0, {n + 1, n + 1}];
    f2 = ConstantArray[0., {n + 1, n + 1}];
    u0 = ConstantArray[0., {n + 1, n + 1}];
    v0 = ConstantArray[0., {n + 1, n + 1}];
    u = ConstantArray[0., {n + 1, n + 1}];
    v = ConstantArray[0., {n + 1, n + 1}];
    u1 = ConstantArray[0., {n + 1, n + 1}];
    v1 = ConstantArray[0., {n + 1, n + 1}]; u0 = uin; v0 = vin; Do[
     f2[[i, j]] = 
       1/Fr^2 (den[[i, j]] - rhoWater20C/2) (rhoWater20C - 
          den[[i, j]]) (den[[i, j]] - rhoAir20C)/(rhoWater20C - 
            rhoAir20C)/rhoWater20C^2;, {i, 2, n}, {j, 2, n}];
    v0 += f2 dt;
    u0 = advect[n, u0, u0, v0, dt]; v0 = advect[n, v0, u0, v0, dt]; 
    mnV = 0;
    u0 = periodic[n, U0, 0, u0]; v0 = periodic[n, mnV, mnV, v0];
    u0 = diffuse[n, r, a, c, u0]; v0 = diffuse[n, r, a, c, v0];
    u0 = periodic[n, U0, 0, u0]; v0 = periodic[n, mnV, mnV, v0];
    {u1, v1} = project[n, r, u0, v0, u, v];
    u0 = periodic[n, U0, 0, u1]; 
    v0 = periodic[n, mnV, mnV, v1]; {u0, v0}]

  cf = With[{cg = Compile`GetElement, hp = HoldPattern, 
      dv = DownValues}, 
     Hold@Compile[{{u0argu, _Real, 2}, {v0argu, _Real, 
               2}, {denargu, _Real, 
               2}, {sm, _Integer}, {n, _Integer}, {r, _Integer}, dif, 
              pec, F0}, 
             Module[{u0 = u0argu, v0 = v0argu, uu, vv, dd, 
               den = denargu, c = Table[0., {n + 1}, {n + 1}], 
               dt = 40./n^2, a, dnup = den[[1, n + 1]], 
               dnd = den[[1, 1]]}, a = dt dif n n;
              uu = vv = dd = Table[0., {sm + 1}, {n + 1}, {n + 1}];
              Do[

               den = advect[n, den, u0, v0, dt];
               den = periodic[n, dnup, dnd, den]; 
               den = diffuse[n, r, a/pec, c, den];
               den = periodic[n, dnup, dnd, den];

               dd[[step + 1]] = den; {u0, v0} = 
                onestep[n, step, r, a, u0, v0, dt, c];
               uu[[step + 1]] = u0;
               vv[[step + 1]] = v0;, {step, 0, sm}]; {uu, vv, dd}], 
             CompilationTarget -> C, RuntimeOptions -> "Speed"] /. 
           dv@onestep /. 
          Flatten[dv /@ {Fx, Fy, advect, diffuse, periodic, project}] /. 
         hp@ConstantArray[c_, {i_, j_}] :> Table[0., {i}, {j}] /. 
        hp@Part[a__] :> cg[a] /. hp[cg[a__] = rhs_] :> (Part[a] = rhs) //
       ReleaseHold];

Visualization

rst = cf[u0, v0, den, sm, n, r, dif, pec, F0];Do[lstu[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, rst[[1, k, i, j]]}, {i, 
      n1}, {j, n1}], 1]; 
  lstv[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, rst[[2, k, i, j]]}, {i, 
      n1}, {j, n1}], 1];, {k, sm}];
Do[Uvel[i] = Interpolation[lstu[i], InterpolationOrder -> 3];, {i, 1, 
   sm}];
Do[Vvel[i] = Interpolation[lstv[i], InterpolationOrder -> 3];, {i, 1, 
  sm}]; Do[lst4[k] = 
   Flatten[Table[{{(i - 1)/n, (j - 1)/n}, rst[[3, k, i, j]]}, {i, 
      n1}, {j, n1}], 1];, {k, sm}];
Do[rh[k] = Interpolation[lst4[k], InterpolationOrder -> 3];, {k, sm}];{ContourPlot[Uvel[10][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 15, ContourStyle -> Yellow, PlotLegends -> Automatic], 
 ContourPlot[Vvel[10][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 15, ContourStyle -> Yellow, PlotLegends -> Automatic], 
 ContourPlot[1 - rh[10][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 15, ContourStyle -> Yellow, PlotLegends -> Automatic], 
 ContourPlot[Uvel[sm][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 15, ContourStyle -> Yellow, PlotLegends -> Automatic], 
 ContourPlot[Vvel[sm][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 15, ContourStyle -> Yellow, PlotLegends -> Automatic], 
 ContourPlot[1 - rh[sm][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 8, ContourStyle -> Yellow, PlotLegends -> Automatic]}

In this Figure are shown velocity components $u_x$ (left), $u_y$ (middle) and 1- density (right) on step 10 (upper line) and on final step $sm=600$ Figure 2

UPDATE

Next step is to animate wave transformation using interface function. We can show velocity field over density as follows

Show[ContourPlot[1 - rh[sm][x, y]/rhoWater20C, {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "BlueGreenYellow", Contours -> 8,
   ContourStyle -> Yellow, Frame -> False, PlotLegends -> Automatic], 
 StreamPlot[{Uvel[sm][x, y], Vvel[sm][x, y]}, {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, StreamColorFunction -> None, VectorPoints -> Fine,
   VectorColorFunction -> Hue, PlotLegends -> Automatic, 
  StreamColorFunctionScaling -> True, StreamStyle -> LightGray]]

Figure 3

We can also compute frames for animation with using interface function

frames=Table[ContourPlot[1 - rh[i][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 8, ContourStyle -> Yellow, PlotLabel -> i], {i, 20, sm, 
  20}]; Animate[frames]

Figure 4

References

Jos Stam. Stable fluids. In Computer Graphics Proceedings Annual Conference Series, Los Angeles, Aug. 3–8, 199.

Jos Stam. Interacting with smoke and fire in real time. Communications of the ACM, 43(7):77–83, July 2000.

6 Replies

@Alexander, first of all, this is a wonderful work, thank you very much for sharing! I gave a shout-out on LinkedIN and the feedback was great:

https://www.linkedin.com/feed/update/urn:li:activity:6866071848321368064

I am curious about stable fluid algorithm FDM you mentioned. Is there a general reference to it or it was developed in your report? I know NDSolve does use FDM

but I assume your FDM is a very custom thing.

POSTED BY: Vitaliy Kaurov

Dear Vitaliy, thank you very much for links to NDSolve tutorial and LinkedIn. Actually I have used NDSolve with MOL, DAE and FEM solver as well in different applications. Stably fluids algorithm is very specific for implementation with Mathematica. Every step 1-9 discussed above can be realized with NDSolve, but we prefer Module, Do and Compile with option CompilationTarget -> C to speedup simulation. First time this algorithm has been developed by Jos Stam in the last century for computer graphics applications. For our problem we rebuild this algorithm to make it more stable.

I note that the code has multiple nest Do loops, and does not use NDSolve at all. I expect that there is a better way...

Also, the article does not reference Experimental Evidence of Hydrodynamic Instantons: The Universal Route to Rogue Waves, which is definitely relevant to this work.

POSTED BY: Paul Abbott

Thank you very much for reference link. Note, that in our code we have to implement the stable fluid algorithm, therefore we use FDM, not FEM. As nice example of this implementation we can visualize flow around square cylinder
enter image description here

Next step is to animate wave transformation using interface function. We can show velocity field over density as follows

Show[ContourPlot[1 - rh[sm][x, y]/rhoWater20C, {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "BlueGreenYellow", Contours -> 8,
   ContourStyle -> Yellow, Frame -> False, PlotLegends -> Automatic], 
 StreamPlot[{Uvel[sm][x, y], Vvel[sm][x, y]}, {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, StreamColorFunction -> None, VectorPoints -> Fine,
   VectorColorFunction -> Hue, PlotLegends -> Automatic, 
  StreamColorFunctionScaling -> True, StreamStyle -> LightGray]]

Figure 3

We can also compute frames for animation with using interface function

frames=Table[ContourPlot[1 - rh[i][x, y], {x, 0, 1}, {y, 0, 1}, 
  ColorFunction -> "BlueGreenYellow", Frame -> False, 
  PlotRange -> All, ImageSize -> Small, MaxRecursion -> 2, 
  Contours -> 8, ContourStyle -> Yellow, PlotLabel -> i], {i, 20, sm, 
  20}]; Animate[frames]

enter image description here -- you have earned Featured Contributor Badge enter image description here Your exceptional post has been selected for our editorial column Staff Picks http://wolfr.am/StaffPicks and Your Profile is now distinguished by a Featured Contributor Badge and is displayed on the Featured Contributor Board. Thank you!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract