Rotational boundary condition

General Discussion
1

Dear community,

I would like to learn how to set an angular velocity as a boundary condition. Hence, I made a simple test case with a clasic problem of a flow around a cylinder and I try to set an angular velocity to the cylinder wall so it will rotate around its own axis with an angular velocity of “1”.

I wrote the boundary condition as: +on(0, u=1cos(2pi), v=1sin(2pi));

I post processed the results in paraview, but I do not understand how the boundary condition works. As you can see in the attached photo there is a velocity field equal to 1 around the cylinder but the velocity glypsh do not show any rotating movement so I am a bit confused.

Any help or suggestion is welcome.

Regards,
Robert

velocity
velocity1612×786 254 KB

2

If the center of the cylinder is at (x,y)=(0,0), the angular velocity is S, and the radius is R, then you can use:

u=-y/sqrt(x^2+y^2)*S*R;
v=x/sqrt(x^2+y^2)*S*R;
3

Can you post the code? I think there is something like Tl or similar for tangent
line similar to normal. That gives you the direction of the boundary at every
point which you can relate to velocity or whatever vector you need. Search
the manual for “Tl” or tangent

4

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

5

If you are trying to create a boundary layer with zero relative velcoity betwen
the moving surface and the fluid, r probably the best way to do that is
to use the tangent line function Inidicated above. See the changes below.
Also the absoluate value of p does not really appear here AFACIT
and the examples include a small term to help it converge.
( while you do have a “p” term it can be integrated by parts and it goes away lol ).
Grep the examples
for “convect” and “NS”

I verified this compiled and ran but I did not look at the results as
your mesh in this example looks a lot different from the original picture.

 diff forum10*
2c2
< 	
---
> load "medit"	
52c52,53
< 		  +on(0, u=-y/sqrt(x^2+y^2)*1, v=x/sqrt(x^2+y^2)*1)
---
> 		  //+on(0, u=-y/sqrt(x^2+y^2)*1, v=x/sqrt(x^2+y^2)*1)
> 		  +on(0, u=Tl.x, v=Tl.y)
6

For some reason Tl did not work but apparently the normal is ok,

forum10
forum101067×601 242 KB 

diff forum10*
2c2
< 	
---
> load "medit"	
9c9,10
< border circle(t=0, 2*pi){x=R*cos(t); y=R*sin(t); label=0;}
---
> //border circle(t=0, 2*pi){x=R*cos(t); y=R*sin(t); label=0;}
> border circle(t=0, 2*pi){x=R*cos(t); y=R*sin(t); label=99;}
34c35
< real dt = 1;             
---
> real dt =1;             
40c41,42
< problem NS([u,v,p],[hx,hy,q])
---
> //problem NS([u,v,p],[hx,hy,q],solver=Crout,tgv=-1)
> problem NS([u,v,p],[hx,hy,q],solver=CG,tgv=-1)
42,44c44,47
< 		   (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)) 	   
---
> 		   (1.0/dt)*(u*hx+v*hy)				       
> 		  +(1.0/Re)*(dx(u)*dx(hx)+dy(u)*dy(hx)) 	   
> 		  +(1.0/Re)*(dx(v)*dx(hy)+dy(v)*dy(hy)) 	   
> 			+p*q*1e-10
49,50c52,53
<           - (1/dt)*convect([up,vp],-dt,up)*hx
<           - (1/dt)*convect([up,vp],-dt,vp)*hy
---
>           - (1.0/dt)*convect([up,vp],-dt,up)*hx
>           - (1.0/dt)*convect([up,vp],-dt,vp)*hy
52c55,58
< 		  +on(0, u=-y/sqrt(x^2+y^2)*1, v=x/sqrt(x^2+y^2)*1)
---
> 		 // +on(0, u=-y/sqrt(x^2+y^2)*1, v=x/sqrt(x^2+y^2)*1)
> 		  //+on(99, u=Tl.x, v=Tl.y)
> 		  +on(99, u=-N.y, v=N.x)
> 		  //+on(99, u=1000, v=1000)
60c66,67
<       up[] = u[];
---
> //      up[] = u[];
> 	[up,vp,pp]=[u,v,p];
64c71,73
< 		savevtk("Axial"+i+".vtu", Th,[u,v,0], dataname = DataName, order = Order);
---
> //		savevtk("Axial"+i+".vtu", Th,[u,v,0], dataname = DataName, order = Order);
> //	medit(""+i,Th,p);
> 	plot([u,v],wait=1,fill=1);
1 Like
7

There are a lot of example cases in the Stabfem library (StabFem · GitLab), you can find many solutions and useful tricks there.

1 Like