Not clear which version of Shan Chen model is implemented

[Jendrik Weise]

Background:
———–
The original ShanChen93 paper proposed and implemented the forcing within the Shan-Chen 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 Rayleigh-Taylor testcase, however it may still be important for more rigorously measured numerical accuracy.
2. In ll. 122-125 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 Rayleigh-Taylor test case.

[Jendrik Weise]

I 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 Shan-Chen 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.

[Jendrik Weise]

Having 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 5 months, 3 weeks ago by Jendrik Weise.

[TimBingert]

Hi Jendrik,

there is an actual issue with the rayleighTaylor cases as they use the Shan-Chen post processor, but they use dynamics with the Guo forcing scheme. The model that is meant to be iplemented is the multi-component Shan-Chen model from the Krüger book, where you have the Shan-Chen 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
Tim

[Jendrik Weise]

Hi,
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

[TimBingert]

You 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 Shan-Chen-like 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

[TimBingert]

Excuse 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.

[danial]

Hi, 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 Shan-Chen 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 uses rhoTot = 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 for waterAirflatInterface2d, 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.

[TimBingert]

Hi 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

[Author]

Posts