如果用DPMFoam求解稀相流会怎么样？误差大么

[李东岳]

基于此(我的st数小于1)，

你的St有多小？

[Alvin]

@李东岳 我的模拟工况目前并没有试验支持，只作为研究用。采用DPMFoam求解器是因为个人原因，稀相到密相流均能采用。不过对于稀相流在多大范围内具有适用性（误差）仍需要调查。您对比“欧拉方法下标量传输，Or DPMFoam求解器”的结果怎么样？

[Alvin]

@李东岳 0.01级别

[李东岳]

欧拉方法下标量传输，Or DPMFoam求解器”的结果怎么样？

这是目前我打算做的，但是目前还没做。所以问问你有没有实验数据支撑。
如果St 0.01，粒子速度不应该和连续相有很大差别，至少书中是这样提及的。并且你的相分数远小于1. 应该可以当做是dusty gas，中文暂时不知道是什么称呼。:

[Alvin]

@李东岳 好的，我再调查一下

[chpjz0391]

@Alvin 你好，请问一下你之前用欧拉方法用的是哪个求解器可以说一下么？就是把稀相流当做是气体考虑的。你用的是twophaseeulerfoam？

[马乔]

稀相可以不考虑固相体积分数，如果考虑颗粒碰撞的话，仍然可以使用CollidingCloud，但是把主程序中的控制方程消去alphac（不能不要，计算曳力有用）

[Alvin]

@chpjz0391 最近利用DPMFoam基于of30 自带算例windaroundbuildings，加了两个烟囱（建筑附近一大一小），单纯计算流场。右侧入口面速度2m/s，大烟囱3m/s，小烟囱2m/s。计算的流场结果和simpleFoam（simpleFoam结果从流场数值上看结果合理）做对比：

计算了15.3s后整体流场结果（最大超过百米）：

下面是我的计算文件以及部分log文件（为上传方便，未加建筑）：
0_1522814173864_DPMFoam.zip
0_1522814350821_log.zip

什么原因呢？DPMFoam求解流场有什么限制吗？

[李东岳]

@alvin

https://coding.net/u/dyfluid/p/Solvers_DYFLUID/git/tree/master/DPMFoam_TestCase?public=true

请在此处下载修改后的算例

[Alvin]

@东岳 谢谢老师，下载了修改后的算例，我再测试一下，稍后给您回复

[Alvin]

@东岳 您给的速度noslip边界条件of3.0版本用不了

[李东岳]

将noSlip改为

type  fixedValue;
value uniform (0 0 0);

[Alvin]

@东岳 谢谢老师，流场确实收敛了。我注意到您做的更改中关掉了重力g（设为（0 0 0）），出于什么考虑呢？测试的考虑重力后就不合理了.

[李东岳]

@alvin 在 如果用DPMFoam求解稀相流会怎么样？误差大么 中说：

出于什么考虑呢？测试的考虑重力后就不合理了.

我需要深入研究一下。主要起源于体积力重力的处理，从期刊看到的公式来看，重力项的处理不太一样，比如下面这俩个：

我怀疑还是粒子压力梯度的处理问题，http://www.cfd-china.com/topic/1488 在弄清楚之后，可以从代码上进而在结果上解释原因。

目前的猜测，DPMFoam连续性方程里面的phiForces引起的重力驱动流动。后续我会更新。

[Alvin]

谢谢您耐心的讲解，受教了，可是如果关掉重力，DPmFoam求解器是不是就不能考虑颗粒受重力的情况了，感觉这个求解器的处理不是很理解

[李东岳]

问题是phiForces导致，比如重力向下，那么phiForces将为负值，其产生一个向下的速度。目前我自己的求解器也遇到了这个问题。还需要进一步研究。还不能说是DPMFoam的bug，因为我用这个模拟过气泡流，结果不错。

我要看一下密度为常数的情况下，速度和压力的平衡处理。

[李东岳]

