/* 下放为多孔介质区域，上方为电解液区域

*/
#include "stdio.h"
#include "math.h"
#include "time.h"
#define n 3200
#define m 2400
int kk=0;
void main(){
static double f[9][n+1][m+1];
static double feq[9][n+1][m+1];
static double rho[n+1][m+1];
static double u[n+1][m+1];
static double v[n+1][m+1];//声明变量和函数
double w[9],cx[9],cy[9],uo=0.0000279,sumvelo=0.0,rhoo=5.0,dx=1.0,dy=1.0,dt=1.0,cs,alpha=0.01,Re,omega;
double kv,K,dp,Fx[n+1][m+1],Fy[n+1][m+1],F[9][n+1][m+1];
static int i,j,k,mstep=40000;
static int tag[n+1][m+1];
FILE *fp_0;
char judge;
FILE *fp_tag;
FILE *fp_result_u,*fp_result_v,*fp_rho;
int collesion();
int streaming();
int sfbound();
int rhouv();
void system();//声明库函数，此项实际上可有可无，不声明会有警告
void exit();
clock_t start, finish;
//判定是否进行过计算
if((fp_0=fopen("result_u.txt","r"))!=NULL)
{
for(;;)
{
printf("已计算过流场，是否重新计算？（y/n）\n");
scanf("%s",&judge);
if(judge=='n')
{
fclose(fp_0);
exit(0);
}
else if(judge=='y') break;
}
}

//printf("输入入口速度:\n");scanf("%lf",&uo);
//printf("输入迭代次数:\n");scanf("%d",&mstep);

//打开标记文件
if((fp_tag=fopen("tag.txt","r"))==NULL)
{
printf("无法打开文件\"tag.txt\",请确保文件存在!\n");
system("pause");
exit(0);
}

//创建记录文件
fp_result_u=fopen("result_u.txt","w+");
fp_result_v=fopen("result_v.txt","w+");
fp_rho=fopen("rho.txt","w+");

//读取固相标记
for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fscanf(fp_tag,"%d",&tag[i][j]);
}

//计算基本参数
dp=0.0006;
K=0.0000012;
Re=uo*n/alpha;
omega=1.0/(3.0*alpha+0.5);
kv=((1.0/(3.0*alpha+0.5))-0.5)/3;
cs=(((dx*dx)/(dt*dt))+((dy*dy)/(dt*dt)))/3;
printf("基本参数:\nRe=%f omega=%f alpha=%f cs=%f\n\n",Re,omega,alpha,cs);

//给w赋值
w[0]=4.0/9;for(i=1;i<=8;i++)
{
if(i<=4) w[i]=1.0/9;else w[i]=1.0/36;
}
printf("w0~w8:\nw[6]=%5f w[2]=%5f w[5]=%5f\nw[3]=%5f w[0]=%5f w[1]=%5f\nw[7]=%5f w[4]=%5f w[8]=%5f\n\n",w[6],w[2],w[5],w[3],w[0],w[1],w[7],w[4],w[8]);

//给cx、cy赋值
cx[0]=0,cx[1]=1,cx[2]=0,cx[3]=-1,cx[4]=0,cx[5]=1,cx[6]=-1,cx[7]=-1,cx[8]=1;
cy[0]=0,cy[1]=0,cy[2]=1,cy[3]=0,cy[4]=-1,cy[5]=1,cy[6]=1,cy[7]=-1,cy[8]=-1;

//初始值
for(j=0;j<=m;j++)
{
for(i=0;i<=n;i++)
{
rho[i][j]=rhoo;
u[i][j]=0.0;
v[i][j]=0.0;
}
}
for(j=601;j<=m-1;j++) u[0][j]=uo;//左边界进口速度

``````printf("开始迭代。。。\n");

start=clock();
for(kk=1;kk<=mstep;kk++)//循环迭代
{
collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m,cs,Fx,Fy,F,kv,K);
streaming(f,n,m);
sfbound(f,n,m,uo,tag);
rhouv(f,rho,u,v,cx,cy,n,m,tag);
if(kk%100==0) printf("总迭代次数为%5d 次，已迭代%5d 次。。。\n",mstep,kk);
if(kk%100==0) printf("观察值:u[150][70]=%7.11f v[150][70]=%7.11f\n\n",u[150][70],v[150][70]);
if(kk%100==0) printf("观察值:u[150][220]=%7.11f v[150][220]=%7.11f\n\n",u[150][220],v[150][220]);
if(kk==mstep) printf("迭代完毕！ \n");
}
finish=clock();
printf("\n迭代耗时%.3fs\n",(double)(finish - start) / CLOCKS_PER_SEC);

