@李东岳
@coolhhh
感谢两位老师，我提出这个问题的主要原因是在现有的CFD-DEM流固耦合计算中，一种基于FD/IBM的细网格方法是需要对速度场U进行改变，并针对连续性方程进行修正的。而这个连续性修正的过程就需要再次显式构建连续性条件：div(U)=0。当然这是理论上的条件。
$\nabla \cdot (\nabla \phi) = \nabla \cdot U$
而下面代码中的连续性修正因子的phiIB理论上只在颗粒存在的条件下才有意义，同样才有值。但是如果div(U)在颗粒不存在得情况就不为0了，那么这样的修正就是错误的。虽然当前这个方面的主要算法就是下面展示的这样。
目前看来在PISO外部，U和phi是分离的，但是CFD-DEM的修正是针对U的，这样操作应该是有问题的，误差大小有待考察。而我原本的目标是希望能够在PISO外，得到满足连续性的U或者phi，不过目前看来是很难实现的。
上述过程具体代码如下：
void Foam::cfdemCloud::calcVelocityCorrection
(
volScalarField& p,
volVectorField& U,
volScalarField& phiIB,
volScalarField& voidfraction//颗粒体积分数（颗粒内部为1
)
{
void Foam::cfdemCloudIB::calcVelocityCorrection
(
volScalarField& p,
volVectorField& U,
volScalarField& phiIB,
volScalarField& voidfraction
)
{
setParticleVelocity(U);//改变颗粒所在区域流体的速度
// 修改速度会使其不满足连续性，因此使用一个phiIB进行修正
fvScalarMatrix phiIBEqn
(
fvm::laplacian(phiIB) == fvc::div(U) + fvc::ddt(voidfraction)
);
if(phiIB.needReference())
{
phiIBEqn.setReference(pRefCell_, pRefValue_);
}
phiIBEqn.solve();
U=U-fvc::grad(phiIB);
U.correctBoundaryConditions();
// correct the pressure as well
p=p+phiIB/U.mesh().time().deltaT();
p.correctBoundaryConditions();
}
