Welcome! Log In Create A New Profile

Advanced

Doubt in calculating the drag coefficient

Posted by muzhiyumuzhiyu  
Doubt in calculating the drag coefficient
May 23, 2017 02:12PM
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:
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
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.