for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fprintf(fp_result_u,"%7.11f ",u[i][j]);
fprintf(fp_result_u,"\n");
}
for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fprintf(fp_result_v,"%7.11f ",v[i][j]);
fprintf(fp_result_v,"\n");
}
for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fprintf(fp_rho,"%9.6f ",rho[i][j]);
fprintf(fp_rho,"\n");
}

printf("\n观察值:u[2][m-2]=%7.11f v[2][m-2]=%7.11f\n",u[2][m-2],v[2][m-2]);//观察值

fclose(fp_result_u);
fclose(fp_result_v);
fclose(fp_rho);
fclose(fp_tag);//关闭文件
printf("\n程序运行完毕,");
system("pause");
``````

}//主函数结束

//碰撞过程
int collesion(double u[n+1][m+1],double v[n+1][m+1],double f[9][n+1][m+1],double feq[9][n+1][m+1],double rho[n+1][m+1],
double Fx[n+1][m+1],double Fy[n+1][m+1],double F[9][n+1][m+1],double kv,double K,double omega,double w[9],
double cx[9],double cy[9],int a,int b,double cs){
int i,j,k;
double t1,t2,t3,c1,c2;
t3=cs*cs;
for(i=0;i<=a;i++)
{
for(j=0;j<=b;j++)
{
c1=Fx[i][j]*cx[k]+Fy[i][j]*cy[k];
t1=u[i][j]*u[i][j]+v[i][j]*v[i][j];
for(k=0;k<=8;k++)
{
c2=u[i][j]*Fx[i][j]+v[i][j]*Fy[i][j];
t2=u[i][j]*cx[k]+v[i][j]*cy[k];
Fx[i][j]=-(kv/0.0000012)*u[i][j];
Fy[i][j]=-(kv/0.0000012)*v[i][j];
F[k][i][j]=rho[i][j]*(1-omega/2.0)*(3.0*c1+20.0*c1*t2-20.0*c2/3);
feq[k][i][j]=rho[i][j]*w[k]*(1.0+3.0*t2+10.0*t2*t2-10.0*t1/3);//当ck=1时，cs*cs=1/3
f[k][i][j]=omega*feq[k][i][j]+(1.0-omega)*f[k][i][j]+F[k][i][j];
//if(i==2&&j==2&&k==1) printf("第%d次迭代结束:\n坐标(%d,%d)处t1=%9.6f t2=%9.6f t3=%9.6f u=%9.6f v=%9.6f\n",i,j,kk,t1,t2,t3,u[2][m-2],v[2][m-2]);
}
}
}
return(0);}

//迁移过程
int streaming(double f[9][n+1][m+1],int a,int b){
int i,j;
for(j=0;j<=b;j++)
{
for(i=a;i>=1;i--) f[1][i][j]=f[1][i-1][j];
for(i=0;i<=a-1;i++) f[3][i][j]=f[3][i+1][j];
}
for(j=b;j>=1;j--)
{
for(i=0;i<=a;i++) f[2][i][j]=f[2][i][j-1];
for(i=a;i>=1;i--) f[5][i][j]=f[5][i-1][j-1];
for(i=0;i<=a-1;i++) f[6][i][j]=f[6][i+1][j-1];
}
for(j=0;j<=b-1;j++)
{
for(i=0;i<=a;i++) f[4][i][j]=f[4][i][j+1];
for(i=0;i<=a-1;i++) f[7][i][j]=f[7][i+1][j+1];
for(i=a;i>=1;i--) f[8][i][j]=f[8][i-1][j+1];
}
return(0);
}

