Recently, I learn the case of flow around circular cylinder. But when I calculate the drag coefficient of the cylinder, there is large deviation with reference.
Please help me
My code are as follw:
[code=“cpp”]
#include
#include
#include
#include
#include
#include
#include
using namespace std;
const int Q = 9;
const int NX = 1000;
const int NY = 160;
const double dx = 0.5;
const double dy = 0.5;
const double dt = dx;
const double U = 0.1;
const double D = 10/dx;
const double cx = 0.25NX;
const double cy = 0.5NY;
const double cr = D / 2.0;
int e[Q][2] = { { 0,0 },{ 1,0 },{ 0,1 },{ -1,0 },{ 0,-1 },{ 1,1 },{ -1,1 },{ -1,-1 },{ 1,-1 } };
int e1[Q] = { 0,3,4,1,2,7,8,5,6 };
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], u[NX + 1][NY + 1][2], u0[NX + 1][NY + 1][2], f[NX + 1][NY + 1][Q], F[NX + 1][NY + 1][Q];
int i, j, k, ip, jp, n;
double c, CD, CL, Lx, Ly, Re, niu, force1, force2, forcex, forcey, P, rho0, tau_f, error, Umax, q, a1, b1, c1, x1, x2, L;
void init();
double feq(int k, double rho, double u[2]);
void evolution();
void output(int m);
void Error();
double distance(double i, double j);
int main()
{
using namespace std;
init();
for (n = 1; ; n++)
{
evolution();
if (n % 1000 == 0)
{
Error();
cout << "The" << ends << n << "th computation result:" << endl << "The u,v of point (NX/2,NY/2)is:" << setprecision(6) << u[NX / 2][NY / 2][0] << "," << u[NX / 2][NY / 2][1] << endl;
cout << "The max relative error of uv is:"
<< setiosflags(ios::scientific) << error << endl;
cout << "CD=2*forcex/(rho0*Umax*Umax*D)=:" << 2 * forcex / (rho0*U*U *D) << " " << "CL=2*forcey/(rho0*Umax*Umax*D)=" << 2 * forcey / (rho0*U*U *D) << endl;
if (n >= 1000)
{
if (n % 1000 == 0) output(n);
if (error<1.0e-8) break;
}
}
}
return 0;
}
void init()
{
Lx = dx*double(NX);
Ly = dy*double(NY);
c = dx / dt;//1.0
rho0 = 1.0;
Re = 20;
L = D*D / 4.0;
niu = U*D *dx/ Re;
tau_f = 3 * niu / dt + 0.5;
std::cout << "tau_f=" << tau_f << ends << "niu=" << niu << endl;
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;
u[0][j][0] = U;
u[0][0][0] = 0;
u[0][NY][0] = 0;
for (k = 0; k<Q; k++)
{
f[i][j][k] = feq(k, rho[i][j], u[i][j]);
}
}
}
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 + 3.0eu + 4.5eueu - 1.5uv);
return feq;
}
double distance(double i, double j)
{
double as, bs, cs;
as = i;
bs = j;
cs = (i - cx)(i - cx) + (j - cy)(j - cy);
return cs;
}
void evolution()
{
for (i = 0; i <= NX; i++)
for (j = 0; j <= NY; j++)
for (k = 0; k < Q; k++)
{
F[i][j][k] = f[i][j][k] + (feq(k, rho[i][j], u[i][j]) - f[i][j][k]) / tau_f;
}
for (i = 0; i <= NX; i++)
for (j = 0; j <= NY; j++)
for (k = 0; k<Q; k++)
{
ip = i - e[k][0];
jp = j - e[k][1];
if (ip >= 0 && ip <= NX&&jp >= 0 && jp <= NY)
f[i][j][k] = F[ip][jp][k];
}
for (i = cx-cr; i <= cx+cr; i++)
for (j = cy-cr; j <= cy+cr; j++)
for (k = 1; k<Q; k++)
{
if (distance(double(i), double(j)) <=L&&distance(double(i + e[k][0]), double(j + e[k][1])) > L)
{
a1 = e[k][0] * e[k][0] + e[k][1] * e[k][1];
b1 = 2 * (e[k][0] * (i - cx) + e[k][1] * (j - cy));
c1 = (i - cx)*(i - cx) + (j - cy)*(j - cy) - L;
x1 = (-b1 + sqrt(b1*b1 - 4 * a1*c1)) / (2 * a1);
x2 = (-b1 - sqrt(b1*b1 - 4 * a1*c1)) / (2 * a1);
if (x1 > 0)
{
q = 1 - x1;
}
else
{
q = 1 - x2;
}
if (q >= 0.5)
{
f[i + e[k][0]][j + e[k][1]][k] = 1 / (q*(1 + 2 * q))*F[i + e[k][0]][j + e[k][1]][e1[k]] + (2 * q - 1) / q*f[i + 2 * e[k][0]][j + 2 * e[k][1]][k] - (2 * q - 1) / (2 * q + 1)*f[i + 3 * e[k][0]][j + 3 * e[k][1]][k];
}
else if (q<0.5)
{
f[i + e[k][0]][j + e[k][1]][k] = q*(1 + 2 * q)*F[i + e[k][0]][j + e[k][1]][e1[k]] + (1 - 4 * q*q)*F[i + 2 * e[k][0]][j + 2 * e[k][1]][e1[k]] - q*(1 - 2 * q)*F[i + 3 * e[k][0]][j + 3 * e[k][1]][e1[k]];
}
}
}
for (i = 1; i<NX; i++)
for (j = 1; j<NY; j++)
{
u0[i][j][0] = u[i][j][0];
u0[i][j][1] = u[i][j][1];
rho[i][j] = f[i][j][0] + f[i][j][1] + f[i][j][2] + f[i][j][3] + f[i][j][4] + f[i][j][5] + f[i][j][6] + f[i][j][7] + f[i][j][8];
u[i][j][0] = (f[i][j][1] + f[i][j][5] + f[i][j][8] - f[i][j][3] - f[i][j][6] - f[i][j][7]) / rho[i][j];
u[i][j][1] = (f[i][j][2] + f[i][j][5] + f[i][j][6] - f[i][j][4] - f[i][j][7] - f[i][j][8]) / rho[i][j];
if (distance(double(i), double(j)) <= L)
{
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
u0[i][j][0] = 0.0;
u0[i][j][1] = 0.0;
}
}
forcex = 0;
forcey = 0;
for (i = cx - cr; i <= cx + cr; i++)
for (j = cy - cr; j <= cy + cr; j++)
for (k = 1; k < Q; k++)
{
if (distance(double(i), double(j)) <= L&&distance(double(i + e[k][0]), double(j + e[k][1])) > L)
{
force1 = e[e1[k]][0] * f[i + e[k][0]][j + e[k][1]][e1[k]] - e[k][0] * f[i][j][k];
forcex += force1;
force2 = e[e1[k]][1] * f[i + e[k][0]][j + e[k][1]][e1[k]] - e[k][1] * f[i][j][k];
forcey += force2;
}
}
i = 0;
for (j = 1; j <= NY - 1; j++)
for (k = 0; k<Q; k++)
{
u[i][j][0] = U;
u[i][j][1] = 0.0;
rho[i][j] = 1.0 / (1.0 - u[i][j][0])*(f[i][j][0] + f[i][j][2] + f[i][j][4] + 2.0*(f[i][j][3] + f[i][j][6] + f[i][j][7]));
f[i][j][1] = f[i][j][3] + 2.0 / 3.0*rho[i][j] * u[i][j][0];
f[i][j][5] = f[i][j][7] + 1.0 / 6.0*rho[i][j] * u[i][j][0] + 0.5*(f[i][j][4] - f[i][j][2]) + 1.0 / 2.0*rho[i][j] * u[i][j][1];
f[i][j][8] = f[i][j][6] + 1.0 / 6.0*rho[i][j] * u[i][j][0] - 0.5*(f[i][j][4] - f[i][j][2]) - 1.0 / 2.0*rho[i][j] * u[i][j][1];
}
// Outlet
i = NX;
for (j = 1; j <= NY - 1; j++)
for (k = 0; k<Q; k++)
{
u[i][j][0] = u[i - 1][j][0];
u[i][j][1] = u[i - 1][j][1];
f[i][j][3] = feq(3, rho[i][j], u[i][j]);
f[i][j][6] = feq(6, rho[i][j], u[i][j]);
f[i][j][7] = feq(7, rho[i][j], u[i][j]);
}
// Bottom
j = 0;
for (i = 0; i <= NX; i++)
for (k = 0; k<Q; k++)
{
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
rho[i][j] = 1.0 / (1 - u[i][j][1])*(f[i][j][0] + f[i][j][1] + f[i][j][3] + 2 * (f[i][j][4] + f[i][j][7] + f[i][j][8]));
f[i][j][2] = f[i][j][4] + 2.0 / 3.0*rho[i][j] * u[i][j][1];
f[i][j][5] = f[i][j][7] + 1.0 / 6.0*rho[i][j] * u[i][j][1] - 0.5*(f[i][j][1] - f[i][j][3]) + 1.0 / 2.0*rho[i][j] * u[i][j][0];
f[i][j][6] = f[i][j][8] + 1.0 / 6.0*rho[i][j] * u[i][j][1] + 0.5*(f[i][j][1] - f[i][j][3]) - 1.0 / 2.0*rho[i][j] * u[i][j][0];
}
// Top
j = NY;
for (i = 0; i <= NX; i++)
for (k = 0; k<Q; k++)
{
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
rho[i][j] = 1.0 / (1 + u[i][j][1])*(f[i][j][0] + f[i][j][1] + f[i][j][3] + 2.0*(f[i][j][2] + f[i][j][5] + f[i][j][6]));
f[i][j][4] = f[i][j][2] - 2.0 / 3.0*rho[i][j] * u[i][j][1];
f[i][j][7] = f[i][j][5] - 1.0 / 6.0*rho[i][j] * u[i][j][1] + 0.5*(f[i][j][1] - f[i][j][3]) - 1.0 / 2.0*rho[i][j] * u[i][j][0];
f[i][j][8] = f[i][j][6] - 1.0 / 6.0*rho[i][j] * u[i][j][1] - 0.5*(f[i][j][1] - f[i][j][3]) + 1.0 / 2.0*rho[i][j] * u[i][j][0];
}
// Corner Nodes
// Bottom Left
i = 0; j = 0;
rho[i][j] = rho[i][j + 1];
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
f[i][j][1] = f[i][j][3];
f[i][j][2] = f[i][j][4];
f[i][j][5] = f[i][j][7];
f[i][j][6] = 0.5*(rho[i][j] - f[i][j][1] - f[i][j][2] - f[i][j][3] - f[i][j][4] - f[i][j][5] - f[i][j][0] - f[i][j][7]);
f[i][j][8] = f[i][j][6];
// Top Left
i = 0; j = NY;
rho[i][j] = rho[i][j - 1];
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
f[i][j][1] = f[i][j][3];
f[i][j][4] = f[i][j][2];
f[i][j][8] = f[i][j][6];
f[i][j][5] = 0.5*(rho[i][j] - f[i][j][1] - f[i][j][2] - f[i][j][3] - f[i][j][4] - f[i][j][6] - f[i][j][0] - f[i][j][8]);
f[i][j][7] = f[i][j][5];
// Bottom Right
i = NX; j = 0;
rho[i][j] = rho[i][j + 1];
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
f[i][j][2] = f[i][j][4];
f[i][j][3] = f[i][j][1];
f[i][j][6] = f[i][j][8];
f[i][j][5] = 0.5*(rho[i][j] - f[i][j][1] - f[i][j][2] - f[i][j][3] - f[i][j][4] - f[i][j][6] - f[i][j][0] - f[i][j][8]);
f[i][j][7] = f[i][j][5];
// Top Right
i = NX; j = NY;
rho[i][j] = rho[i][j - 1];
u[i][j][0] = 0.0;
u[i][j][1] = 0.0;
f[i][j][3] = f[i][j][1];
f[i][j][4] = f[i][j][2];
f[i][j][7] = f[i][j][5];
f[i][j][6] = 0.5*(rho[i][j] - f[i][j][1] - f[i][j][2] - f[i][j][3] - f[i][j][4] - f[i][j][5] - f[i][j][0] - f[i][j][7]);
f[i][j][8] = f[i][j][6];
}
void output(int m)
{
ostringstream name;
name << “cavity_” << m << “.dat”;
ofstream out(name.str().c_str());
out << “Title=“LBM Lid Driven Flow”\n” << “VARIABLES=“X”,“Y”,“U”,“V”,“P”\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 << i << " "
<< j << " " << u[i][j][0] << " " << u[i][j][1] << " " << rho[i][j] << endl;
}
}
void Error()
{
double temp1, temp2;
temp1 = 0;
temp2 = 0;
for (i = 1; i<NX; i++)
for (j = 1; 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;
}