Question about the IB force f from version 2.2

[AaaaronLee]

Dear developers,

Thank you for sharing the code.

I noticed that the IB force f has changed from version 2.1 to version 2.2 as follows:

Version 2.1

f = surface*(UA*Ui - UEqn.H() + fvc::grad(p));

Version 2.2

f += 0.8*surface*(Ui - U)/runTime.deltaT();

The previous version made sense to me, as the IB force cancels out the external force and maintains the IB condition.

Would you mind explaining the mechanism behind the new implementation?

I suspect it may act as a damping force based on the current velocity difference. However, such a scheme could potentially lead to sinusoidal velocity or IB force responses.

Does it still contribute stably to the IB?

It seems that the pEqn.H file is the only one manipulating the f.

Please correct me if I'm missing some code about f or have misunderstood the changes.

Thank you for your time.

Best regards,
Aaron

[techMathGroup]

Hi, we switched to the formulation used for example in Blais et al. (2016) [doi: 10.1016/j.compchemeng.2015.10.019], see equation (13). In transient cases, this proved to be slightly more stable (due to the under-relaxation). However, in steady-state, we still use the original definition, see e.g. https://github.com/techMathGroup/openHFDIBRANS/blob/main/applications/solvers/simpleHFDIBRANSFoam/UEqn.H and we are repeatedly re-evaluating and comparing the formulations. Altogether, at convergence to steady-state, these provide almost identical results (as can be expected). In transient and CFD-DEM, the newer formulation is easier to handle when it comes to potential updating of the forcing term during DEM loop (but this is not at the moment active in the public version as it is required only by "ugly" physical scenarios).

Best,
Martin