//边界
int sfbound(double f[9][n+1][m+1],int a,int b,double uo,int tag[n+1][m+1]){
int i=0,j=0;double rhow;
for(j=1;j<=b-1;j++)
{
for(i=1;i<=a-1;i++)
{
if(tag[i][j]==0)
{
if(tag[i-1][j]==1) f[1][i-1][j]=f[3][i-1][j];
if(tag[i][j-1]==1) f[2][i][j-1]=f[4][i][j-1];
if(tag[i+1][j]==1) f[3][i+1][j]=f[1][i+1][j];
if(tag[i][j+1]==1) f[4][i][j+1]=f[2][i][j+1];
if(tag[i-1][j-1]==1) f[5][i-1][j-1]=f[7][i-1][j-1];
if(tag[i+1][j-1]==1) f[6][i+1][j-1]=f[8][i+1][j-1];
if(tag[i+1][j+1]==1) f[7][i+1][j+1]=f[5][i+1][j+1];
if(tag[i-1][j+1]==1) f[8][i-1][j+1]=f[6][i-1][j+1];
}
}
}
for(j=201;j<=b-1;j++)
{
rhow=(f[0][0][j]+f[2][0][j]+f[4][0][j]+2*(f[3][0][j]+f[6][0][j]+f[7][0][j]))/(1.0-uo);
f[1][0][j]=f[3][0][j]+2.0*rhow*uo/3.0;
f[5][0][j]=f[7][0][j]+rhow*uo/6.0;
f[8][0][j]=f[6][0][j]+rhow*uo/6.0;//左边界
}

``````//for(j=0;j<=b;j++) {   f[3][a][j]=2*f[3][a-1][j]-f[3][a-2][j]; f[6][a][j]=2*f[6][a-1][j]-f[6][a-2][j]; f[7][a][j]=2*f[7][a-1][j]-f[7][a-2][j];}

return(0);
``````

}

//重新求解U、V
int rhouv(double f[9][n+1][m+1],double rho[n+1][m+1],double u[n+1][m+1],double v[n+1][m+1],double cx[9],double cy[9],int a,int b,int tag[n+1][m+1]){
int i,j,k;
double ssum;
double usum;
double vsum;
for(j=0;j<=b;j++)
{
for(i=0;i<=a;i++)
{
ssum=0.0;
for(k=0;k<=8;k++)
{
ssum=ssum+f[k][i][j];
}
rho[i][j]=ssum;
}
}
for(i=1;i<=a;i++)
{
for(j=1;j<=b-1;j++)
{
usum=0.0;
vsum=0.0;
for(k=0;k<=8;k++)
{
usum=usum+f[k][i][j]*cx[k];
vsum=vsum+f[k][i][j]*cy[k];
}
u[i][j]=usum/rho[i][j];
v[i][j]=vsum/rho[i][j];
if(tag[i][j]==1)
{
u[i][j]=0.0;
v[i][j]=0.0;
}
}
}
for(j=0;j<=b;j++) v[a][j]=0.0;
return(0);
}
/* 下放为多孔介质区域，上方为电解液区域

*/
#include "stdio.h"
#include "math.h"
#include "time.h"
#define n 3200
#define m 2400
int kk=0;
void main(){
static double f[9][n+1][m+1];
static double feq[9][n+1][m+1];
static double rho[n+1][m+1];
static double u[n+1][m+1];
static double v[n+1][m+1];//声明变量和函数
double w[9],cx[9],cy[9],uo=0.0000279,sumvelo=0.0,rhoo=5.0,dx=1.0,dy=1.0,dt=1.0,cs,alpha=0.01,Re,omega;
double kv,K,dp,Fx[n+1][m+1],Fy[n+1][m+1],F[9][n+1][m+1];
static int i,j,k,mstep=40000;
static int tag[n+1][m+1];
FILE *fp_0;
char judge;
FILE *fp_tag;
FILE *fp_result_u,*fp_result_v,*fp_rho;
int collesion();
int streaming();
int sfbound();
int rhouv();
void system();//声明库函数，此项实际上可有可无，不声明会有警告
void exit();
clock_t start, finish;
//判定是否进行过计算
if((fp_0=fopen("result_u.txt","r"))!=NULL)
{
for(;;)
{
printf("已计算过流场，是否重新计算？（y/n）\n");
scanf("%s",&judge);
if(judge=='n')
{
fclose(fp_0);
exit(0);
}
else if(judge=='y') break;
}
}

//printf("输入入口速度:\n");scanf("%lf",&uo);
//printf("输入迭代次数:\n");scanf("%d",&mstep);

//打开标记文件
if((fp_tag=fopen("tag.txt","r"))==NULL)
{
printf("无法打开文件\"tag.txt\",请确保文件存在!\n");
system("pause");
exit(0);
}

//创建记录文件
fp_result_u=fopen("result_u.txt","w+");
fp_result_v=fopen("result_v.txt","w+");
fp_rho=fopen("rho.txt","w+");

//读取固相标记
for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fscanf(fp_tag,"%d",&tag[i][j]);
}

//计算基本参数
dp=0.0006;
K=0.0000012;
Re=uo*n/alpha;
omega=1.0/(3.0*alpha+0.5);
kv=((1.0/(3.0*alpha+0.5))-0.5)/3;
cs=(((dx*dx)/(dt*dt))+((dy*dy)/(dt*dt)))/3;
printf("基本参数:\nRe=%f omega=%f alpha=%f cs=%f\n\n",Re,omega,alpha,cs);

