Velocity based error in secondary lattice

[Julius]

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[1] = omega_ * cell[1] + omega_2 * ( rho_ - jx );
Where omega_ = (1-omega), omega_2 = omega / 2, rho_ = rho - 1 and jx = rho * u[0]. Substituting those, one gets:
cell[1] = (1-omega) * cell[1] + omega * 1/2 * ( rho_*cs2 - rho * u[0] )
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] = (1-omega) * cell[1] + omega * w / cs2 * ( rho_*cs2 - rho * u[0] )
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] = (1-omega) * cell[1] + omega * w / cs2 * ( T_0*cs2 - T_1 * u_1[0] )
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[0]
cell[1] = (1-omega) * cell[1] + omega * w / cs2 * T_0 * ( cs2 - u_0[0]  )
cell[1] = (1-omega) * cell[1] + omega * w / cs2 * ( T_0 * cs2 - T_0 * u_0[0]   )
Now let’s compare the two:
OpenLB: cell[1] = (1-omega) * cell[1] + omega * w / cs2 * ( T_0*cs2 - T_1 * u_1[0] )
BGK : cell[1] = (1-omega) * cell[1] + omega * w / cs2 * ( T_0 * cs2 - T_0 * u_0[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 5 years, 4 months ago by Julius.
• This topic was modified 5 years, 4 months ago by Julius. Reason: Formatting code

[mgaedtke]

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

[Julius]

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

[Author]

Posts