Doubt in calculating the drag coefficient

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.5
NY;
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;

}

Hi,muzhiyu,

Would you please tell me the meaning from line130 to line164 in your code?
I am a novice programmer and trying to use the LBM on the classic circular cylinder unsteady problem,
and I’ll be grateful for your replied.

Regards
Mia