Pressure gradients with periodic boundary conditions

Lattice Boltzmann Method LBM: Theory
#1

Dear all,

I’m a physics PhD student, working on a problem that requires a flow field as input, and so far our idea is to get this flow field by simulating it using LBM. To get acquainted with the method, I equipped myself with the C sample code from the website (Big Thank You, Jonas), the Shen/Doolen review, Jonas’ thesis (I really don’t get the regularized part yet :D) and ’ book… and while I understand a lot more than before now, there are still a quite a few places where I’m either not sure if I got the information right or still totally clueless. So, I thought I just write down what my problem is and what I tried to understand it so far, and if anyone has the time to give me some pointers it’ll be most appreciated.

Generally, we need the flow through a device that consist of a huge - basically ‘infinite’ - array of microscopic obstacles. Dimensions will be in the thousands of obstacles in both dimensions, so we plan to basically simulate one obstacle with periodic BCs only. Flow through the device will be driven by a pressure gradient over the entire device, so we’d expect to have a similar pressure gradient over our cell, too. Hence, we’d sort of have to merge the pressure gradient with periodic boundary conditions. So, first question - I read here in the forum Timm writing that “pressure gradient is the same as body force”. Now I’ve no idea what exactly a “body force” is, but Sukop talks on page 54 about implementing “forces”, such as gravity, as drivers for flows. Now I understand the concept for gravity, so my question: 1) Is this what Timm meant with “body force”? And directly, 2) is this indeed the same? Especially, a force like gravity would act on every liquid molecule (or, in LBM, node), while the pressure gradient is only defined over the entire length of the device - I’d expect that at least under some geometries some parts of the system would be less affected under by the pressure flow than under gravity.

Ok, onwards to playing around with the code: Jonas’ code implements all sorts of BCs, but actually uses a Poiseuille velocity condition at the inlet and a pressure BC on the outlet, which - as far as I understand it - has varying pressure in order to eliminate the pressure gradient. ( another question, but not really important: 3) I reckon this “zero gradient” corresponds to having the cylinder ending in some huge “reservoir”. Now, if I have two such reservoirs – how would I use two zero gradients of different value??)

Now, my first idea was to implement just a constant pressure BC at the inlet and an equivalent constant pressure condition at the outlet, to get an idea how such a pressure gradient works anyway. Well, that sort of blew up in my face :), I changed the output to rho and it looked as if the sound waves coming from the pressure were resonating till pressure somewhere exceeded 1.2, which seems to be the point where the Taylor expansion is breaking down. So, 4) what would be the standard approach to implement a pressure gradient?

Anyway, since this is not what I wanted to have anyway, I thought I just take periodic BCs, but with a prefactor:

for (iY=1; iY<=ly; ++iY) {
lat[1][iY].fPop[1] = lat[lx+1][iY].fPop[1]*presGrad;
lat[1][iY].fPop[5] = lat[lx+1][iY].fPop[5]*presGrad;
lat[1][iY].fPop[8] = lat[lx+1][iY].fPop[8]*presGrad;

    lat[lx][iY].fPop[3] = lat[0][iY].fPop[3]/presGrad;
    lat[lx][iY].fPop[6] = lat[0][iY].fPop[6]/presGrad;
    lat[lx][iY].fPop[7] = lat[0][iY].fPop[7]/presGrad;
}

I’ve set presGrad to 1.01 so far. Since this will generate a net flow to the right, I’m also re-normalising the flows such that I always have the same amount of ‘fluid’ in the system. I was under the impression that this was necessary in order to keep rho close to 1 for numeric stability, however after reading the forum I’m not too sure anymore?
This generally seems to work o.k. for a while, and I get the Poisseuille flow and the fluctuations behind the obstacle, but after a while it blows up – pressure at the obstacle increases further and further and then explodes quickly after it reaches 1.1. So my last question for now would be, 5) has anyone done this so far, and if so, are there any example codes/guidelines/suggestions available?

OK, thanks everyone for reading, I appreciate any input. Can be very short, i.e. “read this paper” or “can be done with the C++ lib by using xyz” will suffice.

Thanks a lot,

Oliver

1 Like
#2

Dear Oliver,

actually, my statement was a bit unclear. If you have a look at the Navier-Stokes equation, you will see that a pressure gradient as well as the external force appear. From a lattice Boltzmann point of view, you can use a pressure gradient to drive the flow (at inlet and outlet), or an external force density (I call it body foce) or both. Right now I am also simulating a system which is periodic in x-direction. And I find it much easier to use a body force to drive the flow. How do I do it? I add a source term to the collision step, adding momentum to the fluid at each time step and each lattice point. There is a long list of literature. The most important article is by Guo, Zheng, and Shi: “Discrete lattice effects on the forcing term in the lattice Boltzmann method” (2002). You should study it, it may be useful for you.
However, in a complex geometry, the results obtained from both methods (body force or pressure gradient) will be different. The reason is that the pressure gradient develops in the flow, and it will have a shape which also depends on the geometry. This cannot be predicted ahead of the simulation in general cases, where the body force has to be defined at any point at any time step explicitly. So, if you want to simulate porous media, a pressure gradient driven flow will be more accurate (I think).