一个解决办法是将浮力项和重力项进行下面的转换（参考其他求解器）：$\mathrm{\nabla }p-\rho \mathbf{g}=\mathrm{\nabla }{p}_{\mathrm{r}\mathrm{g}\mathrm{h}}+\mathbf{g}\mathbf{h}\mathrm{\nabla }\rho$，我植入看看。

[Alvin]

真是太谢谢您了，不过排除了案例本身设置的问题外，从物理直观上不好理解同样求解不可压缩流场，simpleFoam不考虑重力，DPMFoam流体动量方程考虑重力就不合理。

[李东岳]

将pEqn.H改为p_rghEq.H之后，初步看起来是合理的，上图左侧为一个渐进的，下面越来越大的水压，右侧速度为0（1e-10）。这只一种单相附加重力的静态模拟。不过还需要进一步测试。

http://dyfluid.com/buoyantPimpleFoam.html 我在这里添加了点内容，跟这个有关

[Alvin]

@东岳 参考您提供的网址 http://dyfluid.com/buoyantPimpleFoam.html 中的部分解析，下面这个方程替换若应用于DPMFoam求解器中不可压缩流场的求解，动量方程中等效于不考虑重力

您提到“OpenFOAM中并没有附加重力的单相流求解器”，显然“buoyantPimpleFoam是OpenFOAM的传热求解器之一，其用于求解瞬态的、由于温度变化导致的密度变化、浮力驱动流动。”，受重力(浮力)驱动。
物理上讲，不管流体可压与否，重力做功会对竖直方向上的流动产生影响。最好还是能够了解到DPMFoam求解器中流场代码植入之所以是现在这个样子，它的设计思想及应用范围，目前测试来看，在求解不可压缩单相流场时，至少它不是通用的。

[李东岳]

谢谢。

OpenFOAM中并没有附加重力的单相流求解器

已更改为

OpenFOAM中并没有附加重力的恒密度单相流求解器。例如在恒密度的icoFoam和simpleFoam中，是压力差引起的流体流动。

最好还是能够了解到DPMFoam求解器中流场代码植入之所以是现在这个样子，它的设计思想及应用范围，目前测试来看，在求解不可压缩单相流场时，至少它不是通用的。

DPMFoam主要用于气固流动，连续相密度远小于固相，空气的密度无法引致这么大的压力差。这不同于气液（连续相密度大于离散相）。

[李东岳]

DPMFoam主要用于气固流动，连续相密度远小于固相，空气的密度无法引致这么大的压力差。这不同于气液（连续相密度大于离散相）。

一个解决办法是将浮力项和重力项进行下面的转换（参考其他求解器）：$\mathrm{\nabla }p-\rho \mathbf{g}=\mathrm{\nabla }{p}_{\mathrm{r}\mathrm{g}\mathrm{h}}+\mathbf{g}\mathbf{h}\mathrm{\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.

[Alvin]

@东岳 谢谢您的耐心回复！

您回复中

提到"......后面的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()
        )
    );
}
//

[李东岳]

@alvin
非常感谢，确实如此，我再重新植入一下。

最近在植入自己的算法，算气固非常完美（连续相密度为1），算气液就是耦合作用太小（连续相密度为1000），差点劲，感觉是密度方面的原因导致耦合部分的效果太小了。需要重新研究一下。

[Alvin]

@东岳 再次感谢老师的回复。诚如您所说的那样，连续相动量方程中不考虑耦合的颗粒动量传递外，其它各项确实没有考虑连续相密度，默认为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），差点劲，感觉是密度方面的原因导致耦合部分的效果太小了。需要重新研究一下。”算气液的话，耦合部分相对于右端的其它项不就提高了吗？我需要测试一下。希望和您保持联系。

[李东岳]

回帖很复杂。。。

我这几天做一个带密度的可压缩的DPMFoam试试。

[李东岳]

将pEqn.H改为p_rghEq.H之后，初步看起来是合理的，上图左侧为一个渐进的，下面越来越大的水压，右侧速度为0（1e-10）。这只一种单相附加重力的静态模拟。不过还需要进一步测试。

1. 之前的植入，直接将将pEqn.H改为p_rghEq.H，是无密度速度方程+prgh方程。虽然流场预测正确，但这种植入在密度不为1的时候耦合项不正确导致速度过低。

