MPPICInterFoam 的Ueqn中到底該不該再乘以alphac

zhe

大家好,之前在其他的帖子裏也發過類似的問題,但最後也沒有解決掉 (沒有一個標準答案).在MPPICInterFoam的Ueqn中,密度rho應該已經經過alpha1$\cdot$rho1+alpha2$\cdot$rho2得到的整體流體密度.那麼是不是在書寫公式時,只需要一個rho就可以代表alphc$\cdot$rho? 從code裏看,在ddt(alphaRho, U)和divDevRhoReff(U)裏都考慮了alphac$\cdot$rho.而對流項,只是用div(rhoPhi,U)表達了,並沒有再用alphac$\cdot$rho. 我個人比較傾向直接用rho表達,alphac$\cdot$rho就兩次引進了體積分數. 希望能得到個準確的答案. 謝謝大家的解惑.

qingdong_wang

我也很纠结, 我问过好几个人, 他们也不清楚. 不过我的算例中, 乘和不乘对结果影响不大. 我的是dilute flow

李东岳

fvVectorMatrix UEqn
(
    fvm::ddt(alphacRho, U)
  + MRF.DDt(alphacRho, U)
  - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
  + fvm::div(rhoPhi, U)
  + turbulence->divDevRhoReff(U)
  ==
    fvOptions(rho, U)
  + cloudSU
);

之前好像看到过类似的问题，为何不把上面的代码变成

fvVectorMatrix UEqn
(
    fvm::ddt(rho, U)
  + MRF.DDt(rho, U)
  - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
  + fvm::div(rhoPhi, U)
  + turbulence->divDevRhoReff(U)
  ==
    fvOptions(rho, U)
  + cloudSU
);

你说的就是上面这种形式吧？然后计算试试？

目前没有一行一行的推导这个求解器方程，大略看一眼好像方程植入不一致，不过需要更深层次的debug

zhe

@东岳改成您寫的這個形式,相當於全部都只是用rho,不再考慮alphac.但是在應力項中也得把alphac去掉.
您這個就相當於跟interFoam一樣了,但是MPPICInterFoam用的是可壓縮的模型,考慮到了固體粒子的相分數和密度. 那麼是不是代表,可以把流體相(氣體+液體)整體的密度rho通過alpha1rho1+alpha2rho2計算出來,然後再通過alphac+theta=1來跟固相耦合.那麼就需要再考慮每項都需要alphac. 那麼就是說根據現在的公式只要把公式裏的對流相乘以個alphac. 我會把兩種情況都試試,後期跟新上來.如果有覺得我這個觀點有問題的,歡迎隨時更正,謝謝.

yfclark

我觉得应该，但是传统的欧拉朗日方法基本认为粒子不占体积，所以你有没有alphac基本没有影响

zhe

@yfclark 謝謝回覆.那麼高濃度粒子的情況呢?像MPPIC可以允許粒子佔體積65%,這樣的情況是不是會有很大的影響了?

李东岳

MPPICInterFoam虽然考虑密度，但是是不可压缩。另外，这里面的$\rho$是$\alpha_1\rho_1+\alpha_2\rho_2$，所以原代码应该是$\nabla\cdot((\alpha_1\rho_1+\alpha_2\rho_2)\bfU)$，如果你要加$\alpha$进去，是要换成$\nabla\cdot(\alpha_1\rho_1\bfU)$么？还是$\nabla\cdot(\alpha_1(\alpha_1\rho_1+\alpha_2\rho_2)\bfU)$。

我想推一下这个。不过最近时间太少了如果你们能分享下最终的方程就太好啦

zhe

考慮固體粒子也佔有一部分體積,並且將氣體和液體聯合求解,我推了一下氣液混合的密度公式,並不是rho=rho1*alpha1+rho2*alpha2.而是alphac*rho=alpha1*rho1+alpha2*rho2. 所以在原來creatfile.H裏的定義應該是有問題的,請查證.那也就是說,凡是涉及到通過密度求質量的公式,在MPPICInterFoam裏都應該用alphac*rho. 那麼在對流項,我覺得應該改成div(alphacRho, U). 有錯的話,希望指出,謝謝.

zhe

@zhe 不好意思,打錯了.最後的對流項應該是div(alphaRhoPhic, U).

李东岳

