Thank you for the solution and ideas!

```
load "iovtk"
//CREATE THE MESH//
real R=0.5;
real l=80;
real L=120;
real t;
border circle(t=0, 2*pi){x=R*cos(t); y=R*sin(t); label=0;}
border inlet(t=120, 40) {x=-40; y=l-t; label=1;}
border outlet(t=120, 40) {x=L-40; y=l-t; label=3;}
border upper(t=0, L) {x=t-40; y=40; label=4;}
border lower(t=0, L) {x=t-40; y=-40; label=4;}
border internalleft(t=-8, 8) {x=-4; y=t;}
border internalright(t=-8, 8) {x=40; y=t;}
border internalup (t=-4, 40) {x=t; y=-8;}
border internaldown(t=-4, 40) {x=t; y=8;}
mesh Th=buildmesh(inlet(-25)+outlet(25)+upper(-25)+lower(25)+circle(-75)+internalleft(-75)+internalright(75)+internalup(-75)+internaldown(75));
fespace Xh(Th,P2);
fespace Mh(Th,P1);
fespace XXXFES(Th,[P2,P2,P1]);
XXXFES[u,v,p];
XXXFES[hx,hy,q];
XXXFES[up,vp,pp];
int[int] Order = [1];
string DataName = "Vector";
real Re = 60 ;
real uinlet = 1 ;
real dt = 1;
real time = 0.;
int iter = 50;
[u,v,p]=[0, 0, 0];
problem NS([u,v,p],[hx,hy,q])
=int2d(Th) (
(1/dt)*(u*hx+v*hy)
+(1/Re)*(dx(u)*dx(hx)+dy(u)*dy(hx))
+(1/Re)*(dx(v)*dx(hy)+dy(v)*dy(hy))
- p*(dx(hx)+dy(hy))
-(dx(u)+dy(v))*q
)
+int2d(Th) (
- (1/dt)*convect([up,vp],-dt,up)*hx
- (1/dt)*convect([up,vp],-dt,vp)*hy
)
+on(0, u=-y/sqrt(x^2+y^2)*1, v=x/sqrt(x^2+y^2)*1)
+on(1, u=0, v=0);
int i;
for (i = 1; i <= iter; i++)
{
time += dt;
cout << "Time step " << i << " - t = " << time << endl;
up[] = u[];
NS;
if( i >= 1 )
{
savevtk("Axial"+i+".vtu", Th,[u,v,0], dataname = DataName, order = Order);
}
}
```

I attach the code. I didn’t write anymore the “Radius” in the velocity equation as it seems that it literally multiplies the angular velocity by the radius value (instead of 1 I get 0.5).

I am not sure if I should create a new topic but I would have one more question, just for the curiosity:

“If the center of the cylinder is at `(x,y)=(0,0)`

”…how can we move the coordinate system if the center is not at (x,y)=(0,0)?

For example if the circle is at (x,y)=(10,5)

```
border circle(t=0, 2*pi){x=10+R*cos(t); y=5+R*sin(t); label=0;}
```

I guess the BC should be

```
u=-ynew/sqrt(xnew^2+ynew^2)*S*R;
v=xnew/sqrt(xnew^2+ynew^2)*S*R;
```

So a new pair of coordinate is needed, something like:

real xnew = 10+x;

real ynew =5+y;

How should be this defined?

Regards,

Robert