Ideal gas equation of state

Dear all,

LBM uses the equation of state of an ideal gas, P = c_s^2 rho, to obtain pressure P from density rho. Still, the LBM is capable of simulating liquids and fluids other than an ideal gas. I do not really understand how this is possible. The thread at http://www.palabos.org/forum/read.php?3,899 discusses this topic also somewhat, but does not really give a satisfying answer to my question. At http://www.scholarpedia.org/article/Lattice_Boltzmann_Method#Idiosyncrasies.2C_lucky_strikes_and_some_misbeliefs they also try to explain this. They state that for incompressible flows this equation of state does not play a role, which is true. But for compressible flows it does play a role and they do not really explain why the LBM can then still be used for fluids other than an ideal gas. Or can be we only recover the Navier-Stokes equations from the LBM (via CE) by using the ideal gas law? Can somebody explain this to me?

Best regards,
Bart

Hi,

The most commonly encountered isothermal LBM can actually not even properly simulate an ideal gas with an equation of state p = \rho R T. In an ideal gas, an adiabatic increase of the pressure p would increase both the density \rho and the temperature T (how much each of them increases depends on the heat capacity ratio \gamma of the fluid). However, since the temperature is held constant in isothermal LB simulations, only \rho increases with p.

The speed of sound in a fluid links the increase in pressure with only the increase in density; speed of sound in any fluid is given as

c_s^2 = (\partial p / \partial \rho)

evaluated at the “rest state” of the fluid (e.g. atmospheric pressure and room temperature, or whatever else). For the isothermal LB model, where R T is constant, this simply gives c_s^2 = R*T. For an ideal gas, we get c_s^2 = \gamma R T[sub]0[/sub], where T[sub]0[/sub] is the rest state temperature.

In any case, this equation for the speed of sound can linearized, giving an equation of state for small-ish changes around the rest state,

\Delta p = c_s^2 \Delta \rho,

where \Delta p = p - p[sub]0[/sub] and \Delta \rho = \rho - \rho[sub]0[/sub], the zero subscripts indicating the rest state. This equation is actually valid for any equation of state, though only in a region around the rest state. You just have to know the speed of sound of the fluid. The isothermal LB equation of state can also easily be put into the form above.

The point I’m trying to get to is that while the isothermal LB equation of state is not correct for any physically realizable gas (no isothermal ideal gases exist), it can be used to simulate small changes around an rest state for any type of fluid. The ratio between the rest pressure p[sub]0[/sub] and the rest density \rho[sub]0[/sub] will be wrong in the simulation, but this should not be a problem as the absolute value of the pressure is not relevant in the Navier-Stokes equation; only the pressure gradient \nabla p is relevant there.

I hope that helps a bit!
Erlend

Hi Erlend,

Thanks for your quick reply. The statement that P = c_s^2 rho corresponds to an ideal gas confused me, since it can be used for any state.

So, let me verify your explanation for myself. The speed of sound c_s is expressed as
c_s^2 = \partial P / \partial \rho
Next we assume that the partial derivative of P and \rho can be approximated as
\partial P = P - P_0 = \Delta P,
\partial \rho = \rho - \rho_0 = \Delta \rho,
where P_0 and \rho_0 correspond to the pressure and density of the fluid in rest state. Substituting the above approximations in the first equation we obtain
c_s^2 = \Delta P / \Delta \rho = P / \rho and hence
P = c_s^2 \rho.
Perfect, I think I get it now. So our equation of state, P = c_s^2 \rho, is only valid if P and \rho do not deviate too much from rest pressure P_0 and rest density \rho_0.

This statement confuses me again. I know the Navier-Stokes equation uses \nabla p, but can we now use P = c_s^2 \rho to obtain the pressure from the density? Or can we only use it to determine \nabla P?

Best regards,
Bart

Hi,

You’ve nearly got it. :slight_smile: I think it’s easiest if I answer your last question first:

In the lattice Boltzmann method, the equation of state is exactly p = c_s^2 \rho. This means that the pressure gradient term in Navier-Stokes is:

\nabla p = \nabla (p - p[sub]0[/sub]) = \nabla (\Delta p) = c_s^2 \nabla (\Delta \rho)

The first equality holds because p[sub]0[/sub] is constant so that \nabla p[sub]0[/sub] = 0. The last holds for a smallish region around the rest state which the partial derivative was linearized around.

