About the drag coefficient calculation based on Moment-Exchange Scheme

Lattice Boltzmann Method LBM: Theory
#1

Hi guys,

I got a problem during the calculation of drag coefficient of a circular cylinder as in the example problem of Cylinder2D.

For Re=100, the computation yields drag coefficient Cd = 12.5 (Cd = Fd/ u \rho^2 r). I cannot find out what wrong with my codes. For each boundary lattice nodes, I first sum up all distribution functions which are non zero as follows,

   T Fd = T();
	typedef descriptors::D2Q9DescriptorBase<T> L;
	if (grid[iX][iY][0]==0){
	for (int iPop=1; iPop<=8; ++iPop) {
		if ( grid[iX][iY][iPop] != 0 ){
			int nextX = iX + L::c[iPop][0];
            int nextY = iY + L::c[iPop][1];
            if (iPop>4)
            	Fd -= ( ( grid[iX][iY][iPop] + grid[nextX][nextY][iPop-4] ) * L::c[iPop][0] );
            else
            	Fd -= ( ( grid[iX][iY][iPop] + grid[nextX][nextY][iPop+4] ) * L::c[iPop][0] );
    }

Then, I sum up all the Fds on boundary nodes, and divided by (\rho u^2 r). As in the example, rho = 1, u = 0.02, r = 7.

The final Cd is approximately 12 which is apparently wrong. can anybody here give me a help plealse? Thank you.

#2

Hi Dear Friend,
I used the momentum exchange method like this, may be it help you.

// calculate force on body
void get_force(int wall[NX][NY], int wb[NX][NY], double *f_x, double *f_y)
{
int i,j; *f_x = 0.0; *f_y = 0.0;
for(i=1;i<NX-1;i++){ int ip = i+1; int im = i-1;
for(j=1;j<NY-1;j++){ int jp = j+1; int jm = j-1;
if(wb[i][j]==1){

	        *f_x = *f_x - (f_1[im][j][3]+f_1[i][j][1])* (1-wall[im][j])       + (f_1[ip][j][1]+f_1[i][j][3])* (1-wall[ip][j])
   		            - (f_1[im][jm][7]+f_1[i][j][5])*(1-wall[im][jm]) 	    + (f_1[ip][jp][5]+f_1[i][j][7])*(1-wall[ip][jp])
        	        + (f_1[ip][jm][8]+f_1[i][j][6])*(1-wall[ip][jm])  	            - (f_1[im][jp][6]+f_1[i][j][8])*(1-wall[im][jp]);

	        *f_y = *f_y - (f_1[i][jm][4]+f_1[i][j][2])*(1-wall[i][jm])        	+ (f_1[i][jp][2]+f_1[i][j][4])*(1-wall[i][jp])
    		     	- (f_1[im][jm][7]+f_1[i][j][5])*(1-wall[im][jm]) 		 	+ (f_1[ip][jp][5]+f_1[i][j][7])*(1-wall[ip][jp])
       		        - (f_1[ip][jm][8]+f_1[i][j][6])*(1-wall[ip][jm]) 		        + (f_1[im][jp][6]+f_1[i][j][8])*(1-wall[im][jp]);

		}
	}
}

}

where wall means the top and bottom boundaries. wb means surface of the object.

Then you can calculate the drag and lift

       cdx = -fx/rho[2][(NY-1)/2]/ui/ui/radius;
   cdy = -fy/rho[2][(NY-1)/2]/ui/ui/radius;
#3

thank you very much! Khan
by the way, is rho[2][(NY-1)/2] a constant and equal to 1.0?
what are the values of fx and cdx for the cylinder in the channel ( Re =100 and uMax=0.02)?
thanks again

#4

Hi Dear Friend,
rho[2][(NY-1)/2 my friend not constant, it’s varies for few thousand iterations and then it’s becomes constant. But the variation is very small almost negligible. From the above function we access the values fx and fy in the main code. If you don’t mind send me your code i will check it for you.

void main()
{

int t;
int i,j;
int wall[NX][NY], wb[NX][NY];
double fx, fy;
double cdx, cdy;
get_force(wall, wb, &fx, &fy); // force on fluid by body
cdx = -fx/rho[2][(NY-1)/2]/ui/ui/radius;
cdy = -fy/rho[2][(NY-1)/2]/ui/ui/radius;
}

#5

Dear khan
Is it possible for you more detaile explanation of this nice and short code such as wath is the diffrence between wall and wb?

#6

Hi Dear srtmlb,
My friend in my code i used wall for wall boundaries of the computational domian. And wb i used for the solid object surface like circular cylinder.
Thanks

#7

Dear khan
I transfered this code to fortran code and used it.i think its lift calculation is acceptable but for drag its work wrong.
is this problem happen for you too?
thanks.

#8

Hi Dear srtmlb,
Dear srtmlb check this one also and also read the private message.
The implementation is good the other way to find forces.
void drag_lift()
{
int i,j,k,id,jd;
double rrho;
FD=FL=0.0;
for(j=1;j<NY;j++) for(i=1;i<NX;i++) for(k=1;k<9;k++)
{
jd=j-e[k][1]; id=i-e[k][0];
if(wall[j][i]=1)
{
FD-=(e[k][0](f[j][i][k]+ff[j][i][r[k]]));
FL-=(e[k][1](f[j][i][k]+ff[j][i][r[k]]));
}
}
rrho=2./(Lu0u0);
FD*=rrho; FL*=rrho;
}

where e[k][0] for lattice velocities in axes and e[k][1] for lattice velocities along diagonal directions.
For example e[1][0]=1; e[0][1]=0;
r[k] for opposite direction, for example r[0]=0; r[1]=3; r[3]=1;
rrho i used because i try this for incompressible model proposed by Guo et.al.

You can try this one no problem in this one .

#9

Dear Khan,

Could you please tell me what is the incompressible model you are talking about. You said it is the one proposed by Guo et al… Is it that from the paper: Guo Z et al. J. Comput. Phys. 165: 288 (2000) ?

If so, did you try to implement the He and Luo incompressible model (He and Luo, J. Stat. Phys. 88: 927 (1997))? Did you find any big difference between both?

Thanks in advance for your answer.

Regards

Goncalo

#10

Hi Dear Friend GoncaloSilva,
Yeah you are right i used that proposed Guo et al model. I not do a comparison between these two models. I do the comparison of this proposed Guo model and MRT incompressible model proposed by Rui Du et al. I tested these models for circular cylinders in tandem arrangement. I find some results and some differences using these models. Actually i am an average student not so intelligent but work hard. I facing problems to write a function for streamlines which shows some more features about the flow behaviour. I just have vortex formation. Almost one friend through this side ask me to send a code and he will help me to make it more better and find some mistakes. Once i send the code then he not contact any more. If possible can you help me if you have a function written in c++ for streamlines and send it to me or post here. It will be a great help. My tandem arrangemet code is very good and not find any troubel.
Thanks Friend

#11

Hi Khan,

Thank you for your quick answer. My concern regarding the incompressible models is due to the fact that I am myself studying these models… I found out that when modeling steady incompressible flows they seem more appropriate than standard ones (well in fact this is not new I only achieve to the same conclusion)… However I still do not have enough tests demonstrating that for my study in particular.

Sorry about the C++ streamline function but I do not have any (I use matlab for simple 2d implementations and fortran for heavy stuff and in neither cases I have studied streamlines because they do not bring much information to me, my geometries are planar so streamlines will just follow the walls pattern).
However, I believe that such function is easy to find. Did you look in openlb code? I think they have some coding for that. Nevertheless you can google for streamline functions in c++ libraries… I bet someone has already put such work available in internet.
Matlab comes with a library that computes streamlines…other alternative is to base your code on that…

Good luck for your work

Regards

Goncalo

#12

My new friend khan please check your PM.

#13

Hi everyone,

My code for calculating the drag and lift forces is as follows. It is written in fortran



    INTEGER:: ex(0:8) = (/0,1,0,-1,0,1,-1,-1,1/)  
    INTEGER:: ey(0:8) = (/0,0,1,0,-1,1,1,-1,-1/)  
    INTEGER:: opposite(0:8) = (/0,3,4,1,2,7,8,5,6/)
        F_X=0
        F_Y=0
        DO x=(obstX-(obstR+2)), (obstX+(obstR+2)) ! reduced domain 
            DO y=(obstY-(obstR+2)), (obstY+(obstR+2)) ! reduced domain
                fx=0
                fy=0
                IF (image(x,y)==cylbdry) THEN  ! Boundary points on cylinder
                    DO j=1,8
                        xnew = x + ex(opposite(j))
                        ynew = y + ey(opposite(j))
                        IF (image(xnew,ynew)/=cylwall) THEN ! interior cylinder nodes.
                            fx = fx + ex(opposite(j))*( fplus(x,y,j)  + fplus(xnew, ynew, opposite(j)) )
                            fy = fy + ey(opposite(j))*( fplus(x,y,j)  + fplus(xnew, ynew, opposite(j)) ) 
                        END IF
                    END DO
                END IF
                F_X=F_X+fx
                F_Y=F_Y+fy
            END DO
        END DO

I don’t know what is the mistake, my values for drag are differing. Please suggest.

Thanks,
Narender

#14

I’m sure this is going to be obvious so I hope I don’t sound too silly and/or patronising. The value of the drag coefficient depends strongly on the aspect ratio of the channel (that is, the ratio of the width to the length). For example, at Re=100 the drag should be about 1.92 if the aspect ratio is 0.2, and about 10.8 if the ratio is 0.7. If the flow if periodic in the vertical direction then the drag is about 1.2. If you’ve implemented the boundary conditions correctly (eg, no slip on all solid walls and appropriate inflow and outflow conditions), and if the flow is fully developed by the time it reaches the cylinder, then the momentum exchange method should give reasonable results.

Regarding incompressibility, for steady flows the model of He and Luo has a compressibility error term an order in Mach number better than “standard” models (that is, it approximates the incompressible Navier-Stokes equations with a Ma^3 error). But for unsteady flows the compressibility error is of the same order for both models, I think.

#15

Oh yes, some packages come with the ability to plot streamlines from the velocity field (I know Tecplot can do this, for example). However, it is not always clear how they are doing their integration and they are not always appropriate for extracting data (such as the vales of the function at specific places). As well as any streamline finding libraries, Matlab also has a fast solver for Poisson’s equation, which can be used to compute the streamfunction from the vorticity (remember that the streamfunction psi satisfies d^2\dx^2 phi = -omega, where omega is the vorticity - found by taking the curl of your velocity). Or you can use Simpson’s rule or the Trapezium rule to find the streamfunction from the velocity.

#16

Hi Pleb01,

Thanks for the reply. I am using a parabolic velocity profile (Fully developed) as inlet condition. I am getting the required drag coefficient but it is oscillating too much without any regular pattern between the maximum and the minimum. I think I have some problem in calculating the drag coefficients.

Thanks,
Narender

#17

Have a look at the flow field around the cylinder to make sure nothing weird is happening. If you can get high enough Re then look at the period of vortex shedding too. The drag can oscillate for Re>100 and a small aspect ratio but there should be a max and a min. It’s hard to understand how the momentum exchange method can give big oscilliations unless the flow field is oscillating (either that or you’re taking gradients from strange places or your cylinder is nor correctly defined)

#18

Hi Pleb01,

Thanks for your comments. I have checked my cylinder boundary marking algorithm and found a bug I fixed it and running the test case again. Hopefully it will be fixed.

Thanks,
Narender

#19

Hi Everyone,

I am getting the drag to oscillate but it has two maximums and two minimums. Did anyone face the similar issue? Please suggest. I am working on asymmetrically placed cylinder in a channel.

Thanks a lot.
Narender

#20

Hi
Dear All
I Read the mentioned code about calculation of Drag Coefficient I use Half-way bounceback on wll and I want to calculate Drag coeffiecient in the case of flow over squre please help me in this area. I don’t understand above code. is this possible to calculate this coefficient for my problem fluid flow due to body force over square and using bounceback on wall with periodic boundary ? and if yes, How?

Tnx In Advance
Sajjad