OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB Bug Reports Velocity based error in secondary lattice

Viewing 3 posts - 1 through 3 (of 3 total)
• Author
Posts
• #5038
Julius
Participant

Dear developers,

I’ve been trying to get my thermal multiphase flow to work, but the thermal lattice always came up with large errors. After digging through the source code, I’ve found where the origin of it lies.

This bug is relevant to any simulation where velocity is passed from the fluid lattice to another lattice.

I’ll walk through the code iteration step by step to explain where the error is coming from. To couple the thermal lattice to the fluid lattice, a post processor is needed, such as `NavierStokesAdvectionDiffusionCouplingGenerator3D`.

In the step of `executeCoupling()`, in the `navierStokesAdvectionDiffusionCouplingPostProcessor` , the velocity is passed on like this:

``````T *u = tPartner->get(iX,iY,iZ).template getFieldPointer<descriptors::VELOCITY>();
blockLattice.get(iX,iY,iZ).computeU(u);``````

The `blockLattice` calls `dynamics::computeU` which calls `dynamics::computeRhoU` which calls `Momenta::computeRhoU` which calls `lbHelpers::computeRhoU`. There the definition of the density is the sum of all cell population plus 1. I’ll name this `rho_1`. I’ll notate the density without the +1 as `rho_0`. Thus `rho_1 = rho_0 + 1`
Then the velocity is calculated in OLB via the sum of direction * cell population, which is `j`, and then divided by `rho_1`. Thus `u_1 = j / rho_1`. In literature, the velocity is calculated as follows: `u_0 = j / rho_0`

The next time `collideAndStream()` is called for the thermal lattice, `lbHelpers::bgkCollision` is called from within `AdvectionDiffusionBGKdynamics::collide()` with the arguments of cell, temperature, the computed velocity, and omega.

The temperature is the density of the thermal lattice, and is computed the same as with the fluid. Using the same notation as with density, I’ll write the computed temperature as `T_1` and `T_1 = T_0 +1`.

In the BGK collision for the thermal lattice, I’ll use D3Q7 as an example, the non-rest velocity directions are calculated as follows, using the 1st direction as example:
`cell = omega_ * cell + omega_2 * ( rho_ - jx );`
Where `omega_ = (1-omega)`, `omega_2 = omega / 2`, `rho_ = rho - 1` and `jx = rho * u`. Substituting those, one gets:
`cell = (1-omega) * cell + omega * 1/2 * ( rho_*cs2 - rho * u )`
In the case of D3Q7 cs2 is 1/4 and the weights for non-rest directions is 1/8, thus the weight can also be said to be cs2/2, or `1/2 = w / cs2`:
`cell = (1-omega) * cell + omega * w / cs2 * ( rho_*cs2 - rho * u )`
Since this bgk collision is called via the thermal lattice, the `rho` in the above is `T_1` and `u` is the computed velocity `u_1`
`cell = (1-omega) * cell + omega * w / cs2 * ( T_0*cs2 - T_1 * u_1 )`
The above equation is the bgk implementation of OLB for thermal lattice. This can be found for example in `rayleighBenard`.

Starting from the BGK step, with first order expansion for the equilibrium population, I arrive to the below equation:
`cell[i] = (1-omega) * cell[i] + omega * w * T_0 * ( 1 + e_i * u_0 / cs2  )`
when factoring out cs2
`cell[i] = (1-omega) * cell[i] + omega * w / cs2 * T_0 * ( cs2 + e_i * u_0  )`
and for `i=1` the following holds `e_i * u_0 = -1 * u_0`
`cell = (1-omega) * cell + omega * w / cs2 * T_0 * ( cs2 - u_0  )`
`cell = (1-omega) * cell + omega * w / cs2 * ( T_0 * cs2 - T_0 * u_0   )`
Now let’s compare the two:
OpenLB: `cell = (1-omega) * cell + omega * w / cs2 * ( T_0*cs2 - T_1 * u_1 )`
BGK : `cell = (1-omega) * cell + omega * w / cs2 * ( T_0 * cs2 - T_0 * u_0   )`

And here is where the bug is: For the fluid lattice `rho_1 * u_1` is equal to `rho_0 * u _0` since `rho_1 * u_1 = rho_1 * (j / rho_1) = j = rho_0 * (j / rho_0) = rho_0 * u _0`
However, since we now are in the thermal lattice, `T_0 * u_0 =/= T_1 * u_1`
since `T_0 * u_0 = T_0 * (j / rho_0) =/= T_1 * (j / rho_1) = rho_1 * u _1`
Thus the simulations do not give the expected results.

This can be ‘fixed’ in two steps:
In the coupling `navierStokesAdvectionDiffusionCouplingPostProcessor` , the velocity should be passed on like this (note the usage of `computeJ` instead of `computeU`):

``````T *u = tPartner->get(iX,iY,iZ).template getFieldPointer<descriptors::VELOCITY>();
blockLattice.get(iX,iY,iZ).computeJ(u);
T localRho = blockLattice.get(iX,iY,iZ).computeRho();
for (int iD=0; iD< DESCRIPTOR::d; ++iD){
u[iD] /= (localRho-1) ;
}``````

And in the `bgkCollision` in the helper file, the lines computing jx, jy, jz should use `rho_` instead of `rho`.
The error that is produced is dependent heavily on the density, and thus becomes quite apparent in multiphase flow.

Regards, Julius

• This topic was modified 2 years, 11 months ago by Julius.
• This topic was modified 2 years, 11 months ago by Julius. Reason: Formatting code
#5047
mgaedtke
Keymaster

Hey Julius,

I can’t reproduce your errors. Please make sure to understand the return values of our cell objects: The distribution functions we save in the data structure are actually f* = f – weight, so that the sum of f* = rho – 1, because the sum of the weights equals to 1. We apply the same shift to all the equilibrium distributions.

See Section V in https://journals.aps.org/pre/abstract/10.1103/PhysRevE.48.4823 for some more information on this topic.

Cheers,
Max

#5048
Julius
Participant

Hello Max,

Thank you for looking into my bug report.
Your comment regarding how the population is stored into the cell values was very helpful. And indeed, with that in mind, my bug report was erroneous.

Since I am relatively new to the LBM scene, I was not aware that such measures are taken to reduce the numerical roundoff error. Perhaps it would be nice to make that clear in the documentation. Technically, it is written in there on page 49 (olb_ug-1.3r0), but only in one sentence, which can be easily read over.

Thank you for your time and help!

Regards, Julius

Viewing 3 posts - 1 through 3 (of 3 total)
• You must be logged in to reply to this topic.