#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define M 63
#define N 7
#define M1 (M+1)
#define N1 (N+1)
#define Ly 1.0
#define Q 9
#define rho0 1.0
#define G 1.0e-4
#define Re 100.0
void mesh(void);
void lbini(void);
void data_read(void);
void Evol(void);
double feq(int k,double U, double V, double RHO);
void datadeal(void);
double dx,dt,rRe,w;
double f[M1][N1][Q],g[M1][N1][Q],u[M1][N1],v[M1][N1],rho[M1][N1],u0;
int e[Q][2],re[Q];
double tp[Q];
double U0;
int main()
{
int m,mmax,TEND=0,M2,N2;
double err;
M2=(M+1)/2;N2=(N+1)/2;
e[0][0]=e[0][1]=0;
e[1][0]=1; e[1][1]=0;
e[2][0]=0; e[2][1]=1;
e[3][0]=-1; e[3][1]=0;
e[4][0]=0; e[4][1]=-1;
e[5][0]=1; e[5][1]=1;
e[6][0]=-1; e[6][1]=1;
e[7][0]=-1; e[7][1]=-1;
e[8][0]=1; e[8][1]=-1;
tp[0]=4.0/9; tp[1]=tp[2]=tp[3]=tp[4]=1.0/9; tp[5]=tp[6]=tp[7]=tp[8]=1.0/36;
re[0]=0;re[1]=3;re[2]=4;re[3]=1;re[4]=2;re[5]=7;re[6]=8;re[7]=5;re[8]=6;
dt=dx=Ly/(M+1);
rRe=sqrt(Ly*G/(8*Re))*Ly;
U0=Re*rRe/Ly;
w=1.0/(3*rRe/dt+0.5);
lbini();
printf("Umax=%f w=%f\n",U0,w);
m=0;
err=1.0;
printf("input mmax:\n");
scanf("%d",&mmax);
TEND+=mmax ;
while((m<TEND)&&err>1.0e-6)
{
m++;
Evol();
if(m%100==0)
{
err=fabs(u[M2][N2]-u0)/fabs(u[M2][N2]);
u0=u[M2][N2];
printf("err=%e u_center=%e m=%d\n",err, u[M2][N2],m);
}
}
printf("Re=%f, Umax=%f\n",Re,U0);
datadeal();
return 0;
}
void lbini()
{
int i,j,k;
for(i=0;i<N1;i++)
for(j=0;j<M1;j++)
{
u[j][i]=v[j][i]=0.0;
rho[j][i]=rho0;
}
for(j=0;j<=M;j++) for(i=0;i<=N;i++) for(k=0;k<Q;k++)
f[j][i][k]=feq(k,u[j][i],v[j][i],rho[j][i]);
}
double feq(int k,double U, double V, double RHO)
{
double uv,eu,x;
eu=e[k][0]*U+e[k][1]*V;
uv=U*U+V*V;
x=tp[k]*RHO*(1.0+3*eu+4.5*eu*eu-1.5*uv);
return x;
}
void Evol()
{
int i,j,k,id,jd;
double FM;
// collison
for(j=0;j<=M;j++) for(i=0;i<=N;i++) for(k=0;k<Q;k++)
{
FM=feq(k,u[j][i],v[j][i],rho[j][i]);
g[j][i][k]=f[j][i][k]-w*(f[j][i][k]-FM)+dt*3.0*tp[k]*G*e[k][0];
}
//streaming
for(j=1;j<M;j++) for(i=0;i<=N;i++) for(k=0;k<Q;k++)
{
jd=j-e[k][1]; id=(N1+i-e[k][0])%N1;
f[j][i][k]=g[jd][id][k];
}
// Bounce-Back
j=0;
for(i=0;i<=N;i++)
{
for(k=0;k<Q;k++)
{
jd=j-e[k][1]; id=(N1+i-e[k][0])%N1;
if(jd>=0) f[j][i][k]=g[jd][id][k];
else f[j][i][k]=g[j][i][re[k]];
}
}
j=M;
for(i=0;i<=N;i++)
{
for(k=0;k<Q;k++)
{
jd=j-e[k][1]; id=(N1+i-e[k][0])%N1;
if(jd<=M) f[j][i][k]=g[jd][id][k];
else f[j][i][k]=g[j][i][re[k]];
}
}
// Flow variables
for(j=0;j<=M;j++) for(i=0;i<=N;i++)
{
rho[j][i]=f[j][i][0]+f[j][i][1]+f[j][i][2]+f[j][i][3]+f[j][i][4]+f[j][i][5]+f[j][i][6]+f[j][i][7]+f[j][i][8];
u[j][i]=(f[j][i][1]+f[j][i][5]+f[j][i][8]-f[j][i][3]-f[j][i][6]-f[j][i][7]);
v[j][i]=(f[j][i][5]+f[j][i][6]+f[j][i][2]-f[j][i][7]-f[j][i][8]-f[j][i][4]);
if(fabs(rho[j][i])>1.0e-9)
{
u[j][i]/=rho[j][i];
v[j][i]/=rho[j][i];
}
else u[j][i]=v[j][i]=0.0;
}
}
void datadeal()
{
int i,j;
double yj;
FILE *fp;
for(j=0;j<=M;j++)
{
yj=(j+0.5)*dx;
u[j][0]=4*U0*yj*(1-yj);
}
if( ( fp=fopen("uh.dat","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<=M;j++) {
for (i=0; i<=N; i++) fprintf(fp,"%e ",u[j][i]/U0);
fprintf(fp,"\n");
}
fclose(fp);
if( ( fp=fopen("vh.dat","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<=M;j++) {
for (i=0; i<=N; i++) fprintf(fp,"%e ",v[j][i]/U0);
fprintf(fp,"\n");
}
fclose(fp);
if( ( fp=fopen("rho.dat","w"))==NULL)
{
printf(" File Open Error\n");exit(1);
}
for(j=0;j<=M;j++) {
for (i=0; i<=N; i++) fprintf(fp,"%e ",rho[j][i]);
fprintf(fp,"\n");
}
fclose(fp);
}
评论0