DPMFoam速度方程代码与实际物理方程如何对应？

mechanicsdog

最近看DPMFoam的求解器源代码，速度方程的代码一直没看懂，速度方程文件中的代码如下

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevTau(Uc)
 ==
    (1.0/rhoc)*cloudSU
);

UcEqn.relax();

fvOptions.constrain(UcEqn);

volScalarField rAUc(1.0/UcEqn.A());
surfaceScalarField rAUcf("Dp", fvc::interpolate(rAUc));

surfaceScalarField phicForces
(
   fvc::flux(rAUc*cloudVolSUSu/rhoc) + rAUcf*(g & mesh.Sf())
);

if (pimple.momentumPredictor())
{
    solve
    (
        UcEqn
     ==
        fvc::reconstruct
        (
            phicForces/rAUcf - fvc::snGrad(p)*mesh.magSf()
        )
    );

    fvOptions.correct(Uc);
}

而实际的物理方程为
替代文字

问题总结如下

1、代码

fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)

表示物理方程的左边两项好像比较清晰，但是代码

- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)

表示的意义有点难以理解。其中

fvc::ddt(alphac) + fvc::div(alphaPhic)

貌似恰好是连续性方程？

2、应力项代码

+ continuousPhaseTurbulence->divDevTau(Uc)

为何与孔隙率alphac无关，好像和实际物理方程无法对应（左边第三项）？

求解！谢谢。

李东岳

你的第一个理解是对的，就是连续性方程

第二个，则是在湍流模型里面已经把alpha找到了，类似rho的代码也一样

mechanicsdog

谢谢@东岳老师，但我依旧有两个后续问题想请教一下：
1、既然

- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)

表示连续性方程和速度的乘积。
从物理意义上来说，连续性方程恒为0，那么将连续性方程放到速度方程里求解的意义是什么？
是为了更好的数值稳定性吗？还是说为了修正速度值？
2、我在湍流模型里确实找到了divDevTau的代码，如下

template<class BasicMomentumTransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicMomentumTransportModel>::divDevTau
(
    const volScalarField& rho,
    volVectorField& U
) const
{
    return
    (
      - fvc::div((this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
    );
}

代码里面确实包含了密度rho和孔隙率alpha。由于我发现由于是不可压缩流体，速度方程代码两端都是除以密度rho表示的，然而应力求解出来却没有除以rho，这是不是与物理方程相违背
替代文字
是不是速度方程应该改为

    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + (1.0/rhoc)*continuousPhaseTurbulence->divDevTau(Uc)//这一项除以密度rho
 ==
    (1.0/rhoc)*cloudSU

比较正确？

mechanicsdog

奥第二个问题刚刚查代码搞清楚了，对于不可压缩流动的应力计算，密度都设为1了，所以速度方程代码的量纲没有问题。

李东岳

@mechanicsdog 在 DPMFoam速度方程代码与实际物理方程如何对应？中说：

从物理意义上来说，连续性方程恒为0，那么将连续性方程放到速度方程里求解的意义是什么？
是为了更好的数值稳定性吗？还是说为了修正速度值？

更好的有界性。在顺泰迭代的时候，时间步长内是不一定收敛的。有的时候可能会越界。主要起到这个作用。

mechanicsdog

十分感谢！@东岳