2. 下午植入了有密度的速度方程+p方程，在模拟静止的流体的时候压力正确，但有spurious velocity

我测试下带密度的方程+prhg方程。这个应该没问题了吧。目前openfoam这种可压缩动量方程+prgh方程是主流。

模拟静止液体由于数值原因（静水压力和重力平衡）带来的spurious velocity:

[Alvin]

“1.之前的植入，直接将将pEqn.H改为p_rghEq.H，是无密度速度方程+prgh方程。虽然流场预测正确，但这种植入在密度不为1的时候耦合项不正确导致速度过低。”
在求解不可压缩流体时，动量方程等效于不考虑重力（其他单相不可压缩流体求解器如simpleFoam动量方程不考虑重力），无密度的速度方（密度为1）程中的耦合项在颗粒动量传递到连续相中时考虑了密度（密度为连续相的实际密度），这种处理对连续相的种类做了限定，比如说只能是空气。
“下午植入了有密度的速度方程+p方程，在模拟静止的流体的时候压力正确，但有spurious velocity”
您提到在模拟静止的流体的时候压力正确，言外之意，非静止流体还会出现之前速度发散的现象（流速预测过高）？

[Alvin]

@东岳 在 如果用DPMFoam求解稀相流会怎么样？误差大么 中说：
单相可压DPMFoam求解器，应用领域比较局限了。。。

[Alvin]

“...我测试下带密度的方程+prhg方程。这个应该没问题了吧。目前openfoam这种可压缩动量方程+prgh方程是主流”
可压缩动量方程植入，至少还要添加能量方程，这样处理的话更加复杂了......是否单纯将重力项在有密度的连续性方程中去掉就好

[linhan.ge]

@alvin
DPMfoam的压力场给定的时候不是已经除过密度了吗？从0文件夹中p文件的量纲可以看出来。
压力相的相分数是由颗粒相传递过来的压力力梯度导致的。

[Alvin]

@linhan-ge 谢谢您的解答。从量纲上压力确实除以了密度，但是压力的数值处理上，动量方程中除以的密度默认为1（适用于空气），这是最重要的。“压力相的相分数是由颗粒相传递过来的压力力梯度导致的”，若是稀相流，压力梯度力不考虑，压力项中也就不存在相分数了。回归正题，DPMFoam求解的稀相流，误差有多大（我自己用DPMFoam单纯计算流场的算例考虑重力后是发散的）？DPMFoam动量方程中的重力项该如何处理才好（舍掉可好）？

[李东岳]

DPMFoam动量方程中的重力项该如何处理才好（舍掉可好）？

如果连续相给离散相的浮力很重要（比如气液），不可以舍去。如果连续相给离散相的浮力不重要（比如气固），可以舍去。

[Alvin]

@东岳 “如果连续相给离散相的浮力很重要（比如气液），不可以舍去。如果连续相给离散相的浮力不重要（比如气固），可以舍去。”这里不太理解， 连续相对离散相的力的作用（如浮力）不都是通过颗粒拉格朗日方程体现的么？DPMFoam动量方程中的重力项单纯是流体的重力？

[dzw05]

@alvin 没有太理解，“但是压力的数值处理上，动量方程中除以的密度默认为1（适用于空气），这是最重要的”这句话。
从可压缩动量方程出发，假设密度是常数，那可以将密度提出来，就是现在采用的不可压缩动量方程形式啊，为什么会默认密度为1呢？

[dzw05]

@东岳 目前的DPM求解器中，浮力项并不单纯由连续相的重力表示的，还包括连续相的加速度项，可以参考“INFLUENCE OF DRAG LAWS ON SEGREGATION AND BUBBLING BEHAVIOR IN GAS–FLUIDIZED BEDS”这篇博士论文的2.5.1节的内容。
具体来说在Openfoam中，particleForces中的pressureGradient这项才是浮力项。

[dzw05]

@alvin DPMFoam动量方程中的重力项有一部分是流体的重力，另一部分和pressureGradient这一项组合起来作为浮力。具体请参考“INFLUENCE OF DRAG LAWS ON SEGREGATION AND BUBBLING BEHAVIOR IN GAS–FLUIDIZED BEDS”这篇博士论文的2.5.1节。

[李东岳]

particleForces中的pressureGradient这项才是浮力项。

gravity中植入了浮力项。pressureGradient不是浮力项。pressureGradient是粒子上下表面引致的压力差项，粒子若非常小，则可忽略。但若重力gravity若被忽略，则完全没有浮力影响。

可进行类似浮力测试：

particleForces
{
    //constantDrag
    //{
    //	constantValue	0.1;
    //}
 
    gravity; //浮力On or Off
 
     	//pressureGradient
     	//{
     	//    U	U.water;
     	//}
 
    }

DPMFoam动量方程中的重力项单纯是流体的重力？

考虑连续相的密度大于离散相的密度（气液）。重力的作用：连续相植入重力项，导致静水压力的产生。离散相植入重力，同样导致离散相的下降。重力的作用对连续相和离散相都是一个方向的。但重力的植入同时会导致浮力的产生，浮力方向和重力相反，且浮力远大于重力，因此离散相向上移动。

如果连续相给离散相的浮力很重要（比如气液），不可以舍去。如果连续相给离散相的浮力不重要（比如气固），可以舍去

因为浮力和重力相关联。因此如果浮力很重要（气液中的气泡要上浮），连续相的重力的植入必须。如果浮力不重要（气固中的颗粒不会上浮），连续相的重力植入则无关紧要。

[dzw05]

@东岳 没错，在离散相的gravity中，表达式是massg_(1.0 - p.rhoc()/p.rho())，即考虑了连续相和离散相的密度比，这个是浮力的主要来源。
目前dpmFoam中，如果将重力的值设为0，那么就没有考虑浮力影响；如果单纯的将连续相动量方程中的重力g删除，而不是设为0，那么其实并没有影响粒子浮力的计算，因为粒子运动方程中根本没有出现连续相的压力。

[李东岳]

@dzw05 在 如果用DPMFoam求解稀相流会怎么样？误差大么 中说：

粒子运动方程中根本没有出现连续相的压力

这个是问题的根源。因为存在俩套方程：

1. 第一种是你说的，不存在连续相的压力，DEM可以这样模化：

   $$\begin{array}{}\text{(1)}& m_{\mathrm{d}\mathrm{p}\mathrm{m}}\frac{\mathrm{\partial }{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}}{\mathrm{\partial }t}={\mathbf{F}}_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}+\frac{\mathbf{g}({\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}-{\rho }_{\mathrm{c}})}{{\rho }_{\mathrm{c}}}+\mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s}\end{array}$$

   在这个方程中，如果不考虑F，那么就可以模化重力和浮力的效果。这就是DPMFoam中采用的方法。

2. 另一套方程是

   $$\begin{array}{}\text{(2)}& m_{\mathrm{d}\mathrm{p}\mathrm{m}}\frac{\mathrm{\partial }{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}}{\mathrm{\partial }t}={\mathbf{F}}_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}+m_{\mathrm{d}\mathrm{p}\mathrm{m}}\mathbf{g}+\frac{m_{\mathrm{d}\mathrm{p}\mathrm{m}}}{{\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}}\mathrm{\nabla }{p}_{\mathrm{c}}+\mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s}\end{array}$$

   如果采用这一套方程，其中的压力需要为静水压力。

目前我在研究第二种方法。Fluent以及OpenFOAM中选用的第一种方法。参考Fluent理论指南。

[Alvin]

@dzw05 参照 多相流数学模型 http://dyfluid.com/docs/multiphase.html 3.3 连续相模型之欧拉模型（宏观模型）：

“从可压缩动量方程出发，假设密度是常数，那可以将密度提出来，就是现在采用的不可压缩动量方程形式啊，为什么会默认密度为1呢？”
假设密度是常数，那么即使将密度提出来，DPMFoam求解器中的动量方程中压力项还少除以了一个密度（除非密度是1）