This means that the pressure in Navier-Stokes can instead be expressed by the density change for any eq. of state; at least if you know the speed of sound in the fluid and if heat conduction and such is negligible.

Now, the final step in your derivation above is not entirely correct. You had that

c_s^2 = \Delta p / \Delta \rho.

This is correct for any fluid, for a small deviation of p from p[sub]0[/sub] and \rho from [sub]\rho_0[/sub]. But the simulated LB fluid (and only that), has an equation of state p = c_s^2 \rho, with constant c_s for any p and \rho. We can go the other way to find this equation from the LB eq. of state,

p[sub]0[/sub] = c_s^2 \rho[sub]0[/sub]
=> p - p[sub]0[/sub] = c_s^2 ( \rho - \rho[sub]0[/sub] ) or \Delta p = c_s^2 \Delta \rho
=> c_s^2 = \Delta p / \Delta \rho.

However, this is only exact for the LB equation of state. For any other eq. of state,

\Delta p / \Delta \rho != p / \rho.

Sorry if I’m maybe being a bit confusing. I try not to, but it’s late in my workday and I’m a little tired. :slight_smile:
Erlend M

No more posting for you in the late afternoon :wink:

I’m following you on this part.

I’m not following you here. You say that the LB eq of state is p = c_s^2 \rho. Okay. Then you start your derivation with p0 = c_s^2 \rho0, which is simply one example of p = c_s^2 \rho. Okay, fine. Then the ratio of p and p0 with \rho and \rho0 is equal to the ratio of p0 and \rho0 and equal to the ratio of p and \rho, i.e. c_s^2. But then your last statement says otherwise, by stating that \Delta p / \Delta \rho != p / \rho.

Btw, I tried finding an article explaining this topic, but couldn’t find anything yet. If you know one, I hold myself recommended.

Best regards,
Bart

Hi again!

The point is that you can go from the LB equation of state p = c_s^2 \rho to \Delta p / \Delta \rho = c_s^2. However, you need to use p = c_s^2 \rho itself as an assumption to go between these equalities.

This means that for any other equation of state, where p = c_s^2 \rho is not true, you can’t go from \Delta p / \Delta \rho = c_s^2 (which we saw was a correct approximatino for any eq. of state) to p = c_s^2 \rho. Does that make it clearer?

Yeah, I actually don’t know of any either. What I’m telling you now is stuff that I’ve figured out for myself. For some reason, there has been done very little work on compressible phenomena with LBM compared to how much has been done on incompressible flow.

Erlend M

I’m going to take the risk of confusing things further. I am by no means an expert on this topic so you’ll have to forgive me if I miss the point a bit. Still, here is how I think about (or sometimes think about) the pressure and the equation of state in LB.

There is a linear invertible transformation between the distribution functions and their moments. That is,

m=M*f

where f is the vector of the f_i (so it is 9 dimensional for the D2Q9 lattice), M is the transformation matrix (a linear operator) and the vector m contains the corresponding moments. For example, since the zeroth moment of the distribution function is the density then the first row of M is just 1,1,1,…,1. For a 9 velocity model M will be a 9 by 9 matrix (note that a DmQn lattice as at most n independent moments). The rows of M are combinations of the lattice velocities, and these must be independent, of course.

You think I’m missing the point, don’t you?

We can do the same with the equilibrium (ie with for m^e and f^e) and since the matrix is invertible

f^e=M^-1*m^e

Now we set the equilibrium moments (and M). If we are interested in solving the Navier-Stokes equations then we know that
a) the equilibrium moments of density and velocity are zero because they are conserved (the corresponding rows of M being a row of 1’s, c_ix and c_iy, respectively);
b) the equilibrium momentum flux is \Pi^e_{\alpha\beta}=P+rhou_xu_y, where P is the pressure.

Let’s not worry about the other moments for now as it is these 6 quantities that also appear in the Navier-Stokes equations. So, by having independent rows of M and by specifying the equilibrium moments we can solve for f^e. It is easily checked that the “usual” eqilibrium moments will give you the usual equilibrium function.