@zhe 目前不好评论。需要深入去看代码，实在没时间了，我可以把手头的manuscript写完之后，再处理这个求解器。手头有好多求解器要处理 buoyantPimpleFoam和Osher格式也压了好久了.. 不过我会持续跟进你的讨论！也感谢分享！

zhe

@东岳好的,多謝回覆.我會接着再試試重新定義rho,等您有時間了,希望能多指導.

bestucan

说两个猜想。
如果是alphc(fake)*rho = alpha1*rho1 + alpha2*rho2。
那得是alpha1、alpha2都是未归一化的值。
（如果alpha1 + alpha2 = 1，那么rho = alpha1*rho1 + alpha2*rho2，未找到alpha1和alpha2的关系。）
这种情况下rho = (alpha1*rho1 + alpha2*rho2)/alphac(fake)即密度=质量/体积

实际上代码里的rho(code) = alpha1*rho1 + alpha2*rho2，即上面的质量。
即rho(code)是体积为alphac(fake)的混合物的质量
即rho(code)是单位alphac(fake)的质量

下面开始猜：
alphac(code)是归一化过的，在代码里有显示，是1-xxx。
那么猜测：与alpha1和alpha2类似的数值alphac(real)隐含在alphac(code)中。
即alphac(code)=alphac(real)/alphac(fake)

那么rho(code) * alphac(code)代入以上定义就等于：(alpha1*rho1 + alpha2*rho2)/alphac(fake)*alphac(real)
即rho*alphac(real)
(rho(code) * alphac(code) -> rho*alphac(real))

上面是猜rho可能的定义。下面是猜有的地方乘alphac，有的地方不乘。

上面说了alphac是归一化过的，则rho(code)是单位体积内的，这种单位内的量容易搭积分的便车补上计算。对流项没有考虑alphac，因为对流计算过程涉及体积(速度即单位时间内的体积)。

zhe

@bestucan 不好意思,我確實沒看懂,爲什麼一開始就說了alpha1和alpha2是未歸一化的值.

bestucan

@zhe 我说错了，本来想说alpha1+alpha2不等于1。就是括号里的对立面

（如果alpha1 + alpha2 = 1，那么rho = alpha1*rho1 + alpha2*rho2，未找到alpha1和alpha2的关系。）

这个确实不能说 alpha1 和 alpha2 是未归一化的。
但是这个不影响后面的思路，虽然思路很模糊。
其实绝对正确的平均密度公式是：(alpha1 * rho1 + alpha2 * rho2) / (alpha1 + alpha2)。无论alpha1、alpha2有没有归一化都成立。

还有，刚才又回看了一下。你的第一个帖子说

密度rho應該已經經過alpha1*rho1 + alpha2*rho2得到的整體流體密度

如果这个表达式表达整体流体的密度，那是不是就意味着aplha1+alpha2=1?
如果 alpha1+alpha2\=1，那么剩下的体积分数代表的密度不是没算进来。但是代码中并未有限制alpha1和alpha2关系的地方。alphac倒是限制为1-xxx。

不考虑求解器。多相流的平均密度应该是连续相和离散相的平均alphac * rhoc + alphad * rhod，其中alphac + alphad = 1。

MPPIC方法在OpenFOAM里计算好像用的是particle parcel？所有都没有见 alphad。

我说的太外行就不用管了，我没研究过，就是插一句嘴。
附上一个某个帖子的总结：

忘了哪个帖子了，是这个论坛的。

zhe

@bestucan 谢谢，已经说的很好了。我之前想过rho=(alpha1rho1+alpha2rho2)/(alpha1+alpha2).这个我觉得是比较合适的。虽然这样定义，直接就说明了alphac=alpha1+alpha2，而不是从后面的1-theta来的，但是我觉得这不会引起什么问题，因为宏观上来说，alphac就是等于alpha1+alpha2. 这样也满足了alpha1+alpha2+theta=alphac+theta=1. 都是归一后的了。我也试了自带的算例。计算出来的结果与之前定义rho=alpha1rho1+alpha2rho2有差别。因为没有真实实验数据可以对比，所以不敢说哪个就是对的，但是后改的这个，更能理解。单独定义rho后，我也就可以重新考虑alphac*rho这个项。我还是觉得对流项中是要考虑alphac的。希望，有朋友可以指出错误，谢谢！