📄 d2q9_c.mht
字号:
From: <óé Microsoft Internet Explorer 5 ±￡′?>
Subject:
Date: Fri, 30 Jun 2006 15:45:07 +0800
MIME-Version: 1.0
Content-Type: text/html;
charset="iso-8859-1"
Content-Transfer-Encoding: quoted-printable
Content-Location: http://www.physics.ndsu.nodak.edu/wagner/source/D2Q9.c
X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2900.2869
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<HTML><HEAD>
<META http-equiv=3DContent-Type content=3D"text/html; =
charset=3Diso-8859-1">
<META content=3D"MSHTML 6.00.2900.2912" name=3DGENERATOR></HEAD>
<BODY><PRE>/* This is a standard two dimensional lattice Boltzmann =
program with
periodic boundary conditions that should perform reasonably =
efficiently.
It has been written for teaching purposes by Alexander Wagner
at NDSU (March 2003).=20
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h> /* for the sleep function. this may not be =
standard */
#include <math.h>
#include <mygraph.h>
#define xdim 31
#define ydim 21
#define DEBUG
double f0[xdim][ydim],
f1[xdim][ydim],f2[xdim][ydim],f3[xdim][ydim],f4[xdim][ydim],
f5[xdim][ydim],f6[xdim][ydim],f7[xdim][ydim],f8[xdim][ydim];
#ifdef DEBUG
double b1[xdim][ydim],b3[xdim][ydim],b5[xdim][ydim],b6[xdim][ydim],
bb[xdim][ydim];
#endif
double omega=3D1;
/* Some constants that appear often and don't need to be calculated
over and over again. So we calculate them just once here. */
double
fourOnine=3D4.0/9.0,
oneOnine=3D1.0/9.0,
oneOthirtysix=3D1.0/36.0;
/* Some variables for the suspended particle */
typedef double *doublep;
typedef doublep pair[2];
pair *Bf1=3DNULL,*Bf3=3DNULL,*Bf5=3DNULL,*Bf6=3DNULL;
int Bf1c,Bf1m=3D0,Bf3c,Bf3m=3D0,Bf5c,Bf5m=3D0,Bf6c,Bf6m=3D0;
double =
Cx=3D15,Cy=3D10,Mx=3D0,My=3D0,Zxx=3D0,Zxy=3D0,Zyy=3D0,Ux=3D0,Uy=3D0,R=3D5=
,M=3D100;
int Repeat=3D1,iterations=3D0,FrequencyMeasure=3D100,Graphics=3D1;
int Pause=3D1,sstep=3D0,done=3D0;
double Amp=3D0.001,Amp2=3D0.1,Width=3D10;
/* Some special data types that are required to view the graphics */
int Xdim=3Dxdim,Ydim=3Dydim,noreq;
double rho[xdim][ydim];
int rhoreq=3D0;
double u[xdim][ydim][2];
int ureq=3D0;
void BCmem(pair **p, int c, int *cm){
if (c<*cm-2) return;
*cm+=3D100;
*p=3D(pair *) realloc(*p,*cm*sizeof(pair));
}
void circleBC(){
/* Sets up the links that are cut by a circular object. */
int x,y,Px,Py,Pxp,Pyp,R2;
double dx,dy;
#ifdef DEBUG
for (x=3D0;x<xdim;x++) for (y=3D0;y<ydim;y++)
b1[x][y]=3Db3[x][y]=3Db5[x][y]=3Db6[x][y]=3Dbb[x][y]=3D0;
#endif
Bf1c=3DBf3c=3DBf5c=3DBf6c=3D0;
R2=3DR*R;
for (x=3DCx-R-2;x<Cx+R+2;x++){
Px=3D(x+xdim)%xdim;
Pxp=3D(x+1+xdim)%xdim;
dx=3Dx-Cx;
for (y=3DCy-R-2;y<Cy+R+2;y++){
Py=3D(y+ydim)%ydim;
Pyp=3D(y+1+ydim)%ydim;
dy=3Dy-Cy;
#ifdef DEBUG
if ((dx*dx+dy*dy-R2)<=3D0) bb[Px][Py]=3D1;else bb[Px][Py]=3D-1;
#endif
if ((dx*dx+dy*dy-R2)*((dx+1)*(dx+1)+dy*dy-R2)<=3D0){
BCmem(&Bf1,Bf1c,&Bf1m);
Bf1[Bf1c][0]=3D&(f2[Px][Py]);
Bf1[Bf1c++][1]=3D&(f1[(Pxp)][Py]);
#ifdef DEBUG
b1[Px][Py]=3D1;
} else {
b1[Px][Py]=3D-1;
#endif=20
}
if ((dx*dx+dy*dy-R2)*((dx)*(dx)+(dy+1)*(dy+1)-R2)<=3D0){
BCmem(&Bf3,Bf3c,&Bf3m);
Bf3[Bf3c][0]=3D&(f4[Px][Py]);
Bf3[Bf3c++][1]=3D&(f3[Px][Pyp]);
#ifdef DEBUG
b3[Px][Py]=3D1;
} else {
b3[Px][Py]=3D-1;
#endif=20
}
if ((dx*dx+dy*dy-R2)*((dx+1)*(dx+1)+(dy+1)*(dy+1)-R2)<=3D0){
BCmem(&Bf5,Bf5c,&Bf5m);
Bf5[Bf5c][0]=3D&(f7[Px][Py]);
Bf5[Bf5c++][1]=3D&(f5[Pxp][Pyp]);
#ifdef DEBUG
b5[Px][Py]=3D1;
} else {
b5[Px][Py]=3D-1;
#endif=20
}
if =
(((dx)*(dx)+(dy+1)*(dy+1)-R2)*((dx+1)*(dx+1)+(dy)*(dy)-R2)<=3D0){
BCmem(&Bf6,Bf6c,&Bf6m);
Bf6[Bf6c][0]=3D&(f8[Pxp][Py]);
Bf6[Bf6c++][1]=3D&(f6[Px][Pyp]);
#ifdef DEBUG
b6[Px][Py]=3D1;
} else {
b6[Px][Py]=3D-1;
#endif=20
}
}
}
}
void initParticle(){
Cx=3D10;Cy=3D10;Mx=3D0;My=3D0;Ux=3D0;Uy=3D0;R=3D5;M=3D100;
}
void init(){
int x,y;
double n,ux,uy,uxx,uyy,uxy,usq;
printf("initialize\n");
iterations=3D0;
for (x=3D0;x<xdim;x++)
for (y=3D0;y<ydim;y++){
/* here we can define the macroscopic properties of the=20
initial condition. */
n=3D1+Amp2*exp(-(pow(x-xdim/2,2)+pow(y-ydim/2,2))/Width);
ux=3D0/*.05*sin((y-ydim/2)*2*M_PI/xdim)*/;
uy=3D0;
/*n=3D1;ux=3D0;uy=3D0;*/
/* The following code initializes the f to be the local =
equilibrium
values associated with the density and velocity defined above.*/
uxx=3Dux*ux;
uyy=3Duy*uy;
uxy=3D2*ux*uy;
usq=3Duxx+uyy;
f0[x][y]=3DfourOnine*n*(1-1.5*usq);
f1[x][y]=3DoneOnine*n*(1+3*ux+4.5*uxx-1.5*usq);
f2[x][y]=3DoneOnine*n*(1-3*ux+4.5*uxx-1.5*usq);
f3[x][y]=3DoneOnine*n*(1+3*uy+4.5*uyy-1.5*usq);
f4[x][y]=3DoneOnine*n*(1-3*uy+4.5*uyy-1.5*usq);
=
f5[x][y]=3DoneOthirtysix*n*(1+3*(ux+uy)+4.5*(uxx+uxy+uyy)-1.5*usq);
=
f6[x][y]=3DoneOthirtysix*n*(1+3*(-ux+uy)+4.5*(uxx-uxy+uyy)-1.5*usq);
=
f7[x][y]=3DoneOthirtysix*n*(1+3*(-ux-uy)+4.5*(uxx+uxy+uyy)-1.5*usq);
=
f8[x][y]=3DoneOthirtysix*n*(1+3*(ux-uy)+4.5*(uxx-uxy+uyy)-1.5*usq);
=20
}
}
void TotMomentum(){
int x,y;
double Momx,Momy;
Momx=3DMomy=3D0;
for (x=3D0;x<xdim;x++)
for (y=3D0;y<ydim;y++){
Momx+=3Df1[x][y]-f2[x][y]+f5[x][y]-f6[x][y]-f7[x][y]+f8[x][y];
Momy+=3Df3[x][y]-f4[x][y]+f5[x][y]+f6[x][y]-f7[x][y]-f8[x][y];
}
printf("MomF =3D (%e,%e), MomP =3D (%e,%e), MomT =3D (%e,%e)\n",
Momx,Momy,M*Ux,M*Uy,Momx+M*Ux,Momy+M*Uy);
}
void iterationColloid(){
#define alp 1 /* the degree of implicitness */
double zxx,zxy,zyy; /* The inverse Z tensor */
Cx+=3DUx+xdim;
Cx=3Dfmod(Cx,xdim);
Cy+=3DUy+ydim;
Cy=3Dfmod(Cy,ydim);
=20
zxx=3D(M+alp*Zyy)/((M+alp*Zxx)*(M+alp*Zyy)+alp*alp*Zxy*Zxy);
zxy=3D alp*Zxy /((M+alp*Zxx)*(M+alp*Zyy)+alp*alp*Zxy*Zxy);
zyy=3D(M+alp*Zxx)/((M+alp*Zxx)*(M+alp*Zyy)+alp*alp*Zxy*Zxy);
Ux+=3Dzxx*(Mx-Zxx*Ux-Zxy*Uy)+zxy*(My-Zxy*Ux-Zyy*Uy);
Uy+=3Dzxy*(Mx-Zxx*Ux-Zxy*Uy)+zyy*(My-Zxy*Ux-Zyy*Uy);
#undef alp
}
void iteration(){
int x,y,i;
register double tmp1,tmp2;
register double n,ux,uy,uxx,uyy,uxy,usq,Fx,Fy,Fxx,Fyy,Fxy,Fsq;
double f1y[ydim],f2y[ydim],f5y[ydim],f6y[ydim],f7y[ydim],f8y[ydim];
double f3x[xdim],f4x[xdim],f5x[xdim],f6x[xdim],f7x[xdim],f8x[xdim];
/* first we perform the collision step */
=20
for (x=3D0;x<xdim;x++)
for (y=3D0;y<ydim;y++){
n=3Df0[x][y]+f1[x][y]+f2[x][y]+f3[x][y]+f4[x][y]
+f5[x][y]+f6[x][y]+f7[x][y]+f8[x][y];
ux=3Df1[x][y]-f2[x][y]+f5[x][y]-f6[x][y]-f7[x][y]+f8[x][y];
uy=3Df3[x][y]-f4[x][y]+f5[x][y]+f6[x][y]-f7[x][y]-f8[x][y];
ux/=3Dn;
uy/=3Dn;
uxx=3Dux*ux;
uyy=3Duy*uy;
uxy=3D2*ux*uy;
usq=3Duxx+uyy;
/* We now need to implement any forcing terms. The term included
here is just an example. You should not calculate a constant
force term in the interation routine, but outside where it only
gets calculated once.*/
Fx=3DAmp*sin(y*2*M_PI/ydim);
Fy=3D0;
Fxx=3D2*n*Fx*ux;
Fyy=3D2*n*Fy*uy;
Fxy=3D2*n*(Fx*uy+Fy*ux);
Fsq=3DFxx+Fyy;
Fx*=3Dn;
Fy*=3Dn;
f0[x][y]+=3Domega*(fourOnine*n*(1-1.5*usq)-f0[x][y])
-fourOnine*1.5*Fsq;
f1[x][y]+=3Domega*(oneOnine*n*(1+3*ux+4.5*uxx -1.5*usq)-f1[x][y])
+oneOnine*(3*Fx+4.5*Fxx-1.5*Fsq);
f2[x][y]+=3Domega*(oneOnine*n*(1-3*ux+4.5*uxx -1.5*usq)-f2[x][y])
+oneOnine*(-3*Fx+4.5*Fxx-1.5*Fsq);
f3[x][y]+=3Domega*(oneOnine*n*(1+3*uy+4.5*uyy -1.5*usq)-f3[x][y])
+oneOnine*(3*Fy+4.5*Fyy-1.5*Fsq);
f4[x][y]+=3Domega*(oneOnine*n*(1-3*uy+4.5*uyy -1.5*usq)-f4[x][y])
+oneOnine*(-3*Fy+4.5*Fyy-1.5*Fsq);
f5[x][y]+=3Domega*(oneOthirtysix*n*(1+3*(ux+uy)+4.5*(uxx+uxy+uyy)
-1.5*usq)-f5[x][y])
+oneOthirtysix*(3*(Fx+Fy)+4.5*(Fxx+Fxy+Fyy)-1.5*Fsq);
f6[x][y]+=3Domega*(oneOthirtysix*n*(1+3*(-ux+uy)+4.5*(uxx-uxy+uyy)
-1.5*usq)-f6[x][y])
+oneOthirtysix*(3*(-Fx+Fy)+4.5*(Fxx-Fxy+Fyy)-1.5*Fsq);
f7[x][y]+=3Domega*(oneOthirtysix*n*(1+3*(-ux-uy)+4.5*(uxx+uxy+uyy)
-1.5*usq)-f7[x][y])
+oneOthirtysix*(3*(-Fx-Fy)+4.5*(Fxx+Fxy+Fyy)-1.5*Fsq);
f8[x][y]+=3Domega*(oneOthirtysix*n*(1+3*(ux-uy)+4.5*(uxx-uxy+uyy)
-1.5*usq)-f8[x][y])
+oneOthirtysix*(3*(Fx-Fy)+4.5*(Fxx-Fxy+Fyy)-1.5*Fsq);
}
=20
/* now we need to move the densities according to their velocities
we are using periodic boundary conditions */
/* since we are only using one lattice, we need to save some data on=20
the boundaries so that it does not get overwritten */
for (y=3D0;y<ydim;y++){
f1y[y]=3Df1[xdim-1][y];
f2y[y]=3Df2[0][y];
f5y[y]=3Df5[xdim-1][y];
f6y[y]=3Df6[0][y];
f7y[y]=3Df7[0][y];
f8y[y]=3Df8[xdim-1][y];
}
for (x=3D0;x<xdim;x++){
f3x[x]=3Df3[x][ydim-1];
f4x[x]=3Df4[x][0];
f5x[x]=3Df5[x][ydim-1];
f6x[x]=3Df6[x][ydim-1];
f7x[x]=3Df7[x][0];
f8x[x]=3Df8[x][0];
}
=20
/* Now we can move the densities along the lattice.
You can also do this using loops, but that is actually
more complicated */
memmove(&f1[1][0],&f1[0][0],(xdim-1)*ydim*sizeof(double));
memmove(&f2[0][0],&f2[1][0],(xdim-1)*ydim*sizeof(double));
memmove(&f3[0][1],&f3[0][0],(xdim*ydim-1)*sizeof(double));
memmove(&f4[0][0],&f4[0][1],(xdim*ydim-1)*sizeof(double));
memmove(&f5[1][1],&f5[0][0],((xdim-1)*ydim-1)*sizeof(double));
memmove(&f7[0][0],&f7[1][1],((xdim-1)*ydim-1)*sizeof(double));
memmove(&f6[0][1],&f6[1][0],((xdim-1)*ydim)*sizeof(double));
memmove(&f8[1][0],&f8[0][1],((xdim-1)*ydim)*sizeof(double));
/* Now we need to fix the boundaries that have not yet been correctly
updated */
for (y=3D0;y<ydim;y++){
f1[0][y] =3Df1y[y];
f2[xdim-1][y] =3Df2y[y];
f5[0][(y+1)%ydim] =3Df5y[y];
f6[xdim-1][(y+1)%ydim]=3Df6y[y];
f7[xdim-1][y] =3Df7y[(y+1)%ydim];
f8[0][y] =3Df8y[(y+1)%ydim];
}
for (x=3D0;x<xdim;x++){
f3[x][0] =3Df3x[x];
f4[x][ydim-1] =3Df4x[x];
f5[(x+1)%xdim][0] =3Df5x[x];
f6[x][0] =3Df6x[(x+1)%xdim];
f7[x][ydim-1] =3Df7x[(x+1)%xdim];
f8[(x+1)%xdim][ydim-1]=3Df8x[x];
}
/* Objects in flow */
Mx=3DMy=3DZxx=3DZxy=3DZyy=3D0;
/* Really I am missing a factor of n here for the velocity corrections =
*/
for (i=3D0;i<Bf1c;i++){
tmp1=3D*(Bf1[i][0]);
tmp2=3D*(Bf1[i][1]);
*(Bf1[i][0])=3Dtmp2-2./3.*Ux;
*(Bf1[i][1])=3Dtmp1+2./3.*Ux;
Mx+=3D2*(tmp2-tmp1);
Zxx+=3D 2*2./3.;
}
for (i=3D0;i<Bf3c;i++){
tmp1=3D*(Bf3[i][0]);
tmp2=3D*(Bf3[i][1]);
*(Bf3[i][0])=3Dtmp2-2./3.*Uy;
*(Bf3[i][1])=3Dtmp1+2./3.*Uy;
My+=3D2*(tmp2-tmp1);
Zyy+=3D 2*2./3.;
}
for (i=3D0;i<Bf5c;i++){
tmp1=3D*(Bf5[i][0]);
tmp2=3D*(Bf5[i][1]);
*(Bf5[i][0])=3Dtmp2-1./6.*(Ux+Uy);
*(Bf5[i][1])=3Dtmp1+1./6.*(Ux+Uy);
Mx+=3D2*(tmp2-tmp1);
My+=3D2*(tmp2-tmp1);
Zxx+=3D 2*1./6.;
Zxy+=3D 2*1./6.;
Zyy+=3D 2*1./6.;
}
for (i=3D0;i<Bf6c;i++){
tmp1=3D*(Bf6[i][0]);
tmp2=3D*(Bf6[i][1]);
*(Bf6[i][0])=3Dtmp2-1./6.*(-Ux+Uy);
*(Bf6[i][1])=3Dtmp1+1./6.*(-Ux+Uy);
Mx+=3D-2*(tmp2-tmp1);
My+=3D 2*(tmp2-tmp1);
Zxx+=3D 2*1./6.;
Zxy+=3D-2*1./6.;
Zyy+=3D 2*1./6.;
}
}
void analysis(int iterations){
if (FrequencyMeasure%iterations=3D=3D0){
=20
}
}
void GetGraphics(){
double n;
int x,y;
if (rhoreq){
rhoreq=3D0;
for (x=3D0;x<xdim;x++)
for (y=3D0;y<ydim;y++){
rho[x][y]=3Df0[x][y]+f1[x][y]+f2[x][y]+f3[x][y]+f4[x][y]
+f5[x][y]+f6[x][y]+f7[x][y]+f8[x][y];
}
}
if (ureq){
ureq=3D0;
for (x=3D0;x<xdim;x++)
for (y=3D0;y<ydim;y++){
n=3Df0[x][y]+f1[x][y]+f2[x][y]+f3[x][y]+f4[x][y]
+f5[x][y]+f6[x][y]+f7[x][y]+f8[x][y];
u[x][y][0]=3Df1[x][y]-f2[x][y]+f5[x][y]-f6[x][y]-f7[x][y]+f8[x][y];
u[x][y][1]=3Df3[x][y]-f4[x][y]+f5[x][y]+f6[x][y]-f7[x][y]-f8[x][y];
u[x][y][0]/=3Dn;
u[x][y][1]/=3Dn;
}
}=20
}
void GUI(){
=
DefineGraphNxN_R("rho",&rho[0][0],&Xdim,&Ydim,&rhoreq);
#ifdef DEBUG
DefineGraphNxN_R("bb",&bb[0][0],&Xdim,&Ydim,&noreq);
DefineGraphNxN_R("b1",&b1[0][0],&Xdim,&Ydim,&noreq);
DefineGraphNxN_R("b3",&b3[0][0],&Xdim,&Ydim,&noreq);
DefineGraphNxN_R("b5",&b5[0][0],&Xdim,&Ydim,&noreq);
DefineGraphNxN_R("b6",&b6[0][0],&Xdim,&Ydim,&noreq);
#endif
DefineGraphNxN_RxR("u",&u[0][0][0],&Xdim,&Ydim,&ureq);
StartMenu("Lattice Boltzmann",1);
DefineDouble("1/tau",&omega);
DefineInt("Measurement freq.",&FrequencyMeasure);
StartMenu("Reinitialize",0);
DefineDouble("Amplitude",&Amp2);
DefineDouble("Width",&Width);
DefineFunction("reinitialize",&init);
DefineFunction("Particle init",&initParticle);
EndMenu();
DefineDouble("Velocity amplitude",&Amp);
StartMenu("Particle",0);
DefineDouble("R",&R);
DefineDouble("M",&M);
DefineDouble("Ux",&Ux);
DefineDouble("Uy",&Uy);
DefineDouble("Cx",&Cx);
DefineDouble("Cy",&Cy);
DefineDouble("Mx",&Mx);
DefineDouble("My",&My);
EndMenu();
DefineGraph(contour2d_,"Density&Vector plots");
DefineInt("Repeat",&Repeat);
DefineBool("Pause",&Pause);
DefineBool("Single Step",&sstep);
DefineBool("Done",&done);
EndMenu();
}
int main(){
int i,newdata;
=20
if (Graphics) GUI();
init();
while (!done){
if (Graphics){
Events(newdata);
GetGraphics();
DrawGraphs();
} else {done=3D1;Pause=3D0;}
if (!Pause||sstep){
sstep=3D0;
newdata=3D1;
for (i=3D0;i<Repeat;i++){
iterations++;
circleBC();
iteration();
iterationColloid();
TotMomentum();
analysis(iterations);
}
} else sleep(1);
}
return 0;
}
</PRE></BODY></HTML>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -