I found a bug existing in the palabos code.
In “/src/basicDynamics/externalForceDynamics.hh”,
the force in ExternalForceDynamics::computeVelovity means force per unit density (or force per unit mass).
However, in “/src/latticeBoltzmann/externalForceTemplates.h”,
the force term in externalForceTemplatesImpl::addGuoForce means force itself (or force per unit volume).
Thus, when we use GuoExternalForceBGKdynamics to calculate sigle component multiphase coexistence states,
the density will diverge.
After the following modification:
in ExternalForceDynamics::computeVelovity,
change
u[iD] = j[iD]*invRho + force[iD]/(T)2;
to
u[iD] = j[iD]invRho + force[iD]/(T)2invRho;
The problem is solved.
Yujie