Hello,
I am sorry to post on this topic but I am trying to solve the compressible Navier-Stokes equations with Freefem and the continuity and momentum equations are quite the same to the one proposed by M. Guangyuan Wei (except the triple derivative of rho). I have used the following equations :
I have added the equation for energy and the dependance of pressure with respect to temperature and density (ideal gas law), the viscosity as a function of temperature and diffusivity based on a constant Prandtl number. I have tried to write the corresponding weak formulation to solve these equations. In order to try the implementation, I have made a test on the flow around a cylinder at low Reynolds number (212.6).
I would like to know if the finite element space I used are the correct one to solve such kind of equations or if other implementation should be used.
I thank you for your help.
real xmin = -100;
real xmax = 200;
real rmax = 25;
// 1 : axis , 2 : outlet, 3 : cylinder, 4 : inlet, 5 : refinement zone, 6 : top
border a00(t=xmin,-2.5) {x=t; y=0; label=1;}
border a01(t=-2.5,-0.5) {x=t; y=0; label=1;}
border a02(t=pi,0) {x=0.5*cos(t); y=0.5*sin(t); label=3;}
border a03(t=0.5,5.25) {x=t; y=0; label=1;}
border a04(t=0,2) {x=5.25; y=t; label=5;}
border a05(t=5.25,-2.5) {x=t; y=2; label=5;}
border a06(t=2,0) {x=-2.5; y=t; label=5;}
border a07(t=5.25,xmax) {x=t; y=0; label=1;}
border a08(t=0,rmax) {x=xmax; y=t; label=2;}
border a09(t=xmax,xmin) {x=t; y=rmax; label=6;}
border a10(t=rmax,0) {x=xmin; y=t; label=4;}
// Define an initial mesh
int n=3;
mesh Th = buildmesh(a00(20*n)+a01(20*n)+a02(50*n)+a03(20*n)+a04(-20*n)+a05(-40*n)+a06(-20*n)+a07(90*n)+a08(10*n)+a09(60*n)+a10(10*n));
// Save mesh features
{ ofstream file("SOLUT4/coordinates.dat");
for (int j=0;j<Th.nv; j++) {
file << Th(j).x << " " << Th(j).y<< endl;}
}
{ ofstream file("SOLUT4/connectivity.dat");
int nbtriangle = Th.nt;
for (int i=0;i<Th.nt; i++){
file << Th[i][0]+1 << " " << Th[i][1]+1 << " " << Th[i][2]+1 << endl;}
}
real T=1.,dt=0.000005; // total time and time step
real Re=212.6; // Reynolds number (characeristic size =1 and characteristic velocity = 1)
real visco=1.0/Re;
real gamma=1.4;
real rgaz=287.0;
real cp=(gamma*rgaz)/(gamma-1.0); // perfect gas constant
real prandtl=0.7;
//pressure based on the ideal gas law
func real pression(real argrho, real argtemp){
return argrho*rgaz*argtemp;
}
//viscosity as a function of temperature, Sutherland law
func real viscofunc(real argtemp){
return 1.715e-5*(argtemp/273.15)^(3.0/2.0)*((273.15+110.4)/(argtemp+110.4));
}
//conductivity as a function of viscosity considering a constant Prandtl number
func real kappafunc(real argvisco){
return argvisco*cp/prandtl;
}
fespace Rh(Th,P1); Rh rho=1.0,rhot,rhoold,press,temp=300.0,tempt,tempold;
fespace Uh(Th,P1); Uh u=0.0,v=0.0,uu,vv,uold,vold;
problem NSK([rho,u,v,temp],[rhot,uu,vv,tempt])=
int2d(Th)(rho*rhot/dt)-int2d(Th)((convect([uold,vold],-dt,rhoold)*rhot/dt)) +int2d(Th)(rhoold*(dx(u)+dy(v))*rhot) //continuity equation
+int2d(Th)(rhoold*(u*uu+v*vv)/dt)-int2d(Th)(rhoold*(convect([uold,vold],-dt,uold)*uu/dt))-int2d(Th)(convect([uold,vold],-dt,vold)*vv/dt)
-int2d(Th)(pression(rhoold,tempold)*(dx(uu)+dy(vv)))
+int2d(Th)(viscofunc(tempold)*(((4./3.)*dx(u)-(2./3.)*dy(v))*dx(uu)+(dx(v)+dy(u))*(dx(vv)+dy(uu))+((4./3.)*dy(v)-(2./3.)*dx(u))*dy(vv))) // momentum equations
+int2d(Th)(rhoold*cp*(temp*tempt/dt))-int2d(Th)(rhoold*cp*(convect([uold,vold],-dt,tempold)*tempt/dt))
+int2d(Th)(kappafunc(viscofunc(tempold))*(dx(temp)*dx(tempt)+dy(temp)*dy(tempt)))
+int2d(Th)(pression(rhoold,tempold)*dx(u)*tempt+pression(rhoold,tempold)*dy(v)*tempt)
-int2d(Th)(viscofunc(tempold)*(uold*((4./3.)*dx(u)-(2./3.)*dy(v)))*dx(tempt)+viscofunc(tempold)*(vold*(dx(v)+dy(u)))*dx(tempt)+viscofunc(tempold)*(uold*(dx(v)+dy(u)))*dy(tempt)+viscofunc(tempold)*(vold*((4./3.)*dy(v)-(2./3.)*dx(u)))*dy(tempt)) // Temperature equation
+ on(2,4,6,u=1.0,v=0.0,rho=1.0,temp=300.0) + on(3,u=0.0,v=0.0) + on(1,v=0.0); // boundary conditions
int i=0;
for(real t=0.;t<T;t+=dt)
{
uold=u;
vold=v;
rhoold=rho;
tempold=temp;
press=rho*rgaz*temp;
NSK;
if (i%10==0){
{
{ofstream fid("SOLUT4/solut_ux_"+i+".dat");
for (int j=0; j<u[].n; j++)
fid << u[][j] << endl;
}
}
{ofstream fid("SOLUT4/solut_uy_"+i+".dat");
for (int j=0; j<v[].n; j++)
fid << v[][j] << endl;
}
}
{
{ofstream fid("SOLUT4/solut_rho_"+i+".dat");
for (int j=0; j<rho[].n; j++)
fid << rho[][j] << endl;
}
}
{
{ofstream fid("SOLUT4/solut_p_"+i+".dat");
for (int j=0; j<press[].n; j++)
fid << press[][j] << endl;
}
}
cout << "time " << t+dt << " min rho = " << rho[].min << " max rho = " << rho[].max << endl;
i=i+1;
}
