Navier-Stokes equations steady case

\[\begin{eqnarray} \rho\frac{\partial\mathbf{u}}{\partial t} + \rho(\mathbf{u}\cdot\nabla)\mathbf{u} -\mu\Delta\mathbf{u} + \nabla p &= \textbf{f} &\mbox{ in }[0,T]\times\Omega\label{eqn:navierstokes1}\\ \mathrm{div}(\mathbf{u}) &= 0 &\mbox{ in }[0,T]\times\Omega\label{eqn:navierstokes2}\\ \mathbf{u} &= \textbf{g}_1 &\mbox{ on }[0,T]\times\Gamma_D\label{eqn:navierstokes3}\\ \mu\frac{\partial\mathbf{u}}{\partial\mathbf{n}} &= \textbf{g}_2 &\mbox{ on }[0,T]\times\Gamma_N\label{eqn:navierstokes4} \end{eqnarray}\]

With the operator:

\[\mathbf{u}\cdot\nabla\mathbf{u} = \sum_{i=1}^{N}{\mathbf{u}_i\frac{\partial\mathbf{u}}{\partial x_i}}\]

The Navier-Stokes equations reveal a non-linear term that we can manage by different way.

Variational form

The variational form is similar to the unsteady Stokes one, we have:

\[\begin{eqnarray*} \rho\int_{\Omega}{\frac{\partial\mathbf{u}_h}{\partial t}\cdot\mathbf{v}} + \rho\int_{\Omega}{(\mathbf{u}\cdot\nabla)\mathbf{u}\cdot\mathbf{v}} + \mu\int_{\Omega}{\nabla\mathbf{u}_h\cdot\nabla\mathbf{v}_h} - \int_{\Omega}{p_h\mathrm{div}(\mathbf{v}_h)} - \int_{\Gamma_N}{\textbf{g}_2\cdot\mathbf{v}_h} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}_h}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u}_h)q_h} &=& 0 \end{eqnarray*}\]

To manage the non-linearity, the characteristics method is used. For every particles, we write:

\[\begin{eqnarray*} \frac{d\textbf{X}}{dt}(\textbf{x},s;t) &=& \mathbf{u}(t, \textbf{X}(\textbf{x},s;t))\\ \textbf{X}(\textbf{x},x;x) &=& \textbf{x} \end{eqnarray*}\]

Where \(\textbf{X}(\textbf{x},s;t)\) is the particle position at time t who was in x at time s. We can approximate the non-linear term by:

\[\left(\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right)(t^n,\textbf{x}) = \frac{\mathbf{u}(t^{n+1},\textbf{x})-\mathbf{u}(t^n,\textbf{X}^n(\textbf{x})}{\Delta t}+\mathcal{O}(\Delta t)\]

It holds the following discretized variational form:

\[\begin{eqnarray*} \frac{\rho}{\Delta t}\int_{\Omega}{\mathbf{u}_h^{n+1}\cdot\mathbf{v}_h} - \frac{\rho}{\Delta t}\int_{\Omega}{(\mathbf{u}_h^n\circ\textbf{X}_h^n)\cdot\mathbf{v}_h}&&\nonumber\\ + \mu\int_{\Omega}{\nabla\mathbf{u}_h^{n+1}\cdot\nabla\mathbf{v}_h} - \int_{\Omega}{p_h^{n+1}\mathrm{div}(\mathbf{v}_h)} - \int_{\Gamma_N}{\textbf{g}_2^{n+1}\cdot\mathbf{v}_h} &=& \int_{\Omega}{\textbf{f}^{n+1}\cdot\mathbf{v}_h}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u}_h^{n+1})q_h} &=& 0 \end{eqnarray*}\]

FreeFem++ code

load "msh3"

//Mesh
int nn = 8;
mesh3 Th = cube(nn, nn, nn);

//Parameters
real rho = 1.;
real mu = 1.;
real T = 0.05;
real dt = 1.e-2;

func f1 = 0.;
func f2 = 0.;
func f3 = 0.;

//Taylor Hood spaces
fespace Vh(Th,P2);
Vh u1, u2, u3;
Vh u1old, u2old, u3old;
Vh v1, v2, v3;
fespace Qh(Th, P1);
Qh p, q;