I do not know how to implement periodic boundary conditions with a pressure gradient, but I see the problems arising from it.

Concerning the mass increase: Using only a body force, the total mass will always be exactly constant (up to machine accuracy), as long as the code is free of errors. The presence of inlet/outlet boundary conditions violates the mass conservation due to the equation of state of the LB fluid (it is slightly compressible). If your mass increase is significant (say, more than 10% after the entire simulation time), you have the option to renormalize it. However, I do not exactly know, which renormalization method is the best and influences the correct macroscopic solution least. I am not aware of any articles investigating this effect very thoroughly.

Best,
Timm

#3

Hi Timm,

thanks for the quick reply. Read and thought about the paper and decided to - so far - not to use a body force. I read in several papers that “a body force is equivalent to a pressure difference for channels with uniform cross-sectional area” which is precisely what I don’t have :slight_smile: I found a paper (PRE 73, 047702, Zhang and Kwok) that proposes a way to employ a pressure gradient with PBCs and it seems to work for them. However, their pressure difference is much smaller than what I used (factor 10), so maybe my problems stemmed just from having too high a gradient. I’ll try to get their - slightly different - approach to work with the C code and will report back here, maybe someone else will have the same problem.

There were however a number of questions that popped up when further getting into the code of the c-example and the c++ library. I’ll just ask them here, with a little bit of luck someone might have insight:

A) In the C-example, we have a walls that are zeroVelocityBoundary-s. What is the difference in physics between using zero velocity inflow and solid bounce-back nodes for walls? Is ther any?

B) Simulation size - out of your experience, what is the maximum size of a system that I can compute on a reasonably fast quad core machine? Just order of magnitude - 10000, 100000, 10^6 Nodes in total? interested in finding the steady state.

Alright, thanks everyone for their work again!
btw, is there a way to change the username? I could sign up with unicode chars, but looks like the forum doesn’t like it…

Oliver

#4

Dear Oliver,

regarding your question B), there is a rule of thumb. A fast LB simulation has appr. 1 Msu/s or more, i.e. 10[sup]6[/sup] lattice sites can be computed in one second using one CPU (for D3Q15; D3Q19 is much slower). If your simulation requires 10[sup]5[/sup] time steps, a total of 10[sup]6[/sup] lattice nodes using four CPUs will probably require 7 hours. If you use some heavy computations, this will slow down your speed, possibly very significantly. I think, 7 hours will be the absolute minimum.
For example: I have run a simulation with 10[sup]6[/sup] time steps and a lattice size of 10[sup]4[/sup] nodes last Friday. It was completed after 15 hours on 1 CPU, thus I had 0.19 Msu/s. But my D3Q19-LB part of the simulation was only 30% in computing time, and I had additional computations in the LB part, which are often not necessary.

Question A)
Both boundary conditions result in (nearly) zero velocity on the walls, where the zero-velocity boundary condition surely is more accurate. However, it is much more complicated to implement it for non-trivial geometries. If you want to run a high-accuracy simulation of a simple geometry, use velocity-boundary conditions for the walls. If you want to have good results in a complex geometry, try bounce-back.

I hope that helps,
Timm

#5

A question regarding units… The sample code basically takes a maximal inflow velocity as input and then calculates the viscosity and omega from the Reynolds number. On the contrary, we want to simulate a physical system, i.e. one where we know the viscosity already - water has ~10^{-6} m^2/s. Our system is super tiny - 10 micrometers, and the flow would be probably something like 10 \mu m /s, so l_0=10^{-5} and t_0=1. We tried to relate this to the quantities used in the LBM simulation.

For the viscosity, it means that is is expressed as 10^4 l_0^2/s, which is of order of magnitude of 10^6 higher than the default in the c example code. For the pressure gradient, we expect it to be like 10 Pa over this device (0.1 bar over an array of 1000 devices, approx 1cm long) which appears to be quite reasonable to us. However, in terms of dimensionless variables, this would mean that 10Pa= \rho l_0^2/t_0^2 P_d = 10^{-7} Pa p_d ==> p_d=10^8

And this is where we think we messed something up. Is LBM absolutely unsuitable for microscopic flows? Or do we have to set nu=10^4 and the pressure boundaries to 0.95\cdot 10^9 and 1.05\cdot 10^9 respectively and just run the calculations? We followed basically all the steps in the raspberry jam example in the units-howto and I’m a bit clueless…?

