Drag and lift coefficients for turek benchmark

[yevRiz]

Hi Jonathan,

I tried to repeat the benchmark calculations with your code and got different values for the drag and lift coefficients. E.g. for Re = 100 and L_lbm=ny=100 I tried 3 different velocities (original value 0.05; 0.02 and 0.007), i.e. different relaxation rates. The values for the average drag coefficient differ only by around 1% from the 3.5409 that you provide in the README file. However, the average lift coefficient is always around 0.11, falling below the specified value by a factor of 10. Did you perform these calculations with other parameters than in the uploaded code ('magic parameter' etc.)?

I would also appreciate a comment on the calculation of the drag and lift coefficients. According to your code, they are calculated based on first-order moments of the sum of collision and total distribution functions. According to several sources, it should be second-order moments of the non-equilibrium part of the distribution function (https://journals.aps.org/pre/abstract/10.1103/PhysRevE.79.046704).

Thank you in advance for your reply! Best, Yevgeniy

[jviquerat]

Hi @yevRiz,

First, thanks for showing interest in this work. I should underline that this is a side project, and that I'm definitely not an expert in LBM methods, the goal was just to try a few things and have a bit of fun. Just ran the Turek benchmark again, and something definitely seems wrong. Yet, I'm certain I had these results at some point, so I must have broken something in the last commits. I will try to look into it, although I have very limited time to devote to it.

Regarding the computation of drag/lift, I mostly implemented stuff I found in different books/websites, but the details are not exactly clear in my head. If you think the computation needs to be fixed, I would gladly accept a PR.

I'll keep you posted in the following days/weeks.

Best

[jviquerat]

Hey @yevRiz,

As written in the main readme:

Note that for the 2D-2 case, the values correspond to the max drag and lift.

It took me a few tens of minutes to remember that this is how things are defined in the Turek benchmark for Re=100 ;) An option to get the maximum value of a column in a csv-like file:

awk -v col="<column_number>" 'BEGIN {max=-100000} {if ($col>max) max=$col} END {print max}' <file_name>

Yet, if you're interested in implementing what you're suggesting to improve drag/lift computation, I'm still open to a PR ;)

[yevRiz]

Hi Jonathan,

Thank you very much for your response and for taking the time! For me, as an LBM beginner, your code is very helpful for teaching purposes and I am very grateful for sharing this project!

Indeed, I treated the coefficients from the 2D-1 and 2D-2 cases in the same way by considering the moving average. Since the drag values are negative, I took the minimum value now, which matches the one given in the readme pretty well. Also for lift, the minimum value agrees well with the readme. However, the amplitude is much larger and the maximum (positive) has a greater absolute value. Perhaps it would make sense to take this value as the max absolute value for the benchmark comparison (1.324 instead of 1.061).

For the 2D-1 (Re=20) benchmark case, I see a similar problem with the lift coefficient: the minimum of -0.0202 I get is close to the value you provide, but the maximum of 0.0399 is twice as large in absolute value. Besides, this stable system converges to the last value.

Do you obtain similar results?

