Velocity based error in secondary lattice
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › Bug Reports › Velocity based error in secondary lattice
 This topic has 2 replies, 2 voices, and was last updated 4 years ago by Julius.

AuthorPosts

June 24, 2020 at 2:32 pm #5038JuliusParticipant
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 thenavierStokesAdvectionDiffusionCouplingPostProcessor
, 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
callsdynamics::computeU
which callsdynamics::computeRhoU
which callsMomenta::computeRhoU
which callslbHelpers::computeRhoU
. There the definition of the density is the sum of all cell population plus 1. I’ll name thisrho_1
. I’ll notate the density without the +1 asrho_0
. Thusrho_1 = rho_0 + 1
Then the velocity is calculated in OLB via the sum of direction * cell population, which isj
, and then divided byrho_1
. Thusu_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 withinAdvectionDiffusionBGKdynamics::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
andT_1 = T_0 +1
.In the BGK collision for the thermal lattice, I’ll use D3Q7 as an example, the nonrest velocity directions are calculated as follows, using the 1st direction as example:
cell[1] = omega_ * cell[1] + omega_2 * ( rho_  jx );
Whereomega_ = (1omega)
,omega_2 = omega / 2
,rho_ = rho  1
andjx = rho * u[0]
. Substituting those, one gets:
cell[1] = (1omega) * cell[1] + omega * 1/2 * ( rho_*cs2  rho * u[0] )
In the case of D3Q7 cs2 is 1/4 and the weights for nonrest directions is 1/8, thus the weight can also be said to be cs2/2, or1/2 = w / cs2
:
cell[1] = (1omega) * cell[1] + omega * w / cs2 * ( rho_*cs2  rho * u[0] )
Since this bgk collision is called via the thermal lattice, therho
in the above isT_1
andu
is the computed velocityu_1
cell[1] = (1omega) * 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 inrayleighBenard
.Starting from the BGK step, with first order expansion for the equilibrium population, I arrive to the below equation:
cell[i] = (1omega) * cell[i] + omega * w * T_0 * ( 1 + e_i * u_0 / cs2 )
when factoring out cs2
cell[i] = (1omega) * cell[i] + omega * w / cs2 * T_0 * ( cs2 + e_i * u_0 )
and fori=1
the following holdse_i * u_0 = 1 * u_0[0]
cell[1] = (1omega) * cell[1] + omega * w / cs2 * T_0 * ( cs2  u_0[0] )
cell[1] = (1omega) * cell[1] + omega * w / cs2 * ( T_0 * cs2  T_0 * u_0[0] )
Now let’s compare the two:
OpenLB:cell[1] = (1omega) * cell[1] + omega * w / cs2 * ( T_0*cs2  T_1 * u_1[0] )
BGK :cell[1] = (1omega) * 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 torho_0 * u _0
sincerho_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
sinceT_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 couplingnavierStokesAdvectionDiffusionCouplingPostProcessor
, the velocity should be passed on like this (note the usage ofcomputeJ
instead ofcomputeU
):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] /= (localRho1) ; }
And in the
bgkCollision
in the helper file, the lines computing jx, jy, jz should userho_
instead ofrho
.
The error that is produced is dependent heavily on the density, and thus becomes quite apparent in multiphase flow.Regards, Julius
June 30, 2020 at 9:39 am #5047mgaedtkeKeymasterHey 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,
MaxJune 30, 2020 at 3:33 pm #5048JuliusParticipantHello 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_ug1.3r0), but only in one sentence, which can be easily read over.
Thank you for your time and help!
Regards, Julius

AuthorPosts
 You must be logged in to reply to this topic.