//给w赋值
w[0]=4.0/9;for(i=1;i<=8;i++)
{
if(i<=4) w[i]=1.0/9;else w[i]=1.0/36;
}
printf("w0~w8:\nw[6]=%5f w[2]=%5f w[5]=%5f\nw[3]=%5f w[0]=%5f w[1]=%5f\nw[7]=%5f w[4]=%5f w[8]=%5f\n\n",w[6],w[2],w[5],w[3],w[0],w[1],w[7],w[4],w[8]);

//给cx、cy赋值
cx[0]=0,cx[1]=1,cx[2]=0,cx[3]=-1,cx[4]=0,cx[5]=1,cx[6]=-1,cx[7]=-1,cx[8]=1;
cy[0]=0,cy[1]=0,cy[2]=1,cy[3]=0,cy[4]=-1,cy[5]=1,cy[6]=1,cy[7]=-1,cy[8]=-1;

//初始值
for(j=0;j<=m;j++)
{
for(i=0;i<=n;i++)
{
rho[i][j]=rhoo;
u[i][j]=0.0;
v[i][j]=0.0;
}
}
for(j=601;j<=m-1;j++) u[0][j]=uo;//左边界进口速度

``````printf("开始迭代。。。\n");

start=clock();
for(kk=1;kk<=mstep;kk++)//循环迭代
{
collesion(u,v,f,feq,rho,omega,w,cx,cy,n,m,cs,Fx,Fy,F,kv,K);
streaming(f,n,m);
sfbound(f,n,m,uo,tag);
rhouv(f,rho,u,v,cx,cy,n,m,tag);
if(kk%100==0) printf("总迭代次数为%5d 次，已迭代%5d 次。。。\n",mstep,kk);
if(kk%100==0) printf("观察值:u[150][70]=%7.11f v[150][70]=%7.11f\n\n",u[150][70],v[150][70]);
if(kk%100==0) printf("观察值:u[150][220]=%7.11f v[150][220]=%7.11f\n\n",u[150][220],v[150][220]);
if(kk==mstep) printf("迭代完毕！ \n");
}
finish=clock();
printf("\n迭代耗时%.3fs\n",(double)(finish - start) / CLOCKS_PER_SEC);

for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fprintf(fp_result_u,"%7.11f ",u[i][j]);
fprintf(fp_result_u,"\n");
}
for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fprintf(fp_result_v,"%7.11f ",v[i][j]);
fprintf(fp_result_v,"\n");
}
for(j=m;j>=0;j--)
{
for(i=0;i<=n;i++) fprintf(fp_rho,"%9.6f ",rho[i][j]);
fprintf(fp_rho,"\n");
}

printf("\n观察值:u[2][m-2]=%7.11f v[2][m-2]=%7.11f\n",u[2][m-2],v[2][m-2]);//观察值

fclose(fp_result_u);
fclose(fp_result_v);
fclose(fp_rho);
fclose(fp_tag);//关闭文件
printf("\n程序运行完毕,");
system("pause");
``````

}//主函数结束

//碰撞过程
int collesion(double u[n+1][m+1],double v[n+1][m+1],double f[9][n+1][m+1],double feq[9][n+1][m+1],double rho[n+1][m+1],
double Fx[n+1][m+1],double Fy[n+1][m+1],double F[9][n+1][m+1],double kv,double K,double omega,double w[9],
double cx[9],double cy[9],int a,int b,double cs){
int i,j,k;
double t1,t2,t3,c1,c2;
t3=cs*cs;
for(i=0;i<=a;i++)
{
for(j=0;j<=b;j++)
{
c1=Fx[i][j]*cx[k]+Fy[i][j]*cy[k];
t1=u[i][j]*u[i][j]+v[i][j]*v[i][j];
for(k=0;k<=8;k++)
{
c2=u[i][j]*Fx[i][j]+v[i][j]*Fy[i][j];
t2=u[i][j]*cx[k]+v[i][j]*cy[k];
Fx[i][j]=-(kv/0.0000012)*u[i][j];
Fy[i][j]=-(kv/0.0000012)*v[i][j];
F[k][i][j]=rho[i][j]*(1-omega/2.0)*(3.0*c1+20.0*c1*t2-20.0*c2/3);
feq[k][i][j]=rho[i][j]*w[k]*(1.0+3.0*t2+10.0*t2*t2-10.0*t1/3);//当ck=1时，cs*cs=1/3
f[k][i][j]=omega*feq[k][i][j]+(1.0-omega)*f[k][i][j]+F[k][i][j];
//if(i==2&&j==2&&k==1) printf("第%d次迭代结束:\n坐标(%d,%d)处t1=%9.6f t2=%9.6f t3=%9.6f u=%9.6f v=%9.6f\n",i,j,kk,t1,t2,t3,u[2][m-2],v[2][m-2]);
}
}
}
return(0);}