Edit - I looked at Eq. 13 of course, but given that the viscosity was of the same order of magnitude for our discretisation (500 nodes in x-direction, figured we’d also use a \delta t = 10^-5), I ignored it for now, since I didn’t know how or if the discretisations would enter the pressure.

Thanks a LOT,

Oliver

#6

Hello,

this is how I usually do this (p for physical value, L for lattice value):
First collect all relevant data, so nu[sub]p[/sub] = 10[sup]-6[/sup] m[sup]2[/sup]/s, L[sub]p[/sub] = 10[sup]-5[/sup] m, u[sub]p[/sub] = 10[sup]-5[/sup] m/s
From this we know the Reynolds number: Re = L[sub]p[/sub] u[sub]p[/sub] / nu[sub]p[/sub] = 10[sup]-4[/sup], which is rather small.
Then I choose a relaxation time, say tau = 0.9 and thus nu[sub]L[/sub] = (0.9 - 0.5) / 3 = 0.1333. Furthermore I choose a resolution. Assume I want to have 20 lattice nodes across the channel, then dx[sub]p[/sub] = L[sub]p[/sub] / 20 = 5 * 10[sup]-7[/sup] m. The time step comes from nu[sub]p[/sub] = nu[sub]L[/sub] dx[sub]p[/sub][sup]2[/sup] / dt[sub]p[/sub] and thus dt[sub]p[/sub] = 3.33 * 10[sup]-8[/sup] s. This is really small. You may increase it by increasing dx[sub]p[/sub] which would make your spatial resolution worse, or by increasing nu[sub]L[/sub] which would drive the relaxation time beyond 1.
The dimensionless pressure gradient would then be dp[sub]L[/sub] = 4.44 * 10[sup]-5[/sup] according to dp[sub]p[/sub] = dp[sub]L[/sub] rho[sub]p[/sub] dx[sub]p[/sub][sup]2[/sup] / dt[sub]p[/sub][sup]2[/sup] with rho[sub]p[/sub] = 1000 kg/m[sup]3[/sup]. This is quite reasonable. So, there must be an error in your computation.

Timm

#7

Hi Tim,

first of all, thanks for your elaboration. I’ve understood my fault with the unit conversion, this part works very well now. However, when I tried to summarize everything I did I noted something peculiar both in your explanation as well as in the example file. You say you want to keep the relaxation time below one. I reckon you refer to \tau as used in

f_i(t+1)=(1-1/\tau)f_i +1/tauf^{eq}_i

I’d very much expect \tau to be larger than unity in lattice time steps, such that complete relaxation takes more than one time step. If \tau is less than unity, this means the decay goes sort of beyond equilibrium into non-equilibrium (of a different kind), which wouldn’t be what I learned in my Stat Mech :wink: So, first I thought you had somehow confused 0.9 and 1/0.9, but I cross-checked with the C example code, and there indeed \tau=0.5, so I assume I must have misunderstood something?

Best Regards,

Oliver

#8

Hello,

you are right that a larger relaxation time is equivalent to a slower rate of equilibriation. This is true from a physical point of view. However, the important point is how the numerical scheme works. You can check in the literature that usually values smaller than 1 for tau are taken. I principle, you can also set tau = 10, but you will get numerical difficulties and increased inaccuracies, as we have pointed out in our article. I always try to take tau merely as a numerical parameter, although its physical significance cannot be fully ignored.

Timm

#9

Hi Timm,

read through your article. Helped in quite a few points, so thanks for that - but not in this one point where I have the issue :wink:

Let me put it this way: If you have a linear relaxation from one state X into another state Y, at any point of the relaxation you should find some state aX+bY with a+b=1 and, and here is my problem, 0<=a,b<=1. If you rewrite your eqn. 1 you’d get f_i(x+c_i, t+1)=(1-1/tau) f_i(x,t) +1/tau f^eq_i(x,t). Hence, the condition 0<a,b<1 is only fulfilled for \tau>1. Also, a relaxation should proceed exponentially, but with \tau<1 you get an oscillating relaxation. Worse, with \tau<0.5 you blow your whole setup up because the values are probably becoming negative somewhere (The c-code example circumvents this problem by defining \tau in a way such that it \lim_{R\to\infty} \tau=0.5+. Also, if you set \tau=1 then there would be no leftover non-equilibrium at all…

I’m trying to write this down in details and will send you the PDF with a figure, I think phorum is not perfect to discuss something with formulas and figures. Of course simply ignore it if you’re not free, but I think if I got this problem right then it is something that we might add to the unit conversion instructions, since I’m probably not the only one with this problem.