Solid–Fluid Interaction

Overview

The InteractionSolidFluid class handles the coupling between fluid (LBM) and solid (phase-field) regions. It computes the velocities of solid particles induced by fluid forces and distributes the corresponding solid-body velocity field back into the Velocities object, which is then used by the advection and LBM modules.

Defined in openphase/include/FluidDynamics/InteractionSolidFluid.h.

---

Class declaration

struct InteractionSolidFluid : public OPObject
{
    InteractionSolidFluid(Settings& locSettings,
                          std::string InputFileName = DefaultInputFileName);

    void Initialize(Settings& locSettings,
                    std::string ObjectNameSuffix = "") override;
    void ReadInput(std::string InputFileName) override;

    void CalculateSolidVelocities(PhaseField& Phase, Velocities& Vel,
                                  const BoundaryConditions& BC,
                                  const double dt) const;

    static dVector3 LocalGrainVelocity(long i, long j, long k,
                                       const Grain& grain,
                                       const BoundaryConditions& BC,
                                       const PhaseField& Phase);
};

---

Bounce-back boundary condition

The fundamental solid–fluid boundary condition in LBM is bounce-back: a population that would stream into a solid node is reflected back in the opposite direction.

Standard bounce-back

If node (i,j,k) is solid:
    f_{-x,-y,-z}(i,j,k) ← f_{x,y,z}(i,j,k)

This implements a no-slip (zero velocity) condition at the solid surface. In OpenPhase the bounce-back is applied in FlowSolverLBM::BounceBack() and Propagation().

Elastic bounce-back

When Do_BounceBackElastic = true, populations are reflected using the moving boundary correction:

$$\begin{array}{}\text{(1)}& f_{-x,-y,-z}(\mathbf{x},t+1)=f_{xyz}(\mathbf{x},t)-2w_{xyz}\frac{{\mathbf{e}}_{xyz}\cdot {\mathbf{u}}_w}{c_s^2}\rho \end{array}$$

where ${\mathbf{u}}_w$ is the solid wall velocity. This is energy-conserving and essential for mobile solid particles.

---

Drag force

When Do_Drag = true, the solver computes a viscous drag force at the diffuse solid–fluid interface:

void CalculateForceDrag(int i, int j, int k,
                        PhaseField& Phase,
                        const Velocities& Vel);

The drag contribution to the force density is proportional to the velocity difference between the fluid and the solid:

$$\begin{array}{}\text{(2)}& {\mathbf{F}}_{\text{drag}}(\mathbf{x})=h^{\ast }{\varphi }_s(\mathbf{x})[{\mathbf{u}}_s(\mathbf{x})-{\mathbf{u}}_f(\mathbf{x})]\end{array}$$

where $h^{\ast }$ is the dimensionless drag parameter (h_star), ${\varphi }_s$ is the solid fraction, ${\mathbf{u}}_s$ is the solid velocity, and ${\mathbf{u}}_f$ is the fluid velocity.

---

Rigid body dynamics

InteractionSolidFluid computes translational and rotational velocities of solid particles using Newton's law for rigid bodies, driven by the forces that the fluid exerts on them.

With inertia (DoInertia = true)

The grain velocity is updated by integrating:

$$\begin{array}{}\text{(3)}& m\dot{\mathbf{v}}={\mathbf{F}}_{\text{fluid}}\end{array}$$

$$\begin{array}{}\text{(4)}& \mathbf{I}\dot{\mathit{\omega }}={\mathit{\tau }}_{\text{fluid}}\end{array}$$

where $m$ is the grain mass, $\mathbf{I}$ is the moment of inertia, ${\mathbf{F}}_{\text{fluid}}$ is the net fluid force, and ${\mathit{\tau }}_{\text{fluid}}$ is the fluid torque.

Without inertia (DoInertia = false)

A damped (quasi-static) model is used instead:

$$\begin{array}{}\text{(5)}& \mathbf{v}={\tau }_{\text{trans}}{\mathbf{F}}_{\text{fluid}}\end{array}$$

$$\begin{array}{}\text{(6)}& \mathit{\omega }={\tau }_{\text{rot}}{\mathit{\tau }}_{\text{fluid}}\end{array}$$

where ${\tau }_{\text{trans}}$ and ${\tau }_{\text{rot}}$ are the translational and rotational relaxation parameters.

Velocity limit

When DoLimitVelocity = true, solid body velocities are clipped to VelocityLimit to prevent numerical blow-up when solid particles collide.

---

Grain geometry helpers

Internally, the class computes center-of-mass and moment of inertia for each grain using domain-summed integrals:

static void CalculateCenterOfMass(PhaseField& Phase,
                                   const BoundaryConditions& BC);
static void CalculateMomentOfInertia(PhaseField& Phase,
                                      const BoundaryConditions& BC);

These computations include OpenMP and MPI reduction, making them efficient on HPC systems.

---

Obstacle state tracking

FlowSolverLBM tracks changes in obstacle nodes to correctly conserve fluid mass:

State	Meaning
Absent	Fluid node
Present	Solid node
Appeared	Newly solidified (fluid → solid)
Vanished	Newly melted (solid → fluid)
ChangedDensity	Solid changed density (e.g. due to composition change)

When a solid node appears, the fluid it displaced must be redistributed. When a solid node vanishes, new populations are initialized from the neighboring fluid state via SetFluidNodesNearObstacle().

---

Parameters

Parameter	Type	Description
TauTranslation	double	Translational relaxation time (no-inertia model)
TauRotation	double	Rotational relaxation time (no-inertia model)
VelocityLimit	double	Maximum solid body velocity
DoInertia	bool	Enable inertial rigid body dynamics
DoRotations	bool	Enable rotational degrees of freedom
DoLimitVelocity	bool	Enable velocity clamping
DoEnforceZeroTotalMomentum	bool	Conserve total linear momentum
DoEnforceZeroTotalAngularMomentum	bool	Conserve total angular momentum

---

Typical usage

InteractionSolidFluid ISF(OPSettings, InputFile);

// each timestep:
ISF.CalculateSolidVelocities(Phase, Vel, BC, dt);
FL.Solve(Phase, Vel, BC);

CalculateSolidVelocities must be called before FlowSolverLBM::Solve() so that the solid velocities are available for the elastic bounce-back computation.

---

Numerical considerations

• For sintering simulations where grains continuously merge and split, the grain tracking inside InteractionSolidFluid relies on stable grain IDs from the PhaseField module.
• The moment of inertia calculation assumes isotropic grain density; anisotropic grains may require a custom implementation.
• DoEnforceZeroTotalMomentum prevents net drift of the solid ensemble in periodic domains but should not be combined with an imposed flow.