particleforce类内ErgunWenYuDrag的表达式感觉有错误

[李东岳]

你看看这个文章？Drag models for Simulation Gas-Solid Flow in the Bubbling Fluidized Bed of FCC Particles
看起来跟你的方程差不多

[上级]

李老师您这篇文章的公式就是我上面式（4）的公式，但是我好像打错了 。总体比较来看是有细微区别。

[刘雄国]

现在OpenFOAM7 出来了 据我的同事说改变了很多 你可以看一下

[李东岳]

@上级 我最近在弄拉格朗日欧拉耦合，发现了点关于相分数的东西，你的信息量好大我没太看全，不是好像也是关于想分数的？下图会不会给你一点启发

[上级]

李老师@东岳 请问这篇文献出处是？我想研究一下

[上级]

@刘雄国 嗯，我去看一下OF7的代码，谢谢您

[李东岳]

@上级 我写的manuscript还没提交，http://www.cfd-china.com/topic/1488/欧拉及拉格朗日下的动量方程压力项竟然不一样/7 这个我包含了一些推导

[Ant]

@上级 您好，您现在又有新的发现吗，如果Re的定义中（公式6）再多一个alphac,这样是不是就和文献中一样了？公式（8）中mass/density那部分应该是1-alphc（你上面说的），如果Re在MPPIC/DPMfoam的定义是alphac乘以特征项那就应该是一样了. 但是我不知道怎么找Re的定义，如果您这边对代码比较熟麻烦核实下看看？

[Ant]

@东岳 李老师，您好，目前在看一些文献，关于相分数这些搞不懂。看了您下面链接的相关讨论，也看了关于Z. Y. ZhOU J. Fluid Mech. (2010), vol. 661, pp.482-510. 这篇文章里关于相分数的讨论。那篇文章中有关于压力梯度项和应力散度项乘以alpha和不乘两种形式，会导致动力传递项F包含的力不同。您这篇文章里重力也没有乘以alpha，这个是有特殊含义吗，后面那个动力传递项是不是也与传统的几中写法不同？看了您上面关于这两个的解释，也不是很懂？

[李东岳]

关于压力梯度项和应力散度项乘以alpha和不乘两种形式，会导致动力传递项F包含的力不同。

http://www.cfd-china.com/topic/1488/欧拉及拉格朗日下的动量方程压力项竟然不一样/7 应该在这个里面讨论了。乘以相分数之后不需要压力梯度力，不乘相分数需要考虑压力梯度力

[上级]

@东岳 李老师，关于您这张图中的公式（16）我有几个疑问，如果说得不对望见谅。
1.对于连续相动量方程的拉力梯度项前是否需要乘以一个相分数主要决定于$\mathbf{F}$中的压力梯度力项是否提取出来了，如果提取出来了和前面的压力梯度项组合就会出现乘以相分数的情况，如果不提取出来，把颗粒的压力梯度力放在$\mathbf{F}$中，前面就不用再乘以相分数。不知道我这么解释是否正确？这个解释我是参考的文章Zhou et al. 2010
2.对于您式子（16）中最后一项为何$\mathbf{F}$的前面要乘以一个颗粒相分数，并且为何是除以颗粒的体积，按道理不应该是除以流场网格体积吗？
希望能得到李老师的解答。

[李东岳]

1 是的

2 这是每单位体积的粒子所受的力，因此除了颗粒体积。另外，是每个网格内粒子的体平均，因此乘以相分数

[上级]

@东岳 李老师，您这个式子（16）的$\mathbf{g}$前面是不是还少一个$\alpha$?

[李东岳]

@上级 是的，。重力和压力梯度相抗衡，如果压力梯度没有乘以$\alpha$，g也不能乘。我觉得你的Wenyu模型相分数问题也应该类似

[upc_ngh]

@上级
您好，请问您上面的公式 F=Sp(U−Ud)+Su用来计算颗粒受力的，这里面的U和Ud分别表示的是什么东西？和KinematicParcel中的Uc_以及U_是对应关系吗

const parcelType& p = static_cast<const parcelType&>(*this);
 const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
    const forceSuSp Feff = Fcp + Fncp;
    const scalar massEff = forces.massEff(p, mass);
    // Update velocity - treat as 3-D
    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
    const scalar bp = Feff.Sp()/massEff;
 
    Spu = dt*Feff.Sp();
 
    IntegrationScheme<vector>::integrationResult Ures =
        td.cloud().UIntegrator().integrate(U_, dt, abp, bp);
 
    vector Unew = Ures.value();
 
    // note: Feff.Sp() and Fc.Sp() must be the same
    dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());
 
    // Apply correction to velocity and dUTrans for reduced-D cases
    const polyMesh& mesh = td.cloud().pMesh();
    meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
    meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);
 
    return Unew;
}

如果我想要计算颗粒所受的ErgunWenYuDrag这个力的时候是不是就等于

(mass/p.rho() *(150.0*(1.0 - alphac)/alphac +1.75*Re)*muc/(alphac*sqr(p.d()))*（U-Ud）//这样就OK了吗？

[上级]

@upc_ngh
1.是的，$\mathbf{U}$代表的是KinematicParcel中的${\mathbf{U}}_{\mathbf{c}}$，也就是连续相的在颗粒所处位置的速度插值，具体的插值方式OF里面有cell，cellpoint那些，${\mathbf{U}}_{\mathrm{d}}$表示的就是颗粒的速度。
2.是的，系数乘以速度差。