MPPICInterFoam 中如何通过alphac得到alpha1和alpha2?
-
@tens 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说:
通过alphaEqn.H中的206-226行有如下代码,这几行代码应该是求解出来了乘以alphac之后的alpha1,那alpha2怎么是1-alpha1呢?不是应该alphac-alpha1吗?还望高手解答一下啊
else { alphaPhi10 = talphaPhi1Un; MULES::explicitSolve ( alphac, alpha1, phiCN, alphaPhi10, Sp, (Su + divU*min(alpha1(), scalar(1)))(), oneField(), zeroField() ); } alpha2 = 1.0 - alpha1; mixture.correct(); }
-
@jasper-0 在 MPPICInterFoam 中如何通过alphac得到alpha1和alpha2? 中说:
@tens 感谢大佬的解答!也就是说通过以下代码求得的alpha1是占连续相的体积分数?那最后结果跑出来后要想得到总的alpha1是不是要在paraview里面做计算alpha1*alphac?
MULES::explicitSolve ( alphac, alpha1, phiCN, alphaPhi10, Sp, (Su + divU*min(alpha1(), scalar(1)))(), oneField(), zeroField() ); }
我发现interFoam中alpha1(即alpha1)的求解对应下面这个公式,
而MPPICInterFoam中通过以上代码求得的alphan为
公式右边第二项 是用alpha占连续相的体积分数算出来的,
而 和 都是总的相分数,所以此方程右边第二项是不是应该乘以个alphac呢?即
还望大佬解答一下,万分感谢! -
@jasper-0 这个公式是哪里看到的,字母对应的意义不理解,在程序中相体积分数应该是通过控制方程求的吧
fvScalarMatrix alpha1Eqn ( ( LTS ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alphac, alpha1) : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) ) + fv::gaussConvectionScheme<scalar> ( mesh, phiCN, upwind<scalar>(mesh, phiCN) ).fvmDiv(phiCN, alpha1) - fvm::Sp(fvc::ddt(alphac) + fvc::div(phiCN), alpha1) == Su + fvm::Sp(Sp + divU, alpha1) ); alpha1Eqn.solve();
-
@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得到的