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

[李东岳]

欧拉方法下标量传输，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，因为我用这个模拟过气泡流，结果不错。

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

[李东岳]

一个解决办法是将浮力项和重力项进行下面的转换（参考其他求解器）：$\nabla p - \rho \mathbf{g}=\nabla p_{\mathrm{rgh}}+\mathbf{g} \mathbf{h} \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主要用于气固流动，连续相密度远小于固相，空气的密度无法引致这么大的压力差。这不同于气液（连续相密度大于离散相）。

一个解决办法是将浮力项和重力项进行下面的转换（参考其他求解器）：$\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.

[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()
        )
    );
}
//