Navier-Stokes solver documentation¶
Here is a verbatim of the Lorena Barba 12 steps to Naviers Stokes computation, an excellent step-by-step course oriented on scientific computing. The solver is incompressible and deals with the Pressure Poisson equation using a martix-free, pseudo time integration.
Governing equations:¶
Here are our Navier-Stokes equations:
du du du - 1 dp ( d2 u d2 u )
-- + u -- + v -- = --- -- + nu ( ---- + ---- )
dt dx dy rho dx ( dx2 dy2 )
dv dv dv - 1 dp ( d2 v d2 v )
-- + u -- + v -- = --- -- + nu ( ---- + ---- )
dt dx dy rho dy ( dx2 dy2 )
d2 p d2 p ( du du du dv dv dv )
---- + ---- = - rho ( -- -- + 2 -- -- + -- -- )
dx2 dy2 ( dx dx dy dx dy dy )
With patience and care, we write the discretized form of the equations.
The u-momentum equation:
u(n+1,i,j) - u(n,i,j)
---------------------
dt
u(n,i,j) - u(n,i-1,j) u(n,i,j) - u(n,i,j-1)
+ u(n,i,j) --------------------- + v(n,i,j) ---------------------
dx dy
-1 p(n,i+1,j) - p(n,i-1,j)
= --- -----------------------
rho 2 dx
( u(n,i+1,j) - 2 u(n,i,j) + u(n,i-1,j) u(n,i,j+1) - 2 u(n,i,j) + u(n,i,j-1) )
+ nu ( ------------------------------------ + ------------------------------------ )
( dx^2 dy^2 )
+ f(i,j)
The v-momentum equation:
v(n+1,i,j) - v(n,i,j)
---------------------
dt
v(n,i,j) - v(n,i-1,j) v(n,i,j) - v(n,i,j-1)
+ u(n,i,j) --------------------- + v(n,i,j) ---------------------
dx dy
-1 p(n,i+1,j) - p(n,i-1,j)
= --- -----------------------
rho 2 dy
( v(n,i+1,j) - 2 v(n,i,j) + v(n,i-1,j) v(n,i,j+1) - 2 v(n,i,j) + v(n,i,j-1) )
+ nu ( ------------------------------------ + ------------------------------------ )
( dx^2 dy^2 )
And the pressure equation:
p(n,i+1,j) - 2 p(n,i,j) + p(n,i-1,j) p(n,i,j+1) - 2 p(n,i,j) + p(n,i,j-1)
------------------------------------ + ------------------------------------
dx^2 dy^2
( 1 ( u(n,i+1,j) - u(n,i-1,j) v(n,i,j+1) - v(n,i,j-1) )
= rho ( -- ------------------------- + ----------------------- )
( dt 2 dx 2 dy )
u(n,i+1,j) - u(n,i-1,j) u(n,i+1,j) - u(n,i-1,j)
- ----------------------- -----------------------
2 dx 2 dx
u(n,i,j+1) - u(n,i,j-1) v(n,i+1,j) - v(n,i-1,j)
- 2 ----------------------- -----------------------
2 dy 2 dx
v(n,i,j+1) - v(n,i,j-1) v(n,i,j+1) - v(n,i,j-1)
- ----------------------- -----------------------
2 dy 2 dy
Resolution of momentum:¶
As always, we re-arrange these equations to the form we need in the code to make the iterations proceed.
The momentum equation in the u direction:
u(n+1,i,j) = u(n,i,j)
- u(n,i,j) dt/dx ( u(n,i,j) - u(n,i-1,j)
- v(n,i,j) dt/dy ( u(n,i,j) - u(n,i,j-1)
- dt / ( 2 rho dx ) ( p(n,i+1,j) - p(n,i-1,j) )
+ nu dt/dx^2 ( u(n,i+1,j) - 2 u(n,i,j) + u(n,i-1,j) )
+ nu dt/dy^2 ( u(n,i,j+1) - 2 u(n,i,j) + u(n,i,j-1) )
+ dt * f(i,j)
The momentum equation in the v direction:
v(n+1,i,j) = v(n,i,j)
- u(n,i,j) dt/dx ( v(n,i,j) - v(n,i-1,j)
- v(n,i,j) dt/dy ( v(n,i,j) - v(n,i,j-1)
- dt / ( 2 rho dy ) ( p(n,i,j+1) - p(n,i,j-1) )
+ nu dt/dx^2 ( v(n,i+1,j) - 2 v(n,i,j) + v(n,i-1,j) )
+ nu dt/dy^2 ( v(n,i,j+1) - 2 v(n,i,j) + v(n,i,j-1) )
Resolution of Pressure:¶
And for the pressure equation, we isolate the term p(n,i,j) to iterate in pseudo-time:
( p(n,i+1,j) + p(n,i-1,j) ) dy^2 + ( p(n,i,j+1) + p(n,i,j-1) ) dx^2
p(n,i,j) = -------------------------------------------------------------------
2 ( dx^2 + dy^2 )
rho dx^2 dy^2
- -----------------
2 ( dx^2 + dy^2 )
[ 1 ( u(n,i+1,j) - u(n,i-1,j) v(n,i,j+1) - v(n,i,j-1) )
* [ -- ( ----------------------- + ----------------------- )
[ dt ( 2 dx 2 dy )
- u(n,i+1,j) - u(n,i-1,j) u(n,i+1,j) - u(n,i-1,j)
----------------------- -----------------------
2 dx 2 dx
u(n,i,j+1) - u(n,i,j-1) v(n,i+1,j) - v(n,i-1,j)
- 2 ----------------------- -----------------------
2 dy 2 dx
v(n,i,j+1) - v(n,i,j-1) v(n,i,j+1) - v(n,i,j-1) ]
- ----------------------- ----------------------- ]
2 dy 2 dy ]
Present BARBATRUC implementation¶
The present implementation have the following differences with the final step of 12-steps-to-NS.
- The spatial resolution
dxanddyare equaland constant, because we wanted to share the same data-structure as a 2DQ9 LBM solver. This simplifies equations a lot. - We added a scalar transport equation.
- We added a “crude” immersed boundary method.
- There is a convergence criterion in the Pressure solver.
In terms of implementation, the most obvious difference is the separation of derivatives in a dedicated routine. It makes the code more readable, by removing the Finite-Difference clutter from the solver, and allows testing.