📄 filmnew1.cpp
字号:
#include<math.h>
#include<stdlib.h>
#include<stdio.h>
int i,j,k,t,c,ny2;
int e[9][2]={{0,0},{1,0},{0,1},{-1,0},{0,-1},{1,1},{-1,1},{-1,-1},{1,-1}};
//const int ny1=109,nyl1=49;
const int ny=30,nyl=15,nx=30,time=1000;
const double D=4.0,A=0.0005;//w=1.5,
const double u0=0.0,v0=0.0,vss=1.0/3.0,tt=1.0;//the initial velocity and the double of the velosity of sound
const double tr[9]={{4.0/9.0},{1.0/9.0},{1.0/9.0},{1.0/9.0},{1.0/9.0},{1.0/36.0},{1.0/36.0},{1.0/36.0},{1.0/36.0}};
const double Tl=0.5,Tg=5.0;//the coefficient
const double rl=1.0,rg=0.001,p=0.4,kp=0.09992,Bt=0.02;//the initial density
double f[nx][ny][9],g[nx][ny][9],feq[nx][ny][9],geq[nx][ny][9];//the distributiong function
double u[nx][ny],v[nx][ny],uu,vv,fe,ge,pa;//the varieties of velocity and momentum
double den[nx][ny],nden[nx][ny],pr[nx][ny];// the variety of density
double temp,tempp,temp1[nx][ny],tempr2[nx][ny],temprh[nx][ny],temprx[nx][ny],tempry[nx][ny];
double tempdx[nx][ny][9],tempdy[nx][ny][9];
double tbc2[nx],tbc2g[nx],tbc5[nx],tbc5g[nx],tbc6[nx],tbc6g[nx];//the temporary varieties
double bbc4[nx],bbc4g[nx],bbc7[nx],bbc7g[nx],bbc8[nx],bbc8g[nx];//the temporary varieties
double lbc3[ny],lbc3g[ny],lbc6[ny],lbc6g[ny],lbc7[ny],lbc7g[ny];//the temporary varieties
double rbc1[ny],rbc1g[ny],rbc5[ny],rbc5g[ny],rbc8[ny],rbc8g[ny];//the temporary varieties
double summass;
double summon_u,summon_v,summon_p;
double TT,TTT,T,r,rm;//the gene of laxation
double cp[nx][ny],dxy[nx][ny],dx[nx][ny],dy[nx][ny];
double temp1x[nx][ny],temp1y[nx][ny];
double dx2[nx][ny],dy2[nx][ny],dxx[nx][ny],dyy[nx][ny],dxh[nx][ny],dyh[nx][ny];
double dd[nx][ny][9],dr2[nx][ny][9];
double f2[nx][ny][9];
int im[nx],imm[nx],jm[ny],jmm[ny],ip[nx],ipp[nx],jp[ny],jpp[ny];
void solv()
{
for(i=0;i<nx;i++)
{
if(i==0)
{
im[i]=nx-1;
imm[i]=nx-2;
ip[i]=i+1;
ipp[i]=i+2;
}
else if(i==1)
{
im[i]=i-1;
imm[i]=nx-1;
ip[i]=i+1;
ipp[i]=i+2;
}
else if(i==nx-1)
{
ip[i]=0;
ipp[i]=1;
im[i]=i-1;
imm[i]=i-2;
}
else if(i==nx-2)
{
ip[i]=i+1;
ipp[i]=0;
im[i]=i-1;
imm[i]=i-2;
}
else
{
ip[i]=i+1;
ipp[i]=i+2;
im[i]=i-1;
imm[i]=i-2;
}
}
for(j=0;j<ny;j++)
{
if(j==0)
{
jm[j]=1;
jmm[j]=2;
jp[j]=j+1;
jpp[j]=j+2;
}
else if(j==1)
{
jm[j]=j-1;
jmm[j]=1;
jp[j]=j+1;
jpp[j]=j+2;
}
else if(j==ny-1)
{
jp[j]=ny-2;
jpp[j]=ny-3;
jm[j]=j-1;
jmm[j]=j-2;
}
else if(j==ny-2)
{
jp[j]=j+1;
jpp[j]=ny-2;
jm[j]=j-1;
jmm[j]=j-2;
}
else
{
jp[j]=j+1;
jpp[j]=j+2;
jm[j]=j-1;
jmm[j]=j-2;
}
}
}
double cdxfu(double array[nx][ny],int m,int n)
{
return(1/3.0*(array[ip[m]][n]-array[im[m]][n])
+1/12.0*(array[ip[m]][jp[n]]-array[im[m]][jm[n]]
+array[ip[m]][jm[n]]-array[im[m]][jp[n]]));
}
double cdyfu(double array[nx][ny],int m,int n)
{
return(1/3.0*(array[m][jp[n]]-array[m][jm[n]])
+1/12.0*(array[ip[m]][jp[n]]-array[im[m]][jm[n]]
+array[im[m]][jp[n]]-array[ip[m]][jm[n]]));
}
double cdx3fu(double array1[nx][ny][9],int m,int n,int kk)
{
return(1/3.0*(array1[ip[m]][n][kk]-array1[im[m]][n][kk])
+1/12.0*(array1[ip[m]][jp[n]][kk]-array1[im[m]][jm[n]][kk]
+array1[ip[m]][jm[n]][kk]-array1[im[m]][jp[n]][kk]));
}
double cdy3fu(double array1[nx][ny][9],int m,int n,int kk)
{
return(1/3.0*(array1[m][jp[n]][kk]-array1[m][jm[n]][kk])
+1/12.0*(array1[ip[m]][jp[n]][kk]-array1[im[m]][jm[n]][kk]
+array1[im[m]][jp[n]][kk]-array1[ip[m]][jm[n]][kk]));
}
double ddfu(double array[nx][ny],int m,int n,int kk)
{
double array2;
int xpp[nx],xp[nx],ypp[ny],yp[ny],xm[nx],ym[ny];
if(e[kk][0]==0)
{
xpp[m]=m;
xp[m]=m;
xm[m]=m;
}
else if(e[kk][0]==1)
{
xpp[m]=ipp[m];
xp[m]=ip[m];
xm[m]=im[m];
}
else if(e[kk][0]==-1)
{
xpp[m]=imm[m];
xp[m]=im[m];
xm[m]=ip[m];
}
if(e[kk][1]==0)
{
ypp[n]=n;
yp[n]=n;
ym[n]=n;
}
else if(e[kk][1]==1)
{
ypp[n]=jpp[n];
yp[n]=jp[n];
ym[n]=jm[n];
}
else if(e[kk][1]==-1)
{
ypp[n]=jmm[n];
yp[n]=jm[n];
ym[n]=jp[n];
}
tempp=1/2.0*(-array[xpp[m]][ypp[n]]+4*array[xp[m]][yp[n]]-3*array[m][n]);
temp=1/2.0*(array[xpp[m]][ypp[n]]-array[m][n]);
temp=temp*tempp;
if(temp>=0.0)
array2=tempp;
else
array2=1/2.0*(array[xp[m]][yp[n]]-array[xm[m]][ym[n]]);
return(array2);
}
double ddcfu(double array[nx][ny],int m,int n,int kk)
{
double array2;
int xp[nx],yp[ny],xm[nx],ym[ny];
if(e[kk][0]==0)
{
xp[m]=m;
xm[m]=m;
}
else if(e[kk][0]==1)
{
xp[m]=ip[m];
xm[m]=im[m];
}
else if(e[kk][0]==-1)
{
xp[m]=im[m];
xm[m]=ip[m];
}
if(e[kk][1]==0)
{
yp[n]=n;
ym[n]=n;
}
else if(e[kk][1]==1)
{
yp[n]=jp[n];
ym[n]=jm[n];
}
else if(e[kk][1]==-1)
{
yp[n]=jm[n];
ym[n]=jp[n];
}
array2=1/2.0*(array[xp[m]][yp[n]]-array[xm[m]][ym[n]]);
return(array2);
}
void init()
{
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
r=0.5*(rl+rg)+0.5*(rl-rg)*tanh(-2.0*(j-nyl+1)/D);
uu=0.05;
vv=0;
pa=p;
temp=1.0-1.5*(uu*uu+vv*vv)+3.0*(e[k][0]*uu+e[k][1]*vv)+4.5*(e[k][0]*uu+e[k][1]*vv)*(e[k][0]*uu+e[k][1]*vv);
f[i][j][k]=tr[k]*r*temp;
g[i][j][k]=f[i][j][k]-tr[k]*r+tr[k]*pa/vss;
}
}
}
}
void eval()
{
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
den[i][j]=0;
nden[i][j]=0;
for(k=0;k<9;k++)
{
den[i][j]=den[i][j]+f[i][j][k];
nden[i][j]=nden[i][j]+g[i][j][k];
}
}
}
rm=(rl+rg)/2.0;
summass=0.0;
summon_u=0.0;
summon_v=0.0;
summon_p=0.0;
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
temp=den[i][j];
cp[i][j]=4.0*Bt*(temp-rg)*(temp-rl)*(temp-rm);
dx[i][j]=cdxfu(den,i,j);
dy[i][j]=cdyfu(den,i,j);
dxy[i][j]=1/6.0*(den[ip[i]][jp[j]]+den[im[i]][jp[j]]
+den[ip[i]][jm[j]]+den[im[i]][jm[j]]
+4.0*(den[ip[i]][j]+den[im[i]][j]+den[i][jp[j]]+
den[i][jm[j]])-20.0*den[i][j]);
temp1[i][j]=cp[i][j]-kp*dxy[i][j];
tempr2[i][j]=dx[i][j]*dx[i][j]+dy[i][j]*dy[i][j];
temprx[i][j]=dx[i][j]*dx[i][j];
temprh[i][j]=dx[i][j]*dy[i][j];
tempry[i][j]=dy[i][j]*dy[i][j];
}
}
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
temp1x[i][j]=cdxfu(temp1,i,j);
temp1y[i][j]=cdyfu(temp1,i,j);
dx2[i][j]=cdxfu(tempr2,i,j);
dy2[i][j]=cdyfu(tempr2,i,j);
dxx[i][j]=cdxfu(temprx,i,j);
dyh[i][j]=cdyfu(temprh,i,j);
dxh[i][j]=cdxfu(temprh,i,j);
dyy[i][j]=cdyfu(tempry,i,j);
u[i][j]=tt*kp/2.0*(dx2[i][j]-dxx[i][j]-dyh[i][j]);
v[i][j]=tt*kp/2.0*(dy2[i][j]-dxh[i][j]-dyy[i][j]);//-1/2.0*a*den[i][j];
for(k=0;k<9;k++)
{
u[i][j]=u[i][j]+e[k][0]*g[i][j][k];
v[i][j]=v[i][j]+e[k][1]*g[i][j][k];
}
uu=u[i][j]/den[i][j];
vv=v[i][j]/den[i][j];
temp=uu*dx[i][j]+vv*dy[i][j];
pr[i][j]=vss*(nden[i][j]+tt*vss/2.0*temp);
summass=summass+den[i][j];//the total quality
summon_u=summon_u+u[i][j];//the momentum in x direction
summon_v=summon_v+v[i][j];//the momentum in y direction
summon_p=summon_p+pr[i][j];
pa=pr[i][j];
for(k=0;k<9;k++)
{
temp=1.0-1.5*(uu*uu+vv*vv)+3.0*(e[k][0]*uu+e[k][1]*vv)+4.5*(e[k][0]*uu+e[k][1]*vv)*(e[k][0]*uu+e[k][1]*vv);
feq[i][j][k]=tr[k]*den[i][j]*temp;
geq[i][j][k]=feq[i][j][k]-tr[k]*den[i][j]+pa/vss*tr[k];
}
}
}
}
void prestreaming_collision()
{
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
dd[i][j][k]=ddfu(den,i,j,k);
tempdx[i][j][k]=dd[i][j][k]*dx[i][j];
tempdy[i][j][k]=dd[i][j][k]*dy[i][j];
}
}
}
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
dr2[i][j][k]=ddfu(tempr2,i,j,k);
f2[i][j][k]=ddfu(temp1,i,j,k);
tempdx[i][j][k]=cdx3fu(tempdx,i,j,k);
tempdy[i][j][k]=cdy3fu(tempdy,i,j,k);
uu=u[i][j]/den[i][j];
vv=v[i][j]/den[i][j];
fe=feq[i][j][k];
ge=geq[i][j][k];
if(den[i][j]>=0.731328)//j<=nyl-3&&j>=5)
{
T=Tl;
TT=1.0-1.0/(2.0*T);
}
else if(den[i][j]<=0.269672)//j>=nyl+1&&j<=ny-6)
{
T=Tg;
TT=1.0-1.0/(2.0*T);
}
else
{
T=(den[i][j]-rg)/(rl-rg)*Tl+(1.0-(den[i][j]-rg)/(rl-rg))*Tg;
TT=1.0-1.0/(2.0*T);
}
temp=1.0-1.5*(uu*uu+vv*vv)+3.0*(e[k][0]*uu+e[k][1]*vv)+4.5*(e[k][0]*uu+e[k][1]*vv)*(e[k][0]*uu+e[k][1]*vv);
f[i][j][k]=TT*f[i][j][k]+1.0/(2.0*T)*fe+tt*1.0/(2.0*vss)*tr[k]*temp*(dd[i][j][k]*vss
-den[i][j]*f2[i][j][k]-(uu*dx[i][j]+vv*dy[i][j])*vss
+den[i][j]*(uu*temp1x[i][j]+vv*temp1y[i][j]));//+1/(2.0*vss)*(e[k][1]-vv)*a*fe;
g[i][j][k]=TT*g[i][j][k]+1.0/(2.0*T)*ge+tt*1.0/2.0*(dd[i][j][k]-uu*dx[i][j]-vv*dy[i][j])
*(tr[k]*temp-tr[k])+tt*kp/(2.0*vss)*tr[k]*temp*(dr2[i][j][k]-uu*dx2[i][j]-vv*dy2[i][j]
-tempdx[i][j][k]-tempdy[i][j][k]+uu*(dxx[i][j]+dyh[i][j])
+vv*(dxh[i][j]+dyy[i][j]));//+1/(2.0*vss)*(e[k][1]-vv)*a*fe;
}
}
}
}
void poststreaming_collision()
{
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
dd[i][j][k]=ddcfu(den,i,j,k);
tempdx[i][j][k]=dd[i][j][k]*dx[i][j];
tempdy[i][j][k]=dd[i][j][k]*dy[i][j];
}
}
}
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
dr2[i][j][k]=ddcfu(tempr2,i,j,k);
f2[i][j][k]=ddcfu(temp1,i,j,k);
tempdx[i][j][k]=cdx3fu(tempdx,i,j,k);
tempdy[i][j][k]=cdy3fu(tempdy,i,j,k);
uu=u[i][j]/den[i][j];
vv=v[i][j]/den[i][j];
fe=feq[i][j][k];
ge=geq[i][j][k];
if(den[i][j]>=0.731328)//(j<=nyl-3&&j>=5)
{
T=Tl;
TT=1.0-1/(2.0*T+1.0);
TTT=2.0*T/(2.0*T+1.0);
}
else if(den[i][j]<=0.269672)//(j>=nyl+1&&j<=ny-6)
{
T=Tg;
TT=1.0-1/(2.0*T+1.0);
TTT=2.0*T/(2.0*T+1.0);
}
else
{
T=(den[i][j]-rg)/(rl-rg)*Tl+(1.0-(den[i][j]-rg)/(rl-rg))*Tg;
TT=1.0-1/(2.0*T+1.0);
TTT=2.0*T/(2.0*T+1.0);
}
temp=1.0-1.5*(uu*uu+vv*vv)+3.0*(e[k][0]*uu+e[k][1]*vv)+4.5*(e[k][0]*uu+e[k][1]*vv)*(e[k][0]*uu+e[k][1]*vv);
f[i][j][k]=TT*f[i][j][k]+1/(2.0*T+1)*fe+tt*1/(2.0*vss)*TTT*tr[k]*temp*(dd[i][j][k]*vss
-den[i][j]*f2[i][j][k]-(uu*dx[i][j]+vv*dy[i][j])*vss
+den[i][j]*(uu*temp1x[i][j]+vv*temp1y[i][j]))
;//+1/(2.0*vss)*(e[k][1]-vv)*a*fe*TTT;
g[i][j][k]=TT*g[i][j][k]+1/(2.0*T+1)*ge+tt*1/2.0*TTT*(dd[i][j][k]-uu*dx[i][j]-vv*dy[i][j])
*(tr[k]*temp-tr[k])+tt*kp/(2.0*vss)*TTT*tr[k]*temp*(dr2[i][j][k]-uu*dx2[i][j]-vv*dy2[i][j]
-tempdx[i][j][k]-tempdy[i][j][k]+uu*(dxx[i][j]+dyh[i][j])
+vv*(dxh[i][j]+dyy[i][j]));//+1/(2.0*vss)*(e[k][1]-vv)*a*fe*TTT;
}
}
}
}
void streaming()
{
for(i=0;i<nx;i++)
{
tbc2[i]=f[i][ny-1][2];
tbc6[i]=f[i][ny-1][6];
tbc5[i]=f[i][ny-1][5];
tbc2g[i]=g[i][ny-1][2];
tbc6g[i]=g[i][ny-1][6];
tbc5g[i]=g[i][ny-1][5];
bbc4[i]=f[i][0][4];
bbc7[i]=f[i][0][7];
bbc8[i]=f[i][0][8];
bbc4g[i]=g[i][0][4];
bbc7g[i]=g[i][0][7];
bbc8g[i]=g[i][0][8];
}
for(j=0;j<ny;j++)
{
lbc3[j]=f[0][j][3];
lbc6[j]=f[0][j][6];
lbc7[j]=f[0][j][7];
lbc3g[j]=g[0][j][3];
lbc6g[j]=g[0][j][6];
lbc7g[j]=g[0][j][7];
rbc1[j]=f[nx-1][j][1];
rbc5[j]=f[nx-1][j][5];
rbc8[j]=f[nx-1][j][8];
rbc1g[j]=g[nx-1][j][1];
rbc5g[j]=g[nx-1][j][5];
rbc8g[j]=g[nx-1][j][8];
}
for(i=nx-1;i>0;i--) //f1,g1
{
for(j=0;j<ny;j++)
{
f[i][j][1]=f[i-1][j][1];
g[i][j][1]=g[i-1][j][1];
}
}
for(j=ny-1;j>0;j--) //f2,g2
{
for(i=0;i<nx;i++)
{
f[i][j][2]=f[i][j-1][2];
g[i][j][2]=g[i][j-1][2];
}
}
for(i=0;i<nx-1;i++) //f3,g3
{
for(j=0;j<ny;j++)
{
f[i][j][3]=f[i+1][j][3];
g[i][j][3]=g[i+1][j][3];
}
}
for(j=0;j<ny-1;j++) //f4,g4
{
for(i=0;i<nx;i++)
{
f[i][j][4]=f[i][j+1][4];
g[i][j][4]=g[i][j+1][4];
}
}
for(i=nx-1;i>0;i--) //f5,g5
{
for(j=ny-1;j>0;j--)
{
f[i][j][5]=f[i-1][j-1][5];
g[i][j][5]=g[i-1][j-1][5];
}
}
for(i=0;i<nx-1;i++) //f6,g6
{
for(j=ny-1;j>0;j--)
{
f[i][j][6]=f[i+1][j-1][6];
g[i][j][6]=g[i+1][j-1][6];
}
}
for(i=0;i<nx-1;i++) //f7,g7
{
for(j=0;j<ny-1;j++)
{
f[i][j][7]=f[i+1][j+1][7];
g[i][j][7]=g[i+1][j+1][7];
}
}
for(i=nx-1;i>0;i--) //f8,g8
{
for(j=0;j<ny-1;j++)
{
f[i][j][8]=f[i-1][j+1][8];
g[i][j][8]=g[i-1][j+1][8];
}
}
}
void BL_condition()
{
for(i=0;i<nx;i++)
{
f[i][0][2]=bbc4[i];
g[i][0][2]=bbc4g[i];
f[i][ny-1][4]=tbc2[i];
g[i][ny-1][4]=tbc2g[i];
}
for(i=0;i<nx;i++)
{
f[i][ny-1][7]=tbc6[i];
g[i][ny-1][7]=tbc6g[i];
f[i][0][6]=bbc7[i];
g[i][0][6]=bbc7g[i];
}
for(i=0;i<nx;i++)
{
f[i][ny-1][8]=tbc5[i];
g[i][ny-1][8]=tbc5g[i];
f[i][0][5]=bbc8[i];
g[i][0][5]=bbc8g[i];
}
for(j=0;j<ny;j++)
{
f[0][j][1]=rbc1[j];
f[nx-1][j][3]=lbc3[j];
g[0][j][1]=rbc1g[j];
g[nx-1][j][3]=lbc3g[j];
}
for(j=ny-1;j>0;j--)
{
f[0][j][5]=rbc5[j-1];
g[0][j][5]=rbc5g[j-1];
f[nx-1][j][6]=lbc6[j-1];
g[nx-1][j][6]=lbc6g[j-1];
}
for(j=0;j<ny-1;j++)
{
f[0][j][8]=rbc8[j+1];
g[0][j][8]=rbc8g[j+1];
f[nx-1][j][7]=lbc7[j+1];
g[nx-1][j][7]=lbc7g[j+1];
}
}
void main()
{
printf("input c \n");
scanf("%d",&c);
if(c==1)
{
FILE*aq;
aq=fopen("out1.dat","r");
for(i=3;i<nx-3;i++)
{
for(j=3;j<ny-3;j++)
{
for(k=0;k<9;k++)
{
fscanf(aq,"%lf\n",&f[i][j][k]);
}
}
}
fclose(aq);
FILE*ay;
ay=fopen("out2.dat","r");
for(i=3;i<nx-3;i++)
{
for(j=3;j<ny-3;j++)
{
for(k=0;k<9;k++)
{
fscanf(ay,"%lf\n",&g[i][j][k]);
}
}
}
fclose(ay);
}
else
init();
solv();
for(t=1;t<=time;t++)
{
eval();
poststreaming_collision();
streaming();
BL_condition();
eval();
// if(t%100==0)
printf(" %d %lf,%lf,%lf,%lf\n",t,summass,summon_u,summon_v,summon_p);
prestreaming_collision();
streaming();
BL_condition();
}
eval();
poststreaming_collision();
streaming();
BL_condition();
eval();
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
u[i][j]=u[i][j]/den[i][j];
v[i][j]=v[i][j]/den[i][j];
}
}
FILE*ax;
ax=fopen("out.dat","w");
fprintf(ax,"variable=x y u v p r \n");
fprintf(ax,"zone i=%d j=%d\n",nx,ny);
for(j=0;j<ny;j++)
{
for(i=0;i<nx;i++)
{
fprintf(ax,"%d %d %f %f %f %f\n",i,j,u[i][j],v[i][j],pr[i][j],den[i][j]);
}
}
FILE*as;
as=fopen("out1.dat","w");
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
fprintf(as,"%f\n",f[i][j][k]);
}
}
}
FILE*ar;
ar=fopen("out2.dat","w");
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
for(k=0;k<9;k++)
{
fprintf(ar,"%f\n",g[i][j][k]);
}
}
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -