如果用DPMFoam求解稀相流会怎么样？误差大么
-
@alvin 在 如果用DPMFoam求解稀相流会怎么样？误差大么 中说：
出于什么考虑呢？测试的考虑重力后就不合理了.
我需要深入研究一下。主要起源于体积力重力的处理，从期刊看到的公式来看，重力项的处理不太一样，比如下面这俩个：
我怀疑还是粒子压力梯度的处理问题，http://www.cfd-china.com/topic/1488 在弄清楚之后，可以从代码上进而在结果上解释原因。
目前的猜测，DPMFoam连续性方程里面的
phiForces引起的重力驱动流动。后续我会更新。 -
@东岳 参考您提供的网址 http://dyfluid.com/buoyantPimpleFoam.html 中的部分解析，下面这个方程替换若应用于DPMFoam求解器中不可压缩流场的求解，动量方程中等效于不考虑重力
您提到“OpenFOAM中并没有附加重力的单相流求解器”，显然“buoyantPimpleFoam是OpenFOAM的传热求解器之一，其用于求解瞬态的、由于温度变化导致的密度变化、浮力驱动流动。”，受重力(浮力)驱动。
物理上讲，不管流体可压与否，重力做功会对竖直方向上的流动产生影响。最好还是能够了解到DPMFoam求解器中流场代码植入之所以是现在这个样子，它的设计思想及应用范围，目前测试来看，在求解不可压缩单相流场时，至少它不是通用的。 -
DPMFoam主要用于气固流动，连续相密度远小于固相，空气的密度无法引致这么大的压力差。这不同于气液（连续相密度大于离散相）。
一个解决办法是将浮力项和重力项进行下面的转换（参考其他求解器）：$\nabla p - \rho \mathbf{g}=\nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \nabla \rho$，我植入看看。
更简单的方法是，将DPMFoam代码中进行下面的改动：
surfaceScalarField phicForces ( fvc::flux(rAUc*cloudVolSUSu/rhoc) + rAUcf*(g & mesh.Sf())*1000 );后面的1000是连续相的密度，因为DPMFoam默认求解气固，气体的密度为1. 如果用DPMFoam求解气液，液体的密度假定为1000，乘上去即可。
主要体现在压力求解的正确性，不乘以1000，导致静水压压力大小分布不正确。要进行测试：随便模拟一个1米高的容器，充满水，看内部的压力分布。这种情况压力底部的精确解为101325+9.81*1000*1=111325. 如果不进行代码更正，求解后底部的压力为101325+9.81*1*1=101334.
-
@东岳 谢谢您的耐心回复！
您回复中
提到"......后面的1000是连续相的密度，因为DPMFoam默认求解气固，气体的密度为1. 如果用DPMFoam求解气液，液体的密度假定为1000，乘上去即可。"
既然DPMFoam默认求解气固，气相密度默认为1(适用于我目前求解的算例，就是单纯的空气流场，但结果明显不好)，右端重力项乘以1000中的1000替换为rhoc是否可以？参照 多相流数学模型 http://dyfluid.com/docs/multiphase.html 3.3 连续相模型之欧拉模型（宏观模型）：
对DPMFoam求解器中UcEqn.H进行修正的话,是否方程中各项均应乘或者除以一个rhoc，对连续相的动量方程两端作统一的密度处理，而不仅仅是目前采用的在 (1.0/rhoc)*cloudSU 这一项中由于相间动量耦合而考虑rhoc的影响，(UcEqn.H原代码如下：)// fvVectorMatrix UcEqn ( fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc) - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc) + continuousPhaseTurbulence->divDevRhoReff(Uc) == (1.0/rhoc)*cloudSU ); UcEqn.relax(); volScalarField rAUc(1.0/UcEqn.A()); surfaceScalarField rAUcf("Dp", fvc::interpolate(rAUc)); surfaceScalarField phicForces ( (fvc::interpolate(rAUc*cloudVolSUSu/rhoc) & mesh.Sf()) + rAUcf*(g & mesh.Sf()) ); if (pimple.momentumPredictor()) { solve ( UcEqn == fvc::reconstruct ( phicForces/rAUcf - fvc::snGrad(p)*mesh.magSf() ) ); } // -
@东岳 再次感谢老师的回复。诚如您所说的那样，连续相动量方程中不考虑耦合的颗粒动量传递外，其它各项确实没有考虑连续相密度，默认为1？那么动量耦合项中的密度在数值上也就不起作用了。
您提到"......最近在植入自己的算法，算气固非常完美（连续相密度为1）......"是指“...一个解决办法是将浮力项和重力项进行下面的转换（参考其他求解器）：∇p−ρg=∇prgh+gh∇ρ，我植入看看...”这个解决办法吗？如果是的话，需要了解一下这样转换的逻辑基础，正如interFoam求解器中动量方程右端也是这样的转换，但显然interFoam求解器基于VOF模型，求解两相流采用同一套动量方程（方程中不显示相分数，单独求解相分数方程区分相态）描述两相的状态，两相作为混合物（密度采用均相密度）其在相界面是存在不为0的密度梯度的，∇p−ρg=∇prgh+gh∇ρ 这样的转换有其逻辑基础；DPMFoam求解器的连续相动量方程中左端有相分数（考虑稠密颗粒流，体积分数较大的情况），右端的重力项、压力梯度项中（根据颗粒受不受压力梯度力可考虑在该项中加不加相分数）却没有，和以下公式是否矛盾？（参照 多相流数学模型 http://dyfluid.com/docs/multiphase.html 3.3 连续相模型之欧拉模型（宏观模型）：）
上述方程同除以密度rhoc的话（连续相为不可压缩单相流体，这样处理是可以的），方程右端压力项也要除以rhoc，对比原始DPMFoam动量方程：
//...... surfaceScalarField phicForces ( (fvc::interpolate(rAUc*cloudVolSUSu/rhoc) & mesh.Sf()) + rAUcf*(g & mesh.Sf()) ); if (pimple.momentumPredictor()) { solve ( UcEqn == fvc::reconstruct ( phicForces/rAUcf - fvc::snGrad(p)*mesh.magSf() ) ); }只要让压力项除以rhoc，那么您提到的“最近在植入自己的算法，算气固非常完美（连续相密度为1），算气液就是耦合作用太小（连续相密度为1000），差点劲，感觉是密度方面的原因导致耦合部分的效果太小了。需要重新研究一下。”算气液的话，耦合部分相对于右端的其它项不就提高了吗？我需要测试一下。希望和您保持联系。
-
将pEqn.H改为p_rghEq.H之后，初步看起来是合理的，上图左侧为一个渐进的，下面越来越大的水压，右侧速度为0（1e-10）。这只一种单相附加重力的静态模拟。不过还需要进一步测试。
-
之前的植入，直接将将pEqn.H改为p_rghEq.H，是无密度速度方程+prgh方程。虽然流场预测正确，但这种植入在密度不为1的时候耦合项不正确导致速度过低。
-
下午植入了有密度的速度方程+p方程，在模拟静止的流体的时候压力正确，但有spurious velocity
我测试下带密度的方程+prhg方程。这个应该没问题了吧。目前openfoam这种可压缩动量方程+prgh方程是主流。
模拟静止液体由于数值原因（静水压力和重力平衡）带来的spurious velocity:
-
-
“1.之前的植入，直接将将pEqn.H改为p_rghEq.H，是无密度速度方程+prgh方程。虽然流场预测正确，但这种植入在密度不为1的时候耦合项不正确导致速度过低。”
在求解不可压缩流体时，动量方程等效于不考虑重力（其他单相不可压缩流体求解器如simpleFoam动量方程不考虑重力），无密度的速度方（密度为1）程中的耦合项在颗粒动量传递到连续相中时考虑了密度（密度为连续相的实际密度），这种处理对连续相的种类做了限定，比如说只能是空气。
“下午植入了有密度的速度方程+p方程，在模拟静止的流体的时候压力正确，但有spurious velocity”
您提到在模拟静止的流体的时候压力正确，言外之意，非静止流体还会出现之前速度发散的现象（流速预测过高）？ -
@东岳 在 如果用DPMFoam求解稀相流会怎么样？误差大么 中说：
单相可压DPMFoam求解器，应用领域比较局限了。。。
