Good evening,
I’m writing a code to calculate the unsteady flow around circular cylinder. As a time discretization I’m using backward Euler.
To do so i calculated the following weak form:
F(\mathbf{u}^{n+1})=\int_\Omega \frac{1}{\Delta t}\mathbf{u}^{n+1}\cdot\mathbf{v}\,d\Omega+ \int_\Omega\nu\nabla\mathbf{u}^{n+1}:\nabla\mathbf{v}\,d\Omega+\int_\Omega\mathbf{u}^{n+1}\cdot\nabla\mathbf{u}^{n+1}\cdot\mathbf{v}\,d\Omega-\int_\Omega p^{n+1}(\nabla\cdot\mathbf{v})\,d\Omega-\int_\Omega(\nabla\cdot\mathbf{u}^{n+1}
)\mathbf{q}\,d\Omega-\int_\Omega\frac{1}{\Delta t}\mathbf{u}^n\cdot\mathbf{v}\,d\Omega=0
DF(\mathbf{u}^{n+1})(\delta\mathbf{u},\delta p)=\int_\Omega \frac{1}{\Delta t}\mathbf{\delta\mathbf{u}}\cdot\mathbf{v}\,d\Omega+ \int_\Omega\nu\nabla\delta\mathbf{u}:\nabla\mathbf{v}\,d\Omega+\int_\Omega(\delta\mathbf{u}\cdot\nabla\mathbf{u}^{n+1}\cdot\mathbf{v}+\mathbf{u}^{n+1}\cdot\nabla\delta\mathbf{u}\cdot\mathbf{v})\,d\Omega-\int_\Omega \delta p(\nabla\cdot\mathbf{v})\,d\Omega-\int_\Omega(\nabla\cdot\delta\mathbf{u})\mathbf{q}\,d\Omega=0
and I wrote the following code:
real Re=60.;
real nu=1./Re;
real eps=1e-8;
real tol=1e-6;
real L=30.;
real R=15.;
real M=30;
border bd(t=-R,L){x=t;y=-R;label=2;}
border bl(t=R,-R){x=-R;y=t;label=1;}
border bu(t=L,-R){x=t;y=R;label=2;}
border br(t=-R,R){x=L;y=t;label=3;}
border cc(t=0,-2*pi){x=cos(t)/2.;y=sin(t)/2.;label=0;}
mesh Th=buildmesh(bd(M/2)+bl(M/2)+bu(M/2)+br(M/2)+cc(M));
plot(Th,wait=1);
real dt=0.5;
real Niter=100.;
real alfa=1./dt;
fespace Vh(Th,P2);
Vh u1,u2;
Vh v1,v2;
Vh du1,du2;
Vh un,vn;
fespace Wh(Th,P1);
Wh p;
Wh dp;
Wh q;
Wh pn;
macro Grad(u1,u2) [ dx(u1),dy(u1), dx(u2),dy(u2)]//
macro UgradV(u1,u2,v1,v2) [ [u1,u2]'*[dx(v1),dy(v1)] , [u1,u2]'*[dx(v2),dy(v2)] ]//
macro div(u1,u2) (dx(u1)+dy(u2))//
solve Stokes([u1,u2,p],[v1,v2,q],solver=sparsesolver)=
int2d(Th)(
(dx(u1)*dx(v1)
+ dy(u1)*dy(v1)
+ dx(u2)*dx(v2)
+ dy(u2)*dy(v2)
)
- eps*p*q
-p*div(v1,v2)
-div(u1,u2)*q
)
+ on(0,u1=0,u2=0)
+ on(1,u1=1,u2=0)
;
plot(coef=0.2,cmm="[u1,u2], p Stokes",p,[u1,u2],wait=1,value=true);
un[]=u1[];
vn[]=u2[];
pn[]=p[];
for (int i=0;i<Niter;i++){
real err=0;
real Time=i*dt;
for(int j=0;j<=15;j++){
solve NSnewton([du1,du2,dp],[v1,v2,q])=
int2d(Th)(
alfa*(du1*v1+du2*v2)
+nu*(Grad(du1,du2)'*Grad(v1,v2))
+UgradV(du1,du2,u1,u2)'*[v1,v2]
+UgradV(u1,u2,du1,du2)'*[v1,v2]
-div(v1,v2)*dp
-div(du1,du2)*q
-eps*dp*q
)
+int2d(Th)(
alfa*(u1*v1+u2*v2)
+nu*(Grad(u1,u2)'*Grad(v1,v2))
+UgradV(u1,u2,u1,u2)'*[v1,v2]
-div(v1,v2)*p
-div(u1,u2)*q
)
-int2d(Th)(
alfa*(un*v1+vn*v2)
)
+ on(0,1,du1=0.0,du2=0.0)
;
u1[]=u1[]+du1[];
u2[]=u2[]+du2[];
p[]=p[]+dp[];
real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
err=du1[].linfty/Lu1+du2[].linfty/Lu2+dp[].linfty/Lp;
cout<<"iteration number="<<j+1<<"error="<<err<<endl;
if(err<tol)break;
if(err>10 && j>3)break;
}
un[]=u1[];
vn[]=u2[];
pn[]=p[];
real rey=1./nu;
plot(coef=0.2,cmm="Time="+Time+"s,Re="+rey+"[u1,u2] e p",p,[un,vn],wait=0,value=true);
}
I have achieved quadratic convergence but there is a significant transient from the stoke solution to the true solution (there is region with low velocity behind the cylinder wich dissappears only after +100 iterations) and it does’t reproduce vortex shedding for Re>50.
I also want to adapt my mesh to the time dependent problem and found the lenght of the recirculation zone for Re<40.
Any suggest?