Unsteady

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

ρutμΔu+p=f in [0,T]×Ωdiv(u)=0 in [0,T]×Ωu=g1 on [0,T]×ΓDμun=g2 on [0,T]×ΓN

Where ρ is the density.

Variational form

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

ρΩuhtv+μΩuhvhΩphdiv(vh)ΓNg2vh=ΩfvhΩdiv(uh)qh=0

Time discretization

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

uhtun+1hunhΔt

Where Δt is the time step. It holds the following discretized variational form:

ρΔtΩ(un+1hunh)vh+μΩun+1hvhΩpn+1hdiv(vh)ΓNgn+12vh=Ωfn+1vhΩdiv(un+1h)qh=0

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);
}
freefem

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);
}
freefeem