Navier-Stokes equations unsteady 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
The following program solves the unsteday Navier-Stokes equations on a cube with a given velocity:
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
More on Uzawa method and Cahouet-Chabard preconditionner
As in the Stokes problem, we can use the Uzawa method. A Cahouet-Chabard preconditionner [CC98] is used.
| For more information, see Uzawa method and Cahouet-Chabard preconditionner. |
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;
fespace Qh(Th, P1);
Qh p;
//Weak formulation
macro grad(a) [dx(a), dy(a), dz(a)] //
varf vP(p, q)
= int3d(Th)(
grad(p)' * grad(q)
)
;
varf vM(p, q)
= int3d(Th)(
p * q
)
;
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(6, u=1.)
+ on(1, 2, 3, 4, 5, u=0.)
+ int3d(Th)(
(rho/dt) * convect([u1, u2, u3], -dt, u1)*v
)
+ int3d(Th)(
f1 * v
)
;
varf vOnY(u, v)
= on(6, u=0.)
+ on(1, 2, 3, 4, 5, u=0.)
+ int3d(Th)(
(rho/dt) * convect([u1, u2, u3], -dt, u2)*v
)
+ int3d(Th)(
f2 * v
)
;
varf vOnZ(u, v)
= on(6, u=0.)
+ on(1, 2, 3, 4, 5, u=0.)
+ int3d(Th)(
(rho/dt) * convect([u1, u2, u3], -dt, u3)* v
)
+ int3d(Th)(
f3 * v
)
;
matrix PP = vP(Qh, Qh, solver=sparsesolver);
matrix PM = vM(Qh, Qh, solver=sparsesolver);
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;
}
func real[int] Precon (real[int] & p){
real[int] pp = PP^-1 * p;
real[int] pm = PM^-1 * p;
real[int] ppp = (rho/dt) * pp + mu * pm;
return ppp;
}
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[], precon=Precon, eps=1.e-6);
//Plot
plot(p, fill=true);
}