Help with the velocity shift in ShanChen

[lpaiva]

Hi all!

I am attempting to model some higher density ratio simulations using ShanChen and was looking for a way to test some different solutions for it using e.g. different force models for solid-fluid interaction in the velocity shift equilibrium equation. Although I have seen that the velocity shift was already implemented, I am having trouble understanding where the force comes from and where I could for example write a function to implement these other different force equations.

So far, what I was able to gather was that the velocity shift is implemented in “forcing.h” through

CellStatistic<V> compute(CELL& cell, PARAMETERS& parameters, FEQ& fEq) any_platform {
V rho, u[DESCRIPTOR::d], uSqr{};
MomentaF().computeRhoU(cell, rho, u);
const auto force = cell.template getFieldPointer<descriptors::FORCE>();
for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) {
u[iVel] += force[iVel] / parameters.template get<descriptors::OMEGA>();
uSqr += u[iVel] * u[iVel];
}
for (unsigned iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
fEq[iPop] = EquilibriumF().compute(iPop, rho, u);
}
return {rho, uSqr};
};

but when I try to follow the rabbit hole of finding where the force is specifically defined, I only find it defined as a generic FIELD_BASE, not having quite understood where the actual equation comes in.

Thank you,

Lucas

[stephan]

Dear Lucas,

thanks for posting.
Please note that in general, forces are realized in OpenLB as fields, which probably is the reason for you ending up in the generic FIELD_BASE.
The fields act as cell storage and can be accessed via e.g. getFieldPointer in the code snippet you pasted.

Long story short, to implement different forces you could hardcode them in your example/application .cpp file.
By that, I mean specify the force as a function and then pass it to the field storage.
This is done e.g. in examples/multiComponent/rayleighTaylor3d/rayleighTaylor3d.cpp, where an external force is added via an additional field

using DESCRIPTOR = D3Q19<FORCE,EXTERNAL_FORCE,STATISTIC>;

which is in turn explicitly specified later via

sLatticeOne.defineField<descriptors::EXTERNAL_FORCE>( superGeometry, 1, f );

and hence can be used in the forcing scheme through accessing the EXTERNAL_FORCE field.

BR
Stephan

[GASciences]

Hi Lucas,

I just want to provide additional comments, mostly for other answers seekers like myself.

You said in you question that you are looking for a way to change the interaction force’s formulation and that you looked into the velocity shift implementation and still had no clue on how to do it.

It is implicit in Stephan’s answer, but the velocity shift is not where the interaction force is calculated, it is rather where this force is introduced into the collide step, no matter the equation you use to calculate it. As Stephan said, you rather want to use another postProcessor to calculate it, look for example ShanChenForcedSingleComponentPostProcessor which does (used in examples//multiComponent/phaseSeparation2d/phaseSeparation2d.cpp).

Finally, you said you are interesting in implementing your own solid/liquid interaction force’s formulation? Then maybe considering a boundary condition is what you need ? In that case, have you looked at this answer https://www.openlb.net/forum/topic/wetting-with-free-surface-or-large-density-ratio/#post-7308 ?

Cheers,

Gilles

[Author]

Posts