The point is, we haven’t said what P is. Indeed, we need to define it to close the system. Most people will have P=rho/3, and for good reason (it’s the natural choice from the low Ma number expansion of the Maxwell-Boltzmann function, it turns out to be the most stable, and anything else for single component flow may mean it’s hard/impossible to obtain a truly time-accurate solution if certain restrictions aren’t met - see below). However, as Erlend says, it is somewhat unphysical…

The ideal gas equation of state (in the LB computational scaling) is P=rho(‘RT’), where ‘RT’ is the gas constant RT divided by c^2 (c is the “particle” speed). The standard choice is ‘RT’=1/3 and gives maximum stability and accuracy. It is unphysical because it is a function of the scaling used (ie c=Delta_x/Delta_t, Delta_ being the space and time step). This isn’t really an issue in the incompressible limit but can be important in, for example, multicomponent systems.

We could, in principle at least, use the true ideal equation of state, P=rho(‘R*T’), without the assumption that ‘RT’ must be 1/3. The resulting equilibrium function would have rho’RT’ terms in it. The problems of having ‘RT’>1/3, for example, can become apparent after the Chapman-Enskog analysis (we’d get a negative bulk viscosity). However if you wanted to simulate, say, hydrogen, in Poiseuille flow at a given (constant) temperate then you could adjust P accordingly. Note that to make sure you’d not have a negative bulk viscosity you’d have to chose your mesh parameters (c, and hence space and time steps) so that ‘RT’ is less that or equal to 1/3.

I hope that hasn’t made thing worse!!

Good luck!

ps
P=c_s^2 is just P=1/3/c^2 for D2Q9. I should also say that the ideal gas law is commonly used for the compressible Navier-Stokes when energy conservation isn’t explicitly declared. Naturally this isn’t perfect but is often sufficient (and the energy equation in LB is another issue)

As far as I know we have
\sum_i f_i^eq = \rho,
\sum_i c_i f_i^eq = \rho u,
which are not zero. Or do you mean something different?

I guess this is one of the problems I have with the LBM. A lot of LBM literature seems to do ‘curve-fitting’, we adjust the physics so it fits in our mathematical framework (Chapman-Enskog expansion in this case). I start to become uneasy when I see things like this and it makes wonder in as much the physics isn’t compromised when people do this (which is why I started this topic). The pressure P is set to c_s^2 \rho, otherwise the Chapman-Enskog expansion does not recover the Navier-Stokes equations. But that it not a good reason for me to say that the pressure should be defined like that, without a proper explanation it seems like ‘curve-fitting’ to me.

Alright, this is how I see it right now:
We have c_s^2 = \partial p \ partial \rho. (Eq. 1) This is pure physics, no approximations or whatever here.
Now we start to approximate \partial p and \partial \rho by linearizing these two terms around a rest state with p0 and \rho0:
\partial p = p - p0
\partial \rho = \rho - \rho0
Substituting this in Eq. 1 we obtain
c_s^2 = (p - p0) / (\rho - \rho0) = \Delta p / \Delta \rho (Eq. 2)
I get this, it does not have anything to do yet with the LBM.

But after that I am no longer following things. Erlend says that we can go from the LB equation of state, c_s^2 = p / \rho, to Eq. 2 but I’m not seeing how. Erlend, could you write that out for me?

Btw, what about all these articles that say “density is related to the pressure with the equation of state of an ideal gas P = c_s^2 \rho” and subsequently simulate liquids? Don’t they know what they are talking about? Are they wrong?

The moments of the equilibrium for density and momentum and the same as the moments of density and momentum with respect to f. In other words, they are conserved by collision. Sorry for the typo - keyboards are not my strong point!

I would say that the LBE has no “extra” physics than the macroscopic equations we are trying to simulate. We prescribe the moments so that we get, for example, the Navier-Stokes. The D2Q9 LBE has 9 moments. 6 of these are the moments that are relevant to the Navier-Stokes flow (density, momentum, monentum flux), the other three are “ghost” moments and should be regarded as numerical features. They do not give us correct kinetic behavior so, if you are using BGK collision, then they should be thought of as error terms. If you use an MRT scheme then they can be adjusted to enhance the stability. The important point, in my opinion, is it is the moments of f_i that are the most relevant, not f_i. The f’s are great for computation (and we can some good advantages because of this), but the moments are far better for the modelling and they can be be defined in a mathematically consistent way so that we approximate (for example) the Navier-Stokes equations to a good degree of accuracy.

