Unsteady

Under the same asumptions of the steady case, let \(T\) be the study time, the unsteady Stokes problem reads:

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

Where \(\rho\) is the density.

Variational form

The variational form is calculated similarly to the steady Stokes case:

\[\begin{eqnarray*} \rho\int_{\Omega}{\frac{\partial\mathbf{u}_h}{\partial t}\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*}\]

Time discretization

For the time discretization, a finite difference in time is used:

\[\frac{\partial\mathbf{u}_h}{\partial t}\approx\frac{\mathbf{u}_h^{n+1} - \mathbf{u}_h^n}{\Delta t}\]

Where \(\Delta t\) is the time step. It holds the following discretized variational form:

\[\begin{eqnarray*} \frac{\rho}{\Delta t}\int_{\Omega}{(\mathbf{u}_h^{n+1}-\mathbf{u}_h^n)\cdot\mathbf{v}_h} + \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

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

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
problem Stokes ([u1,u2,u3,p], [v1,v2,v3,q], solver=UMFPACK)
    = 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) * (u1old*v1 + u2old*v2 + u3old*v3)
        + f1*v1 + f2*v2 + f3*v3
    )
    + on(1, 2, 3, 4, 5, 6, u1=uEx, u2=vEx, u3=wEx)
    ;

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

    //Solve
    Stokes;

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

Uzawa method

The method is already described in the steady Stokes part.

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

For more information, see toto 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);
}