Not clear which version of Shan Chen model is implemented
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › Bug Reports › Not clear which version of Shan Chen model is implemented
 This topic has 8 replies, 3 voices, and was last updated 5 months ago by TimBingert.

AuthorPosts

November 6, 2023 at 3:28 pm #7897Jendrik WeiseParticipant
Background:
———–
The original ShanChen93 paper proposed and implemented the forcing within the ShanChen model via the use of a forcing scheme based solely on altering the velocity. As this was later shown to pose problems where large viscosity ratios are concerned, Porter12 proposed the use of a more standard forcing scheme in combination with a only slightly altered equilibrium velocity calculation in order to alleviate such issues. Most commonly these days, the forcing scheme chosen for it is, unlike the original paper which used He forcing, Guo forcing.Potential issues:
—————–
1. The current implementation within OpenLB uses Guo forcing, while not adding the force/2 term from Porter12 to the phase momenta. Changing this does not produce any visible changes, at least in the RayleighTaylor testcase, however it may still be important for more rigorously measured numerical accuracy.
2. In ll. 122125 of shanChenForcedPostProcessor.h, the force contributions from the other fluids are, in addition to being multiplied by G, also divided by rhoField[0] and rhoField[1], respectively. I can’t find any references that would motivate this operation. Removing it causes the phase interface to be more diffuse, but the simulation results, at least visually, look most in line with what for instance Palabos(which implements the traditional version of ShanChen93) produces for its RayleighTaylor test case.November 6, 2023 at 4:44 pm #7901Jendrik WeiseParticipantI believe the issue with 1. may be the following: internally(elements.h#1095), the Guo forcing adds the term force/2 for each lattice individually, however normally this should happen before the averaging(at least for the ShanChen forces), meaning that in a correct explicit forcing implementation all fluids at a given location have the same equilibrium velocity. I have yet to address that potential issue in the code myself so take the above results with a grain of salt.
November 7, 2023 at 10:47 am #7903Jendrik WeiseParticipantHaving investigated this further, point 2 is merely a result of “force” being actually an acceleration and simply confusingly named. Having now tested 1, correct results are indeed obtained by switching to dynamics along the lines of:
template <typename MOMENTA>
struct Identity {
template<typename DESCRIPTOR>
using type = typename MOMENTA::template type<DESCRIPTOR>;
};template <typename T, typename DESCRIPTOR, typename MOMENTA=momenta::BulkTuple>
using IdentityForcedBGKdynamics = dynamics::Tuple<
T, DESCRIPTOR,
MOMENTA,
equilibria::SecondOrder,
collision::BGK,
forcing::Guo<Identity>
>;
and adding the forcing momentum adjustments within the processor. I believe altering the default behavior here and documenting this in more detail may indeed by advisable. This reply was modified 10 months, 1 week ago by Jendrik Weise.
November 8, 2023 at 4:00 pm #7913TimBingertParticipantHi Jendrik,
there is an actual issue with the rayleighTaylor cases as they use the ShanChen post processor, but they use dynamics with the Guo forcing scheme. The model that is meant to be iplemented is the multicomponent ShanChen model from the Krüger book, where you have the ShanChen forcing scheme, see equations 9.125 and 9.126. So one should change the dynamics from “ForcedBGKdynamics” to “ForcedShanChenBGKdynamics” that use the velocity shifted equilibrium, similar to the phase separation cases. This will match the post processor, as there is no addition of F/2 to the velocity in this model which is what I suppose you tried to do with these IdentityForcedBGKdynamics?!
Hope this helps
Best
TimNovember 8, 2023 at 4:08 pm #7914Jendrik WeiseParticipantHi,
that is another potential resolution. I had assumed that, given that Guo forcing was chosen, implementing something along the lines of https://journals.aps.org/pre/abstract/10.1103/PhysRevE.86.036701 had been the intention. To do that correctly though, in accordance with formulas (15) and (16) in that paper, the force/2 addition needs to occur before the averaging step.In any case though, if there was no intention to implement that model, it would likely make sense to simply switch to a ShanChen forcing scheme to avoid future confusions, which would bring the model in line with the original Shan Chen paper.
Jendrik Weise
November 8, 2023 at 5:02 pm #7916TimBingertParticipantYou are correct, the Guo forcing was chosen which does not work in this case as the coupling strictly requires Shan Chen forcing as there is no force to add after the averaging of velocities.
In the next OpenLB release, there will probably be a ShanChenlike model that uses Guo forcing instead of Shan Chen forcing. For now, it would be easiest to change the dynamics from “ForcedBGKdynamics” to “ForcedShanChenBGKdynamics” in the Rayleigh Taylor cases.Tim
November 8, 2023 at 5:59 pm #7917TimBingertParticipantExcuse me, I have to make another comment on your case: There is a reason why the Rayleigh Taylor cases used a Guo forcing which is the existence of another body force as gravity. But as the Shan Chen forcing scheme is also a general forcing scheme, see page 241 in the Krüger book, the gravity will be considered in the velocity shifted equilibrium as the external force is being added to the Shan Chen force as a final step of the post processor.
April 1, 2024 at 7:29 am #8435danialParticipantHi, in the ShanChanForcedProcessor, the part that divides the force by density
cellB.template setField<descriptors::FORCE>(externalPartnerForce  g*rhoBlockContribution/rhoField[1]);
I think the division by density is not needed. If we wanted to directly update the velocity like the ShanChen paper that uses u_phase=u_total+tau*F_ext/rho, yes the division by rho is necessary (Kruger p241, p378), but since we are simply passing the force to whatever forcing scheme is being used (Guo or ShanChen) then this division isn’t necessary. Yes the interface will become more diffuse but the critical value of G=4 (Kruger P375) must be used instead of 3 (I have to do some tests to verify this).
The new
MCMPForcedPostProcessor
is including the external force in the shared velocity
uTot = (u_sum + totalForce/2.) / rhoTot + externalForce1/2
directly but the old ShanChen just usesrhoTot = rhoField[0]*omega1+ rhoField[1]*omega2
but that one is based on a newer modern paper so I was trying to see how the new example works, but I don’t see any change during the simulation for the water column forwaterAirflatInterface2d
, I must have messed up a setting somewhere but I’ll have to double check.Bytheway thank you guys for developing OpenLB and sharing it for free and providing this forum so I can also see the issues encountered by others. Cheers.
April 12, 2024 at 11:11 am #8458TimBingertParticipantHi danial,
the division by the density is always needed, as it is once again multiplied by in the forcing scheme.
Please note, that the pseudopotential models as you mentioned are from different works and hence quite different. The old ShanChen model is using a ShanChen forcing scheme so it shows the relaxation frequencies omega at certain points. This makes the equilibrium depend on the relaxation time and therefore on viscosity which can be problematic as well as the old model only provides low density ratios.
The model used in MCMPForcedPostProcessor is using a Guo forcing scheme and has conservation laws to determine total densities and shared velocities.Kind regards
Tim 
AuthorPosts
 You must be logged in to reply to this topic.