Regarding the computation of drag/lift, so far I have only found some analytical derivations for the BGK-scheme, not for the TRT. So I was just wondering if you derived it yourself and why there is such a big difference between these two schemes: using a function + (sorry, just "+" doesn't want to work) instead of under the integral. And also I am curious about the influence of the bounce-back approach as BC: the lattice points are placed close but not at the surface in contrary to the analytical calculations. As already mentioned, I am a beginner with LBM and the force calculation is not of interest to me for the purpose of my research. But, if I find a way to improve the accuracy in the benchmark cases, I will gladly share my suggestions!

Best, Yevgeniy

P.S.: here is the code I used for the visualization

import numpy as np
from matplotlib import pyplot as plt
plt.rcParams.update({'font.size': 7})

data = np.loadtxt('./lbm_jviquerat/results/turek_Re20u05/drag_lift')
'''
### **********************************************
### 2D-2
#Drag coeffcicent

c_d_max = abs(np.amin(data[:,1]))
plt.figure(dpi=200)
plt.text(0, 2, '$c_{Dmax} = $%s'%('{:.4f}'.format(c_d_max)))
plt.plot(data[:,1], "-b")
plt.axhline(y=-3.24, color='r', linestyle='--')
plt.axhline(y=-3.22, color='r', linestyle='--')
plt.xlabel('timestep')
plt.ylabel('Drag coefficient')
plt.title('2D-2: Re = 100; u = 0.05. Drag coefficient')
plt.show()

#Lift coeffiecient
c_l_min = np.amin(data[:,2])
c_l_max = np.amax(data[:,2])
plt.figure(dpi=200)
plt.text(0, 0.7, '\n'.join(('$c_{Lmax} = %.4f$' % c_l_max,
'$c_{Lmax} = %.4f$' % c_l_min)))
plt.plot(data[:,2], "-b")
plt.axhline(y=-0.99, color='r', linestyle='--')
plt.axhline(y=-1.01, color='r', linestyle='--')
plt.axhline(y=0.99, color='r', linestyle='--')
plt.axhline(y=1.01, color='r', linestyle='--')
plt.xlabel('timestep')
plt.ylabel('Lift coefficient')
plt.title('2D-2: Re = 20; u = 0.05. Lift coefficient')
plt.show()
'''
### **********************************************
### 2D-1
#Drag coeffcicent
c_d_max = abs(np.amin(data[:,1]))
plt.figure(dpi=200)
plt.text(30000, -1, '$c_{Dmax} = $%s'%('{:.4f}'.format(c_d_max)))
plt.plot(data[:,1], "-b")
plt.axhline(y=-5.57, color='r', linestyle='--')
plt.axhline(y=-5.59, color='r', linestyle='--')
plt.xlabel('timestep')
plt.ylabel('Drag coefficient')
plt.title('2D-1: Re = 20; u = 0.05. Drag coefficient')
plt.show()

#Lift coeffiecient
c_l_min = np.amin(data[:,2])
c_l_max = np.amax(data[:,2])
plt.figure(dpi=200)
plt.text(0, 0.02, '\n'.join(('$c_{Lmax} = %.4f$' % c_l_max,
'$c_{Lmax} = %.4f$' % c_l_min)))
plt.plot(data[:,2], "-b")
plt.axhline(y=0.0104, color='r', linestyle='--')
plt.axhline(y=0.011, color='r', linestyle='--')
plt.axhline(y=-0.0104, color='r', linestyle='--')
plt.axhline(y=-0.011, color='r', linestyle='--')
plt.xlabel('timestep')
plt.ylabel('Lift coefficient')
plt.title('2D-1: Re = 20; u = 0.05. Lift coefficient')
plt.show()```

[yevRiz]

Oh, I have a small error in the code and in the diagrams for the lift: the negative values are obviously for

[jviquerat]

Hi @yevRiz,

This is linked to the discretization of the cylinder (obviously the staircasing effect is very strong, and has a large influence on the development of the instability). I suggest you try to reproduce the results with the same parameters as I suggested in the readme. For u=0.03 and ny=100, here is the plot I get for re=100:

While for ny=200, I obtain a symmetric plot with an almost-zero average:

I guess that's the source of the discrepancies you are observing.

Regarding the drag-lift computation, I have not much to say, I won't have time to go back in details into these derivations, and I can't even remember how I got to the way I compute drag and lift currently. If you want to experiment and share, be my guest ;)

Best

[edmald]

Hi Guys, I am just starting out with Lattice Boltzmann Methods (LBM) and I've encountered an issue that I hope you can help me with. In my recent simulations, I am observing negative values for both drag and lift coefficients. Based on my understanding, I was expecting these values to typically be positive. Why do we obtain negative values with this Library? It is some kind of convention? Thanks in advance.

[jviquerat]

Hi Guys, I am just starting out with Lattice Boltzmann Methods (LBM) and I've encountered an issue that I hope you can help me with. In my recent simulations, I am observing negative values for both drag and lift coefficients. Based on my understanding, I was expecting these values to typically be positive. Why do we obtain negative values with this Library? It is some kind of convention? Thanks in advance.

Hi,

Drag and lift forces signs originate from their definitions: they are the surface integrals of the projections of the stress tensor over the axes. Depending on the definition of the axes (which is purely conventional), they can end up being positive or negative.