//Weak formulation
problem NavierStokes ([u1,u2,u3,p], [v1,v2,v3,q], solver=sparsesolver)
    = int3d(Th)(
          (rho/dt) * (u1*v1 + u2*v2 + u3*v3)
        + mu * (
              dx(u1)*dx(v1) + dy(u1)*dy(v1) + dz(u1)*dz(v1)
            + dx(u2)*dx(v2) + dy(u2)*dy(v2) + dz(u2)*dz(v2)
            + dx(u3)*dx(v3) + dy(u3)*dy(v3) + dz(u3)*dz(v3)
        )
        - p * q * 1e-10
        - p*dx(v1) - p*dy(v2) - p*dz(v3)
        - dx(u1)*q - dy(u2)*q - dz(u3)*q
    )
    - int3d(Th)(
          (rho/dt) * (
              convect([u1old, u2old, u3old], -dt, u1old) * v1
            + convect([u1old, u2old, u3old], -dt, u2old) * v2
            + convect([u1old, u2old, u3old], -dt, u3old) * v3
          )
        + f1*v1 + f2*v2 + f3*v3
    )
    + on(6, u1=1., u2=0., u3=0.)
    + on(1, 2, 3, 4, 5, u1=0., u2=0., u3=0.)
    ;

for (int i = 0; i < T/dt; i++){
    //Update
    u1old = u1;
    u2old = u2;
    u3old = u3;

    //Solve
    NavierStokes;

    //Plot
    plot(p, fill=true);
}

Uzawa method

The method is already described in the Stokes steady case part.

The following program solves the unsteady Stokes equations using the Uzawa method on a cube with a given velocity:

For more information, see Uzawa method

load "msh3"

//Mesh
int nn = 8;
mesh T2d = square(nn, nn);
mesh3 Th = buildlayers(T2d, nn, zbound=[0., 1.]);

//Parameters
real rho = 1.;
real mu = 1.;
real T = 0.05;
real dt = 1.e-2;

//Easy ktest
func f1 = 1.;
func f2 = 1.;
func f3 = 1.;
func uEx = x;
func vEx = y;
func wEx = -2*z;
func pEx = x+y+z-3/2.;

//Taylor Hood spaces
fespace Vh(Th,P2);
Vh u1, u2, u3;
Vh u1old, u2old, u3old;
Vh v1, v2, v3;
fespace Qh(Th, P1);
Qh p, q;

//Weak formulation
macro grad(a) [dx(a), dy(a), dz(a)] //

varf vA(u, v)
    = int3d(Th)(
          (rho/dt) * u * v
        + mu * grad(u)' * grad(v)
    )
    + on(1, 2, 3, 4, 5, 6, u=0.)
    ;

varf vBX(q, v)
    = int3d(Th)(
          q * dx(v)
    )
    ;

varf vBY(q, v)
    = int3d(Th)(
          q * dy(v)
    )
    ;

varf vBZ(q, v)
    = int3d(Th)(
          q * dz(v)
    )
    ;

varf vOnX(u, v)
    = on(1, 2, 3, 4, 5, 6, u=uEx)
    + int3d(Th)(
        (rho/dt) * u1*v
    )
    + int3d(Th)(
        f1 * v
    )
    ;

varf vOnY(u, v)
    = on(1, 2, 3, 4, 5, 6, u=vEx)
    + int3d(Th)(
        (rho/dt) * u2*v
    )
    + int3d(Th)(
        f2 * v
    )
    ;

varf vOnZ(u, v)
    = on(1, 2, 3, 4, 5, 6, u=wEx)
    + int3d(Th)(
        (rho/dt) * u3 * v
    )
    + int3d(Th)(
        f3 * v
    )
    ;

matrix A = vA(Vh, Vh, solver=sparsesolver);
matrix BX = vBX(Qh, Vh);
matrix BY = vBY(Qh, Vh);
matrix BZ = vBZ(Qh, Vh);

real[int] OnX = vOnX(0, Vh);
real[int] OnY = vOnY(0, Vh);
real[int] OnZ = vOnZ(0, Vh);

func real[int] Uzawa (real[int] &pp){
    real[int] bX = OnX; bX += BX * pp;
    real[int] bY = OnY; bY += BY * pp;
    real[int] bZ = OnZ; bZ += BZ * pp;

    u1[] = A^-1 * bX;
    u2[] = A^-1 * bY;
    u3[] = A^-1 * bZ;

    pp  = BX' * u1[];
    pp += BY' * u2[];
    pp += BZ' * u3[];
    pp  = -pp;

    return pp;
}

u1 = uEx;
u2 = vEx;
u3 = wEx;
p = pEx;
for (int i = 0; i < T/dt; i++){
    //Update
    OnX = vOnX(0, Vh);
    OnY = vOnY(0, Vh);
    OnZ = vOnZ(0, Vh);

    //Solve
    LinearCG(Uzawa, p[], eps=1.e-6);

    //Plot
    plot(p, fill=true);
}