Foces calculation over stationary solid body (cylinder)

Dear all,
I am trying to simulate the flow over cylinder using D2Q9 algorithm, and I have to calculate the forces.
I have used the Ladd’s method for calculating forces, but the problem is after some time force values starts decreasing.
Can anyone please find the mistake in my algorithm

Steps for simulation :

equilibrium();
collision();
bounceback(); // on upper and lower wall
BC_calc(); // FH boundary condition on surface of cylinder
streaming();
forces = force_calc(); // force calculation using algorithm attached here in question
velocityinit();
err = errorcheck();
terminaloutput(i, err, forces);
vtkoutput(i);

These steps are iterated…

the algorithm I am using for calculation of forces

[code=“cpp”]
double force_calc()
{
static double force[2];
double FX, FY, Xforce, Yforce;
FX = 0.0;
FY = 0.0;
for (int i = 0; i < nx+2; i++)
{
for (int j = 0; j < ny+2; j++)
{
isite = (ny
i) + j;
if (node_type[isite] == 1)
{
Xforce = 0.0;
Yforce = 0.0;
for (int k = 1; k < q; k++)
{
isite = (ny*(i+ex[k])) + (j+ey[k]);
if (node_type[isite] == 0)
{
Xforce += (-2.0ex[k]f[k][i+ex[k]][j+ey[k]]);
Yforce += (-2.0
ey[k]f[k][i+ey[k]][j+ey[k]]);
}
}
FX = FX + Xforce;
FY = FY + Yforce;
}
}
}
force[0] = 2.0
FX/(rho_init
2radu0u0);
force[1] = 2.0
FY/(rho_init2radu0u0);
return force;
}