Put another way, I would argue that the key feature of LBE is not its ability to capture “more physics” but rather its numerical/mathematical and computational features.

The ideal equation of state is quite fine for most liquids. This is used in classical kinetic theory (see Chapman and Cowling, Grobowski, Maxwell, van Kampen, for example). It is also discussed by Landau and Lifshitz, and they mention the difference between hydrodynamic and thermodynamic pressure. Remember also the the LBE is bound by the small Mach number assumption, so we can not use it for fully compressible flow.

I’d agree that there are LB papers that contain a few dubious things, but the Chapman-Enskog expansion isn’t one of them. This is a very solid piece of theory (at least to second order in Knudsen number) and is again standard in classical kinetic theory. It is an improvement on the Hilbert expansion and has a sound mathematical basis (although there is a question of the validity at higher orders in Kn in terms of objectivity - see the work of Woods from the 1970s), and the maths must be obeyed. Used for the LBE or discrete Boltzmann equation, it will show us how close we are to a set of partial differential equations.

…The Chapman-Enskog expansion will recover the Navier-Stokes with the true equation of state if the equilibrium moments (and hence equilibrium function) is defined so that the pressure is said to be, and left as, P=rho(‘RT’), as suggested above. P=rho/3/c^2 is a particular choice.The question then becomes one of stability and ensuring we do not get a negative bulk viscosity. Still, it is certainly true that the pressure is harder to get right than the momentum. Most LB papers that show convergence (which is very important when assessing the credentials of a numerical method, although not so many do this) will say that pressure is an order lower in accuracy that the velocity (but note the agreement with P=rho/3/c^2 still gives results in decent agreement under mesh refinement, which is perhaps to be expected since the Chapman-Enskog expansion shows that we get the pressure in the Navier-Stokes equations, albeit to one order lower is grid spacing). On the plus side, we get the pressure “directly” without needing a Poison solver. There is no curve fitting because the pressure falls out of the maths of the model. If it is incorrect then we can certainly say, and should say, the model needs to be adjusted. If the choice of P=rho/3/c^2 is the best (which it probably is for incompressibe Navier Stokes) but not perfect then we have to decide if a) it is sufficient for our needs and b) if we can tolerate it (do the advantages we get in return - good accuracy, super-linear scaling etc - out weigh this)?

I suppose I think we should remember we are not just doing physics, we are also doing numerics. Numerical analysis has to be respected as much as the physics. After all, the LBE is a numerical algorithm. As it normally appears its numerical accuracy is tuned very specifically to get the Navier-Stokes equations in the macroscopic incompressible limit (or, alternatively, in the asymptotic limit of the discretisation parameters). If it is to be used for other things, more exotic things, then we should model them and discretise them accordingly. LBE is not always the be the best framework for a given problem. All numercal methods are indirect in the sense that they approximate the true equations in some limit, so they all have errors associated with them. The LBE is no different. Picking the right tools for the right job will mean picking the method with (or getting a balance between) the smallest/most manageable errors and its computational efficiency.

Good luck!

That’s an interesting point you have there, Tim. If I understand it correctly, what you’re saying is that if you know what moments you want the equilibrium distribution to have, you can find the equilibrium distribution itself as

f^e=M^-1*m^e,

instead of finding the equilibrium distribution as an expansion of the Maxwell-Boltzmann distribution. This gives you the freedom to set the pressure to something else than p = (1/3) \rho. However, I’m not quite sure about this statement:

The term “bulk viscosity” is kinda tricky, as there are several different definitions. I prefer to use the one that gives a traceless strain rate tensor, so that monatomic gases have \nu_B = 0 and standard LB has \nu_B = 2 \nu / 3. (Here, \nu and \nu_B are kinematic shear and bulk viscosity, respectively.) With this definition, \nu_B = 0 is the stability limit for fluid mechanics. Still, with ‘RT’ = c_s^2 = 1/3, LB has a positive bulk viscosity, and you should be able to have larger values of ‘RT’ than 1/3 and still have a positive (stable) bulk viscosity. I suppose there must be a stability limit somewhere between ‘RT’ = 1/3 and ‘RT’ = 1, though. (c_s > 1 is obviously unstable from CFL condition considerations.)