//迁移过程
int streaming(double f[9][n+1][m+1],int a,int b){
int i,j;
for(j=0;j<=b;j++)
{
for(i=a;i>=1;i--) f[1][i][j]=f[1][i-1][j];
for(i=0;i<=a-1;i++) f[3][i][j]=f[3][i+1][j];
}
for(j=b;j>=1;j--)
{
for(i=0;i<=a;i++) f[2][i][j]=f[2][i][j-1];
for(i=a;i>=1;i--) f[5][i][j]=f[5][i-1][j-1];
for(i=0;i<=a-1;i++) f[6][i][j]=f[6][i+1][j-1];
}
for(j=0;j<=b-1;j++)
{
for(i=0;i<=a;i++) f[4][i][j]=f[4][i][j+1];
for(i=0;i<=a-1;i++) f[7][i][j]=f[7][i+1][j+1];
for(i=a;i>=1;i--) f[8][i][j]=f[8][i-1][j+1];
}
return(0);
}

//边界
int sfbound(double f[9][n+1][m+1],int a,int b,double uo,int tag[n+1][m+1]){
int i=0,j=0;double rhow;
for(j=1;j<=b-1;j++)
{
for(i=1;i<=a-1;i++)
{
if(tag[i][j]==0)
{
if(tag[i-1][j]==1) f[1][i-1][j]=f[3][i-1][j];
if(tag[i][j-1]==1) f[2][i][j-1]=f[4][i][j-1];
if(tag[i+1][j]==1) f[3][i+1][j]=f[1][i+1][j];
if(tag[i][j+1]==1) f[4][i][j+1]=f[2][i][j+1];
if(tag[i-1][j-1]==1) f[5][i-1][j-1]=f[7][i-1][j-1];
if(tag[i+1][j-1]==1) f[6][i+1][j-1]=f[8][i+1][j-1];
if(tag[i+1][j+1]==1) f[7][i+1][j+1]=f[5][i+1][j+1];
if(tag[i-1][j+1]==1) f[8][i-1][j+1]=f[6][i-1][j+1];
}
}
}
for(j=201;j<=b-1;j++)
{
rhow=(f[0][0][j]+f[2][0][j]+f[4][0][j]+2*(f[3][0][j]+f[6][0][j]+f[7][0][j]))/(1.0-uo);
f[1][0][j]=f[3][0][j]+2.0*rhow*uo/3.0;
f[5][0][j]=f[7][0][j]+rhow*uo/6.0;
f[8][0][j]=f[6][0][j]+rhow*uo/6.0;//左边界
}

``````//for(j=0;j<=b;j++) {   f[3][a][j]=2*f[3][a-1][j]-f[3][a-2][j]; f[6][a][j]=2*f[6][a-1][j]-f[6][a-2][j]; f[7][a][j]=2*f[7][a-1][j]-f[7][a-2][j];}

return(0);
``````

}

//重新求解U、V
int rhouv(double f[9][n+1][m+1],double rho[n+1][m+1],double u[n+1][m+1],double v[n+1][m+1],double cx[9],double cy[9],int a,int b,int tag[n+1][m+1]){
int i,j,k;
double ssum;
double usum;
double vsum;
for(j=0;j<=b;j++)
{
for(i=0;i<=a;i++)
{
ssum=0.0;
for(k=0;k<=8;k++)
{
ssum=ssum+f[k][i][j];
}
rho[i][j]=ssum;
}
}
for(i=1;i<=a;i++)
{
for(j=1;j<=b-1;j++)
{
usum=0.0;
vsum=0.0;
for(k=0;k<=8;k++)
{
usum=usum+f[k][i][j]*cx[k];
vsum=vsum+f[k][i][j]*cy[k];
}
u[i][j]=usum/rho[i][j];
v[i][j]=vsum/rho[i][j];
if(tag[i][j]==1)
{
u[i][j]=0.0;
v[i][j]=0.0;
}
}
}
for(j=0;j<=b;j++) v[a][j]=0.0;
return(0);
}