Unsteady
Under the same asumptions of the steady case, let T be the study time, the unsteady Stokes problem reads:
ρ∂u∂t−μΔu+∇p=f in [0,T]×Ωdiv(u)=0 in [0,T]×Ωu=g1 on [0,T]×ΓDμ∂u∂n=g2 on [0,T]×ΓN
Where ρ is the density.
Variational form
The variational form is calculated similarly to the steady Stokes case:
ρ∫Ω∂uh∂t⋅v+μ∫Ω∇uh⋅∇vh−∫Ωphdiv(vh)−∫ΓNg2⋅vh=∫Ωf⋅vh∫Ωdiv(uh)qh=0
Time discretization
For the time discretization, a finite difference in time is used:
∂uh∂t≈un+1h−unhΔt
Where Δt is the time step. It holds the following discretized variational form:
ρΔt∫Ω(un+1h−unh)⋅vh+μ∫Ω∇un+1h⋅∇vh−∫Ωpn+1hdiv(vh)−∫ΓNgn+12⋅vh=∫Ωfn+1⋅vh∫Ω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