“具体来说在Openfoam中，particleForces中的pressureGradient这项才是浮力项。”
颗粒所受的浮力在颗粒所受重力项中已经体现（）： pressureGradient就是单纯的压力梯度力，但是大小和浮力相反（可参考颗粒在流体中的受力分析 https://wenku.baidu.com/view/171125dcd15abe23482f4d93.html、董长银-牛顿流体中的固体颗粒运动模型分析及应用https://wenku.baidu.com/view/8fafedfef705cc175527097a.html）

“DPMFoam动量方程中的重力项有一部分是流体的重力，另一部分和pressureGradient这一项组合起来作为浮力。” 您是说连续相的动量方程中的重力项要被拆分吗？据我调查，颗粒所受重力中已经考虑了流体的浮力影响，那么，颗粒对流体应该有相反的作用力（体积力的源项），那么连续相动量方程里的重力项似乎可以解释了，是作为颗粒的耦合项加入的。

@东岳 “pressureGradient是粒子上下表面引致的压力差项，粒子若非常小，则可忽略。但若重力gravity若被忽略，则完全没有浮力影响。” 颗粒浮力的影响通过连续相的密度以及重力加速度体现出来，关于https://wenku.baidu.com/view/8fafedfef705cc175527097a.html文献中提到的考虑了压力梯度力后，浮力项不起作用了，因为两者大小相等，方向相反，在openfoam中还没有验证，应该根据具体案例和压力分布相关

[Alvin]

@dzw05
“没错，在离散相的gravity中，表达式是massg_(1.0 - p.rhoc()/p.rho())，即考虑了连续相和离散相的密度比，这个是浮力的%(#c43d3d)[主要来源]。” 浮力还有哪些其它来源呢？
“如果单纯的将连续相动量方程中的重力g删除，而不是设为0，那么其实并没有影响粒子浮力的计算，因为粒子运动方程中根本没有出现连续相的压力。” 粒子运动方程中出现的连续相的压力，不是可以通过压力梯度力来体现出来了吗？这样一来，连续相动量方程中的耦合项或者压力项需做相应调整了。

@东岳 您提出存在下面的两套方程，是等效的吗？

[李东岳]

这俩套方程都是文献中看到的，可以参考本帖24楼。我自己的求解器中采用这种形式

$$\begin{array}{}\text{(3)}& m_{\mathrm{d}\mathrm{p}\mathrm{m}}\frac{\mathrm{\partial }{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}}{\mathrm{\partial }t}={\mathbf{F}}_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}+m_{\mathrm{d}\mathrm{p}\mathrm{m}}\mathbf{g}+\frac{m_{\mathrm{d}\mathrm{p}\mathrm{m}}}{{\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}}\mathrm{\nabla }{p}_{\mathrm{c}}+\mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{s}\end{array}$$

俩种形式可以互相转换。依据$p=p_0+{\rho }_{\mathrm{c}}\mathbf{g}\mathbf{h}$即可相互推导

[dzw05]

@alvin @东岳 贴出来的两套方程才是问题的根源所在。目前OpenFOAM中的压力梯度项并没有出现连续相的压力，而是所谓的rhoc*Du/Dt。

@alvin 不可压缩动量方程中的p，其实是p/rho，这个和可压缩方程中是不同的，simpleFoam中也是这样处理的啊。难道咱们讨论的方程不是同一个吗？至少采用DPMFoam计算液-固流动时，连续相的压力是没有问题的。

[李东岳]

