Rayleigh_Benard对流格子Boltzmann代码
// natural convection for closed squre cavity.cpp : Defines the entry point for the console application.
//coupled double distribution fountion
# include "stdafx.h"
# include
# include
# include
# include
# include
# include
# include
using namespace std;
//===============================================Constant declarations==================================================
const int Q = 9; //D2Q9 model
const int NX = 100; //lattice of the x direction
const int NY = 50; //lattice of the y direction
const double T_hot = 2.; //The bottom wall temperature
const double T_cold =1.; //The top wall temperature
double Gg[2] = { 0, -0.005 }; //acceleration of gravity
//================================================Variable declaration====================================================
//The discrete velocity direction of D2Q9 model
double e[Q][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 }, { -1, 0 }, { 0, -1 }, { 1, 1 }, { -1, 1 }, { -1, -1 }, { 1, -1 } };
//weight coefficient of D2Q9 model
double w[Q] = { 4.0 / 9, 1.0 / 9, 1.0 / 9, 1.0 / 9, 1.0 / 9, 1.0 / 36, 1.0 / 36, 1.0 / 36, 1.0 / 36 };
double rho[NX + 1][NY + 1], //density
u[NX + 1][NY + 1][2], //velocity of n+1 time
u0[NX + 1][NY + 1][2], //velocity of n time
f[NX + 1][NY + 1][Q], //velocity distribution funtion before evolution
F[NX + 1][NY + 1][Q], //velocity dietribution funtion after evolution
vel[NX + 1][NY + 1]; //velocity of (i.j)
double T[NX + 1][NY + 1], //temperture of n+1 time
g[NX + 1][NY + 1][Q], //temperture distribution funtion before evolution
G[NX + 1][NY + 1][Q]; //temperture distribution funtion after evolution
double wf, wg;
int i, j, k, ip, jp, n;
double LX,LY,c, Fa,Ra, Prandtl, Rg=8.3145, dx, dy, dt, rho0, P0, tau_f, tau_g, niu, error,beita,Nu,Tav;
//==================================================funtions of this code===================================================
void init();
double feq(int k, double rho, double u[2]);
double geq(int k, double T, double u[2]);
void evolution();
void output(int m);
void Error();
//=======================================================main loop============================================================
int main()
{
using namespace std;
init();
for (n = 0;; n++)
{
evolution()
;
if (n % 5000 == 0)
{
Error();
cout << "The " << n << " error " << setiosflags(ios::scientific) << error << " NU= " << Nu << " Tav= " << Tav << endl;
if (n >= 20000)
{
if (n %20000 == 0)
output(n);
if (error<1.0e-10)
break;
}
}
}
return 0;
}
//======================================================The initialization====================================================
void init()
{
dy = 1;// Lattice step of x direction
dt = dy; // time step
LY = double(NY)*dy;
LX = 2 * LY;
c = dy / dt; // The lattice velocity
rho0 = 1.0; // initial density
Ra = 5.0e3; // Rayleigh number
Prandtl = 0.71; // Prandtl number
niu = 1.e-5;
//Rg = 1 / (3 * (T_hot + T_cold) / 2);// kinematic viscosity
tau_f = 0.1*NY*pow(3 * Prandtl, 0.5) / pow(Ra, 0.5) + 0.5; //relaxation time of velocity field
tau_g = (tau_f-0.5) / Prandtl+0.5; //relaxation time of temperture field
wf = dt / (tau_f + 0.5*dt);
wg = dt / (tau_g + 0.5*dt);
beita = Ra*niu / (abs(Gg[1]) * (T_hot - T_cold)*NY*NY*NY*Prandtl);
for (i = 0; i <= NX; i++)
{
for (j = 0; j <= NY; j++)
{
u[i][j][0] = 0;
u[i][j][1] = 0;
rho[i][j] = rho0;
T[i][j] = (T_hot+T_cold)/2.;
for (k = 0; k
f[i][j][k] = feq(k, rho[i][j], u[i][j]); //initialization distribution function
g[i][j][k] = geq(k, T[i][j], u[i][j]); //velocity field and temperture
}
}
}
}
//The velocity equilibrium distribution function
double feq(int k, double rho, double u[2])
{
double eu, uv, feq;
eu = (e[k][0] * u[0] + e[k][1] * u[1]);
uv = (u[0] * u[0] + u[1] * u[1]);
feq = w[k] * rho*(1.0 + 3.0*eu + 4.5*eu*eu - 1.5*uv);
return feq;
}
//The tempreture equilibrium distribution function
double geq(int k, double T, double u[2])
{
double eu, uv, geq;
eu = (e[k][0] * u[0] + e[k][1] * u[1]);
uv = (u[0] * u[0] + u[1] * u[1]);
geq = w[k] * T*(1.0 + 3.0*eu + 4.5*eu*eu - 1.5*uv);
return geq;
}
//evolution :collision,streaming,macro,boundary;
void evolution()
{
for (i = 1; i
for (j = 1; j
for (k = 0; k
ip = i - int(e[k][0]); //collision and streaming
jp = j - int(e[k][1]); //
//evolution equation of velocity field
//Fa = (1-1/(2*tau_f))*w[k]*(3*(e[k][1]-u[ip][jp][1])+9*e[k][1]*u[ip][jp][1]*e[k][1])*(beita*(T[ip][jp]-(T_hot+T_cold)/2)*Gg[1]);
Fa = dt*Gg[1] * (e[k][1] - u[i
p][jp][1])*feq(k, rho[ip][jp], u[ip][jp]) /Rg/T[ip][jp];
F[i][j][k] = f[ip][jp][k] + (feq(k, rho[ip][jp], u[ip][jp]) - f[ip][jp][k]) / tau_f + dt*Fa;
//evolution equation of temperture field
G[i][j][k] = g[ip][jp][k] + (geq(k, T[ip][jp], u[ip][jp]) - g[ip][jp][k]) / tau_g;
}
}
}
//The microscopic distribution function to the macroscopic physical quantity
for (i = 1; i
for (j = 1; j
u0[i][j][0] = u[i][j][0];
u0[i][j][1] = u[i][j][1];
rho[i][j] = 0;
u[i][j][0] = 0;
u[i][j][1] = 0;
T[i][j] = 0;
for (k = 0; k
f[i][j][k] = F[i][j][k];
rho[i][j] += f[i][j][k];
g[i][j][k] = G[i][j][k];
u[i][j][0] += e[k][0] * f[i][j][k];
u[i][j][1] += e[k][1] * f[i][j][k];
T[i][j] += g[i][j][k];
}
u[i][j][0] /= rho[i][j];
u[i][j][1] /= rho[i][j];
vel[i][j] = sqrt(u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1]);
}
}
Nu = 0.0;
//computational Nu
for (i = 0; i <= NX; i++)
{
Nu=Nu+(T[i][0] - T[i][1]);
}
Tav = 0.0;
for (i = 0; i <= NX; i++)
{
for (j = 0; j <= NY; j++)
{
Tav=Tav+(T[i][j] - 1);
}
}
Tav = Tav / ((NX + 1)*(NY + 1));
//---------------------------------------------------------boundary treatment---------------------------------------------------
//The non-equilibrium extrapolation scheme for boundary treatment and periodic boundary
//left and right periodic
for (j = 0; j <= NY; j++)
{
for (k = 0; k < Q; k++)
{
rho[NX][j] = rho[NX - 1][j];
T[NX][j] = T[NX - 1][j];
if (k == 1 || k == 5 || k == 8)
{
f[0][j][k] = f[NX][j][k];
g[0][j][k] = g[NX][j][k];
}
else
{
f[NX][j][k] = feq(k, rho[NX][j], u[NX][j]) + (f[NX - 1][j][k] - feq(k, rho[NX - 1][j], u[NX - 1][j]));
g[NX][j][k] = geq(k, T[NX][j], u[NX][j]) + (g[NX - 1][j][k] - geq(k, T[NX - 1][j], u[NX - 1][j]));
}
rho[0][j] = rho[1][j];
T[0][j] = T[1][j];
if (k == 3 || k == 6 || k == 7)
{
f[NX+1][j][k] = f[1][j][k];
g[NX+1][j][k] = g[1][j][k];
}
else
{
f[0][j][k] = feq(k, rho[0][j], u[0][j]) + (f[1][j][k] - feq(k, rho[1][j], u[1][j]));
g[0][j][k] = geq(k, T[0][j], u[0][j]) + (g[1][j][k] - geq(k, T[1][j], u[1][j]));
}
}
}
//top and bottom non-equilibrium extrapolation scheme
for (i = 0; i <= NX; i++)
{
T[i][0] = T_hot;
rho[i][0] = rho[i][1];
T[i][NY] = T_cold;
rho[i][NY] = rho[i][NY - 1];
for (k = 0; k < Q; k++)
{
f[i][0][k] = feq(k, rho[i][0], u[i][0]) + f[i][1][k] - feq(k, rho[i][1], u[i][1]);
f[i][NY][k] = feq(k, rho[i][NY], u[i][NY]) + f[i][NY - 1][k] - feq(k, rho[i][NY - 1], u[i][NY - 1]);
g[i][0][k] = geq(k, T[i][0], u[i][0]) + g[i][1][k] - geq(k, T[i][1], u[i][1]);
g[i][NY][k] = geq(k, T[i][NY], u[i][NY]) + g[i][NY - 1][k] - geq(k, T[i][NY - 1], u[i][NY - 1]);
}
}
}
//=====================================================o
uput data of Physical field================================================
void output(int m)
{
ostringstream name;
name << "Rayleigh Benard convection_"<< m << ".dat";
ofstream out(name.str().c_str());
out << "Titel=\"LBM convection\"\n"
<< "VARIABLES=\"X\",\"Y\",\"U\",\"V\",\"vel\",\"rho\",\"T\"\n" << "ZONE T=\"BOX\",I="
<< NX + 1 << ",J=" << NY + 1 << ",F=POINT" << endl;
for (j = 0; j <= NY; j++)
for (i = 0; i <= NX;i++)
{
out <
<< endl;
}
}
//=====================================================The max relative error of uv====================================================
void Error()
{
double temp1, temp2;
temp1 = 0;
temp2 = 0;
for (i = 0; i <= NX; i++)
{
for (j = 0; j <= NY; j++)
{
temp1 += ((u[i][j][0] - u0[i][j][0])*(u[i][j][0] - u0[i][j][0])
+ (u[i][j][1] - u0[i][j][1])*(u[i][j][1] - u0[i][j][1]));
temp2 += (u[i][j][0] * u[i][j][0] + u[i][j][1] * u[i][j][1]);
}
}
temp1 = sqrt(temp1);
temp2 = sqrt(temp2);
error = temp1 / (temp2 + 1e-30);
}