Wind_around_sphere.edp

load "iovtk";

real degree=5.0;//AOA
real radian=-2*pi/360*degree;

border bottom(t=0,10){ x=t;	y=0; 	label=1;}
border out(t=0,5){ x=10;	y=t; 	label=2;}
border top(t=10,0){ x=t;	y=5; 	label=3;}
border in(t=5,0){ x=0;	y=t; 	label=4;}
border circle(t=0,-2*pi){ x=5+cos(t);    y=2.5+sin(t); 	label=5;}
 

int n=1;
mesh Th= buildmesh( bottom(30*n)+out(30*n)+top(30*n)+in(30*n)+circle(30*n));


fespace Vh2(Th,P2);
fespace Vh1(Th,P1);

real beta = 0.1;// beta represents "1/Re"
// Reynolds number:1000, Length; spanwise:1, Velocity: max. vel. on inlet

real dt = 0.005; 
real alpha=1/dt;
int i,niter=10000;
real[int] Drag(niter), Lift(niter); 

Vh2 u1,u2;
Vh1 p;
Vh2 v1,v2;
Vh1 q;
Vh2 up1,up2;

problem  NS ([u1,u2,p],[v1,v2,q],solver=Crout,init=i) =
int2d(Th)(
       alpha*( u1*v1 + u2*v2) 
      + beta * ( dx(u1)*dx(v1) + dy(u1)*dy(v1) +  dx(u2)*dx(v2) + dy(u2)*dy(v2) )
      + p*q*(0.000001) //全周Dirichletなら必要
      - p*dx(v1)- p*dy(v2)
      - dx(u1)*q- dy(u2)*q
     )
+ int2d(Th) ( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 )
+ on(4,u1=-y*(y-5),u2=0) 
+ on(1,3,5,u1=0,u2=0);

for (i=0;i<niter;i++)
{
	
	up1=u1;
	up2=u2;
	
	NS;
	plot(coef=0.2,[u1,u2]);  
	
	Drag[i]=int1d(Th,5)(p*N.x);
	Lift[i]=int1d(Th,5)(p*N.y); 
	
	//plot(coef=0.2,p,[u1,u2], bb=[[0.0,-0.5],[1.2,0.5]]);  
 }
 {
 	real Re=1/beta;
	ofstream ff("Lift-Drag Re="+Re+" AOA="+degree+".csv");
	
	ff << "step" <<","<< "time" <<","<<"Drag"<<","<<"Lift"<<","<<"LiftDrag"<<endl;
	for(int ii1 = 0; ii1 < niter; ii1++)
	{
		real t=ii1*dt;
		real LiftDrag=Lift[ii1]/Drag[ii1];
		
		ff << ii1 <<","<< t <<","<<Drag[ii1]<<","<<Lift[ii1]<<","<<LiftDrag <<endl;
	}
}