$$\begin{array}{}\text{(4)}& m_{\mathrm{d}\mathrm{p}\mathrm{m}}\frac{\mathrm{d}{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}}{\mathrm{d}t}=-\frac{m_{\mathrm{d}\mathrm{p}\mathrm{m}}}{{\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}}\mathrm{\nabla }p+m_{\mathrm{d}\mathrm{p}\mathrm{m}}\mathbf{g}+\frac{m_{\mathrm{d}\mathrm{p}\mathrm{m}}}{{\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}}\frac{3}{4}\frac{C_{\mathrm{D}}{\rho }_{\mathrm{c}}}{d_{\mathrm{d}}}|{\mathbf{U}}_{\mathrm{c}}-{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}|({\mathbf{U}}_{\mathrm{c}}-{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}})+\mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{s}.\end{array}$$
伯努利：
$$\begin{array}{}\text{(5)}& \mathrm{\nabla }p={\rho }_{\mathrm{c}}\mathbf{g}\end{array}$$
代入：
$$\begin{array}{}\text{(6)}& m_{\mathrm{d}\mathrm{p}\mathrm{m}}\frac{\mathrm{d}{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}}{\mathrm{d}t}=m_{\mathrm{d}\mathrm{p}\mathrm{m}}\mathbf{g}(1-\frac{{\rho }_{\mathrm{c}}}{{\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}})+\frac{m_{\mathrm{d}\mathrm{p}\mathrm{m}}}{{\rho }_{\mathrm{d}\mathrm{p}\mathrm{m}}}\frac{3}{4}\frac{C_{\mathrm{D}}{\rho }_{\mathrm{c}}}{d_{\mathrm{d}}}|{\mathbf{U}}_{\mathrm{c}}-{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}}|({\mathbf{U}}_{\mathrm{c}}-{\mathbf{U}}_{\mathrm{d}\mathrm{p}\mathrm{m}})+\mathrm{O}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{F}\mathrm{o}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{s}.\end{array}$$

[Alvin]

@dzw05 抱歉，理解有限，“ 目前OpenFOAM中的压力梯度项并没有出现连续相的压力，而是所谓的rhoc*Du/Dt。” 这句话不是很理解.
“不可压缩动量方程中的p，其实是p/rho，这个和可压缩方程中是不同的，simpleFoam中也是这样处理的啊。”从量纲上确实是这样，查看过simpleFoam的方程，代码中没发现有p/rho这样的处理，您能从源码角度讲一下它的实现过程吗（源码位置？），希望不吝赐教？

@东岳 再次谢谢老师的指导

[dzw05]

@alvin （1）rhoc*Du/Dt，这个是在src\lagrangian\intermediate\submodels\Kinematic \ParticleForces \PressureGradient\PressureGradientForce.T.C代码中出现的，应该对应于东岳写的第一套方程。
（2）至于压力这个，源码确实没有出现，但推导的时候就是这样出来的。simpleFoam计算输出的p也不是真正的压力，而是p/rho，你要看真实的压力的话，应该乘以连续相的密度。

[Alvin]

@东岳 这样说来，利用现有的DPMFoam（3.0.1版本）计算气液两相流，若气泡同时考虑压力梯度力和重力（massg_(1.0 - p.rhoc()/p.rho())），可能会导致气泡浮力失效？

@dzw05 谢谢您，我再看一下。“至于压力这个，源码确实没有出现，但推导的时候就是这样出来的。simpleFoam计算输出的p也不是真正的压力，而是p/rho，你要看真实的压力的话，应该乘以连续相的密度。”
谢谢您的解答。

@东岳 @dzw05 离散相方程中的压力梯度力确实是rhocDu / Dt这种形式，和浮力（rhocg）是两个概念

[李东岳]

@Alvin 不知道你的问题解决没有。我在我自己写的求解器中，通过

1. 将p方程换为p_rgh方程
2. 将g去掉

解决了。因为轻颗粒的浮力在方程$\text{(1)}$中包含了，重颗粒也如此。就像 @dzw05 说的，连续相可以认为是单向不可压缩流，类似simpleFoam。 @Alvin 你试试将DPMFoam中的phig去掉，运算看看？

[Alvin]

@东岳 问题解决了，如您所说，去掉动量方程中的重力加速度项，案例结果就合理了。我在考虑，动力方程中的重力项的物理含义是什么，颗粒受到周围流体对它施加的浮力，反之，根据牛顿第三定律，流体也应该受到对应颗粒浮力的反作用力，方向相反，这个重力项是否就是这个反作用力，尚未查到相关文献，只是一种猜测。