Flow around circular cylinder with Newton method

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

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)=
       + dy(u1)*dy(v1)
       + dx(u2)*dx(v2)
       + dy(u2)*dy(v2)
       - eps*p*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);


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])=
    + on(0,1,du1=0.0,du2=0.0)
  real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
  cout<<"iteration number="<<j+1<<"error="<<err<<endl;
if(err>10 && j>3)break;
  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?

Update: for the calculation of the length of the recirculation bubble i added the following lines of code so as to verify where u1 changes sign.

{ real x0 = 0.5;
  real dx = 0.01;
  int i;
  ofstream f("punti_Re_M.txt");
  for (i=0;i<1000;i++){
    f << x0+i*dx << "    "<< un(x0+i*dx,0)<< endl;

I’ve done many test with different mesh, from M=20 to M=150, different time step, from dt=10 to dt=0.01, in order to find the correlation between the recirculation length and Reynolds number.

Unfortunately the results do not match with what i found in literature: from Coutanceau and Bouard (1977) (I found this one on “Flow around circular cylinders: volume 1, Zdravkovich”) we have a linear empirical relationtship L_w/D=0.05Re, for 4.4<Re<40, while the one i calculated is in the form L_w/D=0.068Re-0.43.
The results are the same for both M=80 and M=150.

The code still does not reproduce the unsteady regime, but rather goes to convergence, i.e. the error of the first newton iteration becomes smaller than 10^{-7}.

What am I doing wrong?