MPPICInterFoam 中如何通过alphac得到alpha1和alpha2?
-
@tens 这一步应该是用upwind格式求出来alpha1之后然后下面再通过MULES::correct进行修正。如果不进行修正的话，就是我之前提到的MULES::explicitSolve，在MULESTemplates.C里面定义了这个函数：
template<class RdeltaTType, class RhoType, class SpType, class SuType> void Foam::MULES::explicitSolve ( const RdeltaTType& rDeltaT, const RhoType& rho, volScalarField& psi, const surfaceScalarField& phiPsi, const SpType& Sp, const SuType& Su ) { Info<< "MULES: Solving for " << psi.name() << endl; const fvMesh& mesh = psi.mesh(); scalarField& psiIf = psi; const scalarField& psi0 = psi.oldTime(); psiIf = 0.0; fvc::surfaceIntegrate(psiIf, phiPsi); if (mesh.moving()) { psiIf = ( mesh.Vsc0()().field()*rho.oldTime().field() *psi0*rDeltaT/mesh.Vsc()().field() + Su.field() - psiIf )/(rho.field()*rDeltaT - Sp.field()); } else { psiIf = ( rho.oldTime().field()*psi0*rDeltaT + Su.field() - psiIf )/(rho.field()*rDeltaT - Sp.field()); } psi.correctBoundaryConditions(); }其中
psiIf = ( rho.oldTime().field()*psi0*rDeltaT + Su.field() - psiIf )/(rho.field()*rDeltaT - Sp.field());就是我上面提到的方程
alphan应该就是最终求得的新的时间步的体积分数 -
@jasper-0
懂了，现在就是说是红色部分的问题对吧，我看了代码中这部分是通过下面代码求的fvc::surfaceIntegrate(psiIf, phiPsi)而里面的phiPsi则对应的这一部分
通过追踪alphaPhi10 = talphaPhi1Un; //又查到 tmp<surfaceScalarField> talphaPhi1Un ( fvc::flux ( phiCN(), cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(), alphaScheme ) + fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) );而其中的phiCN则
tmp<surfaceScalarField> phiCN(alphaPhic); // Calculate the Crank-Nicolson off-centred volumetric flux if (ocCoeff > 0) { phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime(); }如果不采用CN格式phiCN就等于alphaPhic，又查到
alphacf = fvc::interpolate(alphac); alphaPhic = alphacf*phi;所以本质上图一中红色部分还是通过alphac*alpha1得到的
-
-
-
@tens 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说：
@jasper-0
懂了，现在就是说是红色部分的问题对吧，我看了代码中这部分是通过下面代码求的fvc::surfaceIntegrate(psiIf, phiPsi)而里面的phiPsi则对应的这一部分
通过追踪alphaPhi10 = talphaPhi1Un; //又查到 tmp<surfaceScalarField> talphaPhi1Un ( fvc::flux ( phiCN(), cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(), alphaScheme ) + fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) );而其中的phiCN则
tmp<surfaceScalarField> phiCN(alphaPhic); // Calculate the Crank-Nicolson off-centred volumetric flux if (ocCoeff > 0) { phiCN = cnCoeff*alphaPhic + (1.0 - cnCoeff)*alphaPhic.oldTime(); }如果不采用CN格式phiCN就等于alphaPhic，又查到
alphacf = fvc::interpolate(alphac); alphaPhic = alphacf*phi;所以本质上图一中红色部分还是通过alphac*alpha1得到的
@zhe 看这条消息，code中就是按这个公式求的，alphan*alphac是真实的相体积分数，因此求得的alphan就是占连续相的体积分数
-
@jasper-0 不好意思，最近没看到回复。如果像你说的alpha1*alphac就是真实的体积分数，那么alpha1才应该是整个cell里的alpha1啊，否则不需要再去乘以alphac.
我之前怀疑过，MPPICInterFoam并不是像我们想象的那样，alpha1+alpha2+alpha_p = alphac + alpha_p = 1。而是按照interFoam的alpha1 + alpha2 = 1进行了两相流的计算。随后在动量方程和连续方程中，乘以alphac（这个值与lapha1和alpha2无关，而是用1-alpha_p来的）的效果，以保证真实流过的体积是考虑过粒子占比后的，也能解释方程右边加入的momentum transfer的力。所以可以说OpenFOAM目前做的三相流耦合，还是秉承了粒子在流体中不占体积来的（只是有相互作用力）。
以上如果有不对的，希望可以多交流，感谢！ -
-
@jasper-0 MPPICInterFoam中考虑了体积分数是为了保证质量守恒。你举的例子，alpha1=0.5本身就是占据整个cell的体积分数。而乘以alphac后得到的是实际占据cell的体积分数（因为cell中被粒子占据了一部分）。所以，整个过程中，alpha1一直都是0.5（alpha2=1-alpha1=1-0.5=0.5)这样是符合算法的。而上述例子中最后的0.4+0.4+0.2=1。可是OpenFOAM对于alpha2的求解是通过1-alpha1. 0.4+0.4不等于1。
考虑一个简单的例子，只有粒子和水。alphac=1-theta=0.8。那么如果考虑了粒子的占比，那么水在这里就该是alpha1=alphac=0.8.那么在根据alpha1alphac=0.64?这个值是什么？
然而回到我说的没有考虑粒子占比，对于alpha1的值，一直是占据cell的，那么alpha1就该是等于1（alpha2=1-alpha1=1-1=0. 符合假设前提）。而在实际运动过程中，由于粒子占比，alpha1实际的流通量是alpha1alphac=1*0.8=0.8. 这里考虑了alphac的影响，但是跟alpha1本身的值没有任何关系。
我做了个最简单的验证，上下一共5个网格，最下面的放了一颗粒子。下面3个网格是水，上面是空气。得到的数据是
希望能帮助大家。
-
-
@tens 你的意思我明白，但是将alphac看做单位1是混淆了概念。alphac就是alphac,它是个值，并是有时不是1。这个操作在OpenFOAM中没有。在概念上来说，你说的通，我是认同的。但是在存在MPPICInterFoam中，这种把alphac当作1,是不存在的。alpha1一直以来都是占cell的体积分数，你所说的把alphac当作1，也不是不对，而是当作1，其实就是cell本身。
这里有人会迷糊，说那么alphac到底干啥用？它就是用来体现粒子占据一定体积后，对流体总体积造成了影响，不是来梳理alphac, alpha1, alpha2的关系的。