Sure! I’ll try to give a clear summary of my previous stuff, too.

So, you’ve already shown in the text I quoted that
c_s^2 = (p - p[sub]0[/sub]) / (\rho - \rho[sub]0[/sub]) = \Delta p / \Delta \rho (Eq. 2)
is a “pseudo-equation-of-state” which is a valid approximation of the definition of the speed of sound, no matter the physical equation of state for this fluid. I call it “pseudo” because the approximation is only valid if \Delta p and \Delta \rho are small, i.e. there is a small variation in pressure and density from the rest state where p = p[sub]0[/sub] and \rho = \rho[sub]0[/sub].

Now, for the LB eq. of state, p = c_s^2 \rho, Eq. 2 is also valid (even exact!):
p[sub]0[/sub] = c_s^2 \rho[sub]0[/sub]
=> p - p[sub]0[/sub] = c_s^2 ( \rho - \rho[sub]0[/sub] ) or \Delta p = c_s^2 \Delta \rho
=> c_s^2 = \Delta p / \Delta \rho. (Eq. 2)

The result of this is that as long as the variations in pressure and density in the simulation are not that large, you can use LB to emulate a fluid with any equation of state if you’re a bit careful with the scaling. And you have to take into account that p[sub]0[/sub] / \rho[sub]0[/sub] will not be the same in the LB simulation and the physical fluid.

I guess this is comparable to doing experiments in one fluid (e.g. air) and using them to say something about the corresponding behaviour in another fluid with another equation of state (e.g. water). For instance, sound propagates in just the same way in air and water, although the two fluids have a different speed of sound. Also, water is less compressible, which means that a sound wave with the same peak pressure in air and water will have a lower peak density in water than in air. This difference can be seen through good ol’ Eq. 1,

c_s^2 = \partial p \ partial \rho,

as water has a significantly higher speed of sound (about 1500 m/s) than air (about 343 m/s).

However, since the two fluids have different equations of state, they will have different points at which acoustic nonlinearities become important. Nonlinearities occur when the sound wave amplitude becomes really large, i.e. when

\partial p / partial \rho = \Delta p / \Delta \rho

is no longer a good approximation. (There is another source of nonlinearity too, but I’ll spare you by not going into that. :slight_smile: )

I hope that helps!
Erlend M

Yep, that’s what I am saying. The small Mach number assumption still applies as we must enforce U<<c (which effectively amounts to the CFL condition). Prescribing the moments is very similar to defining the equilibrium (it must be - we have a linear transformation between the two) but the reason I (sometimes) like to think about it in this way is we know we can only ever have 9 (independent) moments (for D2Q9) while we only actually care about 6 for Navier-Stokes. The third rank tensor usually called Q=sum{f_ic_ic_i*c_i} is underdetermined - we do not have enough lattice vectors to define all of its components. Luckily we only need the equilibrium values of Q when trying to get Navier-Stokes form the LBE via Chapman-Enskog (incidentally, this gives an Ma^3 error to the Navier-Stokes which can, I think, be eliminated with a larger stencil). This means we have (unphysical) ghost moments which we need to get to grips with. I find that using moments gives a more natural (or more mathematical, perhaps - it has the solid theoretical grounds of linear algebra) way of seeing this and understanding (hopefully eventually improving) the model. This is one of the main motivations behind various MRT schemes.

I completely agree with you about the bulk viscosity and how it should be defined by the traceless tensor. I was being very sloppy above but I meant that if to conduct the Chapman-Enskog expansion on the LBE with ‘RT’ in the equilibrium function (or the moments) then you’ll get a div u type term in the momentum equation that has a coefficient of (1/3-‘RT’).

Perhaps the key thing to remember in all of this is, like Erlend says, we are restricted by the small Ma number assumption so can not possibly hope to do fully compressible flows, and so the LB equation of state is often sufficient. Without considerable changes to the method, the LBE will not produce reliable results for large Ma, just like it won’t give reliable results for large Kn numbers. If we are to get the best out of LBE we must know its weaknesses as well as its strengths.

Thanks Erlend and pleb01, I think I get it now. p0 = c_s^2 \rho0 is valid in the rest state, and c_s^2 = \Delta p / \Delta \rho remains valid in non-rest states, as long as we do not deviate too much from the rest state.

Also thanks for the additional info regarding this topic.