FlowSolverLBM

Overview

FlowSolverLBM is the central class of the FluidDynamics module. It implements a D3Q27 BGK Lattice Boltzmann solver for single- and multi-component fluid flow on the OpenPhase structured grid. It is the only class a typical user needs to instantiate directly; the interaction classes are optional extensions.

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

---

Class declaration

class FlowSolverLBM : public OPObject
{
public:
    FlowSolverLBM(Settings& locSettings,
                  const std::string InputFileName = DefaultInputFileName);

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

    void Solve(PhaseField& Phase, Velocities& Vel,
               const BoundaryConditions& BC);
    void Solve(PhaseField& Phase, const Composition& Cx, Velocities& Vel,
               const BoundaryConditions& BC);
    // ...
};

FlowSolverLBM inherits from OPObject, giving it standard initialization, input-reading, restart, and VTK output capabilities.

---

Data layout

Population arrays

Populations are stored as Storage3D<D3Q27, 1> — a 3D grid where each node holds a D3Q27 object containing 27 double values:

Storage3D<D3Q27, 1> lbPopulations;    ///< Current populations f_i
Storage3D<D3Q27, 1> lbPopulationsTMP; ///< Buffer for propagation step

The D3Q27 class indexes populations with integer offsets $(x,y,z)\in \{-1,0,1\}$:

double& f = lbPopulations(i, j, k)(-1, 0, 1);  // population at e = (-1,0,1)

Internal flat index: $\text{idx}(x,y,z)=13+9x+3y+z$.

Obstacle arrays

Storage3D<ObstacleStates, 0> ObstacleState;  ///< Per-node obstacle state
Storage3D<int, 0>            Obstacle;       ///< 1 if node is solid

ObstacleStates tracks changes: Absent, Present, Appeared, Vanished, ChangedDensity.

Physical field arrays

Storage3D<dVector3, 1> ForceDensity;   ///< Body force density [lb units]
Storage3D<dVector3, 1> MomentumDensity;
Storage3D<double,   1> DensityWetting; ///< Component densities / wetting param

---

Key parameters

Per-component parameters

Member	Type	Description
nu	vector<double>	Kinematic viscosity [m²/s]
lbtau	vector<double>	BGK relaxation time $\tau$
LiquidDensity	vector<double>	Equilibrium liquid density [kg/m³]
VaporDensity	vector<double>	Equilibrium vapor density [kg/m³]
SurfaceTension	vector<double>	Surface tension [kg/s²]
rho_0	vector<double>	Reference density for Benzi force
GasParameter	vector<double>	Kupershtokh gas parameter

Scaling coefficients

Member	Description
dRho	Physical density per lattice density [kg/m³]
dP	Physical pressure per lattice pressure [Pa]
dm	Physical momentum density per lattice value
df	Physical force density per lattice value
dnu	Physical kinematic viscosity per lattice viscosity [m²/s]
dt	Physical time per lattice time step [s]

Feature flags

Flag	Effect
Do_Benzi	Enable Benzi pseudo-potential force
Do_Kupershtokh	Enable Kupershtokh two-phase model
Do_TwoPhase	Two-phase flow with phase separation
Do_BounceBack	Bounce-back solid–fluid BC
Do_BounceBackElastic	Energy-conserving elastic bounce-back
Do_Drag	Surface drag force at interfaces
Do_Gravity	Gravitational body force
Do_GuoForcing	Use Guo force-incorporation scheme
Do_EDForcing	Use exact difference method for forces
Do_FixPopulations	Correct negative populations
Do_ThermalComp	Enable thermal compressibility extension
Do_DI	Treat the interface as diffuse

---

Solver workflow

The high-level per-timestep algorithm called by Solve():

1. DetectObstacles(Phase)              — update solid/fluid node map
2. ApplyForces(Phase, Vel)             — compute body forces (gravity, Benzi, ...)
3. Collision()                         — BGK relaxation
4. Propagation(Phase, BC)             — streaming
5. SetBoundaryConditions(BC)          — apply inlet/outlet/periodic BCs
6. CalculateDensityAndMomentum()      — compute ρ, ρu from populations
7. CalculateFluidVelocities(Vel, ...) — convert to physical velocities

AccountForAndLimitPhaseTransformation() is additionally called when solid/fluid phase changes occur to conserve mass.

---

Key methods

Initialization

void Initialize(Settings& locSettings, std::string ObjectNameSuffix = "");
void ReadInput(const std::string InputFileName);
void FinalizeInitialization(const PhaseField& Phase,
                           const Velocities& Vel,
                           const BoundaryConditions& BC);

FinalizeInitialization must be called after the microstructure is set up to correctly initialize populations and obstacle states based on the phase field.

Time stepping

// Single-component flow
void Solve(PhaseField& Phase, Velocities& Vel, const BoundaryConditions& BC);

// Multi-component flow with composition
void Solve(PhaseField& Phase, const Composition& Cx,
           Velocities& Vel, const BoundaryConditions& BC);

Force computation

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

Diagnostic utilities

size_t CountFluidNodes() const;
size_t CountObstacleNodes() const;
vector<double> CalculateFluidMass() const;
double CalculateLiquidVolume(size_t Fluid_Comp);
double Pressure(int i, int j, int k) const;
dMatrix3x3 PressureTensor(int i, int j, int k) const;
dVector3 FluidVelocity(int i, int j, int k) const;
double FluidDensity(int i, int j, int k) const;

I/O

bool Read(const Settings& locSettings, const BoundaryConditions& BC,
          const int tStep = -1) override;
bool Write(const Settings& locSettings, const int tStep) const override;
void WriteVTK(const Settings& locSettings, const PhaseField& Phase,
              const int tStep, const int precision = 16) const;
void lbWriteVTK(const Settings& locSettings, const PhaseField& Phase,
                const int tStep, const int precision = 16) const;

WriteVTK outputs physical fields (density, velocity, pressure). lbWriteVTK outputs raw lattice populations for debugging.

---

Mass conservation

For moving solid obstacles and phase transformations, explicit mass conservation corrections are applied:

void EnforceMassConservation();        // fluid-only conservation
void EnforceTotalMassConservation(PhaseField& Phase); // solid + fluid
void AccountForAndLimitPhaseTransformation(PhaseField& Phase, int tStep);

These ensure that fluid mass is neither created nor destroyed when solid particles appear, disappear, or change density.

---

Equilibrium distribution

Static helper for computing the Maxwell–Boltzmann equilibrium:

static D3Q27 EquilibriumDistribution(double lbDensity,
    const double weights[3][3][3], dVector3 lbMomentum = {0,0,0});

Used during initialization and in the collision step.

---

Numerical considerations

• The lattice speed of sound is fixed: $c_s^2=1/3$ (hardcoded as lbcs2 = 1.0/3.0).
• The relaxation time $\tau$ must be $>0.5$; values close to $0.5$ lead to instability.
• At high density ratios (two-phase flow), use Do_FixPopulations = true to suppress negative populations near interfaces.
• OptimalParaKuper() provides a recommended Kupershtokh parameter for a given reduced temperature, improving numerical stability of two-phase flows.
• Grid sizes that are powers of two improve memory alignment and, where applicable, FFT-based post-processing performance.

---

Performance

• The D3Q27 stencil has 27 values per node; memory footprint is $2\times 27\times 8\times N_x{N}_y{N}_z$ bytes (two population arrays).
• The collision and streaming steps are $O(N)$ per timestep, $N=N_x{N}_y{N}_z$.
• OpenMP threading is applied across the spatial loop.
• MPI parallelism is supported via domain decomposition; Storage3D handles halo exchange.