Good evening everyone,
this is the code I wrote to study a cylinder that oscillates with an imposed motion in a direction perpendicular to the asymptotic flow. I have developed an ALE method. I put a dimensionless amplitude of oscillations of A = 0.25, a dimensionless frequency F = 1.1 (we should develop the phenomenon of the look in) and a Re = 100. I ask you if you think the code is correct, because the trends in the coefficient of lift and drag that I have seen in literature do not return. Thank you all.
ofstream Drug(“cd.txt”);
ofstream Lift(“cl.txt”);
int n=20, nc=80, M=100;
real Re=100, invRe=1/Re, KC=5. ,U=1, D=1, R=D/2., A=KCD/(2pi), F=1.1;
real L=50D, H=30D;
real K=20, S=10;
real x0=(-1./2./pi)sin(2pi0), y0=0, nu=1./Re, T=20, dt=0.01, a=1/dt, er=1e-8;
border b1(t=-L,L){x=t;y=-H;label=1;}
border b2(t=-H,H){x=L;y=t;label=2;}
border b3(t=L,-L){x=t;y=H;label=1;}
border b4(t=H,-H){x=-L;y=t;label=4;}
border c(t=0,2pi){x=x0+Rcos(t);y=y0+Rsin(t);label=3;}
//border b5(t=-K,K){x=t;y=-S;}
//border b6(t=-S,S){x=K;y=t;}
//border b7(t=K,-K){x=t; y=S;}
//border b8(t=S,-S){x=-K; y=t;}
//border b5(t=0,1){x=-10+50t;y=-5;}
//border b6(t=0,1){x=40;y=-5+10t;}
//border b7(t=0,1){x=40-50t;y=5;}
//border b8(t=0,1){x=-10; y=5-10t;}
//border e(t=0, 2pi){x=2cos(t);y=2sin(t);}
//mesh Th=buildmesh(b1(50)+b2(50)+b3(50)+b4(nH/L)+b5(200)+b6(100)+b7(200)+b8(100)+c(-80));
//mesh Th=buildmesh(b1(80)+b2(80)+b3(80)+b4(80)+e(80)+c(-80));
mesh Th=buildmesh(b1(n)+b2(nH/L)+b3(n)+b4(nH/L)+c(-nc));
plot(Th, wait=1);
fespace Vh(Th,P2);
Vh u, uu, v, vv, uold=0, vold=0,l;
Vh a1, a2, w, s, f1=0, f2=0;
Vh g1=uold-a1, g2=vold-a2;
real[int] tmp(u[].n);
fespace Ph(Th,P1);
Ph p=0, pp;
problem MeshVelocity([a1,a2],[w,s]) //laplce eq: mesh MeshVelocity
=int2d(Th,qforder=3)(dx(a1)dx(w)+dy(a1)dy(w)+dx(a2)dx(s)+dy(a2)dy(s))
+on(1,2,a1=0,a2=0)
+on(3,a1=f1,a2=f2);
problem aleNS ([u,v,p],[uu,vv,pp],solver=UMFPACK)=
int2d(Th,qforder=3)(a(uuu+vvv)
+((uold-a1)dx(u)+(vold-a2)dy(u))uu+((uold-a1)dx(v)+(vold-a2)dy(v))vv
-p(dx(uu)+dy(vv)) //pressione
+invRe(dx(u)dx(uu)+dy(u)dy(uu)+dx(v)dx(vv)+dy(v)dy(vv))//termine diffusivo
+pp(dx(u)+dy(v))+erppp)
-int2d(Th,qforder=3)(a(uolduu+voldvv))
//-int2d(Th,qforder=3)((dx(a1)+dy(a2))(uuu+vvv))
+on(3, u=a1,v=a2)
+on(4,u=1,v=0)
+on(1,v=0);
//+on(b4,u=1,v=0)+on(b1,b3,v=0);
for(real t=0;t<T;t+=dt){
f2=-2pi0.250.161.1cos(2pi0.161.1*t);
f1=0.;
aleNS;
MeshVelocity;
Th=movemesh(Th,[x+dta1,y+dta2]);
l=dx(v)-dy(u)+dx(a2)-dy(a1);
tmp=u[]; uold=0; uold[]=tmp;
tmp=v[]; vold=0; vold[]=tmp;
tmp=a1[]; a1=0; a1[]=tmp;
tmp=a2[]; a2=0; a2[]=tmp;
tmp=g1[]; g1=0; g1[]=tmp;
tmp=g2[]; g2=0; g2[]=tmp;
//plot(coef=0.1,[u,v],value=1,fill=1);
plot(coef=0.1,l,value=1,fill=1);
//plot(coef=0.1,p,value=1,fill=1);
fespace Mh(Th, P0);
Mh h = hTriangle;
cout << "size of mesh = " << h[].max << endl;
cout<< "size of mesh = " << h[].min <<endl;
//Drug of cylinder
real cd=int1d(Th,3,qforder=3)(-4px+4invRe(2dx(g1)x+(dy(g1)+dx(g2))y));
//Lift of cylinder
real cl=int1d(Th,3,qforder=3)(-4py+4invRe*((dy(g1)+dx(g2))x+2dy(g2)*
Drug<<cd<<endl;
Lift<<cl<<endl;
real i=0;
i=i+1;
cout<<i<<endl;
cout<<t<<endl;
}