# Doubt in calculating the drag coefficient

Posted by muzhiyu
 Doubt in calculating the drag coefficient May 23, 2017 02:12PM Registered: 2 years ago Posts: 1
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.
My code are as follw:
```Language: C++#include<iostream>
#include<cmath>
#include<cstdlib>
#include<iomanip>
#include<fstream>
#include<sstream>
#include<string>

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.25*NX;
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.0*eu + 4.5*eu*eu - 1.5*uv);
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;
}```
 Re: Doubt in calculating the drag coefficient August 30, 2017 12:56PM Registered: 2 years ago Posts: 1
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
Sorry, you do not have permission to post/reply in this forum.