of中的多相流求解器
-
@队长别开枪 这个是compressiblemultiphaseFoam的关于求解alpha的源码，从源码看，除了用explicitsolve，在这前还使用了limit以及limitsum 限制，我在想是不是这样子才会使得求解每一相同时还能保证alpha之和为1.
void Foam::multiphaseMixtureThermo::solveAlphas ( const scalar cAlpha ) { static label nSolves=-1; nSolves++; word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); surfaceScalarField phic(mag(phi_/mesh_.magSf())); phic = min(cAlpha*phic, max(phic)); PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size()); int phasei = 0; forAllIter(PtrDictionary<phaseModel>, phases_, phase) { phaseModel& alpha = phase(); phiAlphaCorrs.set ( phasei, new surfaceScalarField ( fvc::flux ( phi_, alpha, alphaScheme ) ) ); surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei]; forAllIter(PtrDictionary<phaseModel>, phases_, phase2) { phaseModel& alpha2 = phase2(); if (&alpha2 == &alpha) continue; surfaceScalarField phir(phic*nHatf(alpha, alpha2)); phiAlphaCorr += fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha, alpharScheme ); } MULES::limit ( 1.0/mesh_.time().deltaT().value(), geometricOneField(), alpha, phi_, phiAlphaCorr, zeroField(), zeroField(), 1, 0, 3, true ); phasei++; } MULES::limitSum(phiAlphaCorrs); rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0); volScalarField sumAlpha ( IOobject ( "sumAlpha", mesh_.time().timeName(), mesh_ ), mesh_, dimensionedScalar("sumAlpha", dimless, 0) ); volScalarField divU(fvc::div(fvc::absolute(phi_, U_))); phasei = 0; forAllIter(PtrDictionary<phaseModel>, phases_, phase) { phaseModel& alpha = phase(); surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei]; phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha); volScalarField::DimensionedInternalField Sp ( IOobject ( "Sp", mesh_.time().timeName(), mesh_ ), mesh_, dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0) ); volScalarField::DimensionedInternalField Su ( IOobject ( "Su", mesh_.time().timeName(), mesh_ ), // Divergence term is handled explicitly to be // consistent with the explicit transport solution divU*min(alpha, scalar(1)) ); { const scalarField& dgdt = alpha.dgdt(); forAll(dgdt, celli) { if (dgdt[celli] < 0.0 && alpha[celli] > 0.0) { Sp[celli] += dgdt[celli]*alpha[celli]; Su[celli] -= dgdt[celli]*alpha[celli]; } else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0) { Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]); } } } forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2) { const phaseModel& alpha2 = phase2(); if (&alpha2 == &alpha) continue; const scalarField& dgdt2 = alpha2.dgdt(); forAll(dgdt2, celli) { if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0) { Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]); Su[celli] += dgdt2[celli]*alpha[celli]; } else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0) { Sp[celli] += dgdt2[celli]*alpha2[celli]; } } } MULES::explicitSolve ( geometricOneField(), alpha, phiAlpha, Sp, Su ); -
@队长别开枪 ，最近做了一些测试，发现多相流求解器为什么求解每一相的同时保证
alpha的总和有界。这里最关键的地方代码是limitsum通量限制器。这个是limitsum的源码，void Foam::MULES::limitSum(UPtrList<scalarField>& phiPsiCorrs) { forAll(phiPsiCorrs[0], facei) { scalar sumPos = 0; scalar sumNeg = 0; for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++) { if (phiPsiCorrs[phasei][facei] > 0) { sumPos += phiPsiCorrs[phasei][facei]; } else { sumNeg += phiPsiCorrs[phasei][facei]; } } scalar sum = sumPos + sumNeg; if (sum > 0 && sumPos > VSMALL) { scalar lambda = -sumNeg/sumPos; for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++) { if (phiPsiCorrs[phasei][facei] > 0) { phiPsiCorrs[phasei][facei] *= lambda; } } } else if (sum < 0 && sumNeg < -VSMALL) { scalar lambda = -sumPos/sumNeg; for (int phasei=0; phasei<phiPsiCorrs.size(); phasei++) { if (phiPsiCorrs[phasei][facei] < 0) { phiPsiCorrs[phasei][facei] *= lambda; } } } } }看完代码后，可以比较清晰的看出，它其实做的工作量就是统计每一个面上的每一相的
alpha的通量，然后进行了修正，得到了修正系数再乘上了alpha的通量，输出。surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei]; phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha)这两行代码是指，最后这个修正的通量加上原来的通量，作为
alpha的通量计算。 -
@mohui OpenFOAM的另一个非官方版本（foam-extend）里alpha方程的源代码比官方的看起来简洁一点，比如求解器
interFoam，extend版本（4.0）中的alphaEqn.H为{ word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); surfaceScalarField phic = mag(phi/mesh.magSf()); phic = min(interface.cAlpha()*phic, max(phic)); surfaceScalarField phir = phic*interface.nHatf(); for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) { surfaceScalarField phiAlpha = fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), alpha1, alpharScheme ); MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; } Info<< "Liquid phase volume fraction = " << alpha1.weightedAverage(mesh.V()).value() << " Min(alpha1) = " << min(alpha1).value() << " Max(alpha1) = " << max(alpha1).value() << endl; }官方版本（4.x）中的为
{ word alphaScheme("div(phi,alpha)"); word alpharScheme("div(phirb,alpha)"); tmp<fv::ddtScheme<scalar>> ddtAlpha ( fv::ddtScheme<scalar>::New ( mesh, mesh.ddtScheme("ddt(alpha)") ) ); // Set the off-centering coefficient according to ddt scheme scalar ocCoeff = 0; if ( isType<fv::EulerDdtScheme<scalar>>(ddtAlpha()) || isType<fv::localEulerDdtScheme<scalar>>(ddtAlpha()) ) { ocCoeff = 0; } else if (isType<fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha())) { if (nAlphaSubCycles > 1) { FatalErrorInFunction << "Sub-cycling is not supported " "with the CrankNicolson ddt scheme" << exit(FatalError); } ocCoeff = refCast<const fv::CrankNicolsonDdtScheme<scalar>>(ddtAlpha()) .ocCoeff(); } else { FatalErrorInFunction << "Only Euler and CrankNicolson ddt schemes are supported" << exit(FatalError); } scalar cnCoeff = 1.0/(1.0 + ocCoeff); // Standard face-flux compression coefficient surfaceScalarField phic(mixture.cAlpha()*mag(phi/mesh.magSf())); // Add the optional isotropic compression contribution if (icAlpha > 0) { phic *= (1.0 - icAlpha); phic += (mixture.cAlpha()*icAlpha)*fvc::interpolate(mag(U)); } surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef(); // Do not compress interface at non-coupled boundary faces // (inlets, outlets etc.) forAll(phic.boundaryField(), patchi) { fvsPatchScalarField& phicp = phicBf[patchi]; if (!phicp.coupled()) { phicp == 0; } } tmp<surfaceScalarField> phiCN(phi); // Calculate the Crank-Nicolson off-centred volumetric flux if (ocCoeff > 0) { phiCN = cnCoeff*phi + (1.0 - cnCoeff)*phi.oldTime(); } if (MULESCorr) { fvScalarMatrix alpha1Eqn ( ( LTS ? fv::localEulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) : fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) ) + fv::gaussConvectionScheme<scalar> ( mesh, phiCN, upwind<scalar>(mesh, phiCN) ).fvmDiv(phiCN, alpha1) ); alpha1Eqn.solve(); Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(" << alpha1.name() << ") = " << min(alpha1).value() << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux()); alphaPhi = talphaPhiUD(); if (alphaApplyPrevCorr && talphaPhiCorr0.valid()) { Info<< "Applying the previous iteration compression flux" << endl; MULES::correct(alpha1, alphaPhi, talphaPhiCorr0.ref(), 1, 0); alphaPhi += talphaPhiCorr0(); } // Cache the upwind-flux talphaPhiCorr0 = talphaPhiUD; alpha2 = 1.0 - alpha1; mixture.correct(); } for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) { surfaceScalarField phir(phic*mixture.nHatf()); tmp<surfaceScalarField> talphaPhiUn ( fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); // Calculate the Crank-Nicolson off-centred alpha flux if (ocCoeff > 0) { talphaPhiUn = cnCoeff*talphaPhiUn + (1.0 - cnCoeff)*alphaPhi.oldTime(); } if (MULESCorr) { tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi); volScalarField alpha10("alpha10", alpha1); MULES::correct(alpha1, talphaPhiUn(), talphaPhiCorr.ref(), 1, 0); // Under-relax the correction for all but the 1st corrector if (aCorr == 0) { alphaPhi += talphaPhiCorr(); } else { alpha1 = 0.5*alpha1 + 0.5*alpha10; alphaPhi += 0.5*talphaPhiCorr(); } } else { alphaPhi = talphaPhiUn; MULES::explicitSolve(alpha1, phiCN, alphaPhi, 1, 0); } alpha2 = 1.0 - alpha1; mixture.correct(); } if (alphaApplyPrevCorr && MULESCorr) { talphaPhiCorr0 = alphaPhi - talphaPhiCorr0; } if ( word(mesh.ddtScheme("ddt(rho,U)")) == fv::EulerDdtScheme<vector>::typeName ) { rhoPhi = alphaPhi*(rho1 - rho2) + phiCN*rho2; } else { if (ocCoeff > 0) { // Calculate the end-of-time-step alpha flux alphaPhi = (alphaPhi - (1.0 - cnCoeff)*alphaPhi.oldTime())/cnCoeff; } // Calculate the end-of-time-step mass flux rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2; } Info<< "Phase-1 volume fraction = " << alpha1.weightedAverage(mesh.Vsc()).value() << " Min(" << alpha1.name() << ") = " << min(alpha1).value() << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; }具体区别我不大清楚，但看起来还是有区别的。
-
@mohui 在 of中的多相流求解器 中说：
我说的错误是这个，里面是没有alpha1，alpha2的
是说 compreisblInterFoam解析 中下面这个方程有误？
\begin{equation}
\mathrm{dgdt}= \left( \frac{\alpha_2}{\rho_2}\frac{\mathrm{D}\rho_2}{\mathrm{D}t} - \frac{\alpha_1}{\rho_1}\frac{\mathrm{D}\rho_1}{\mathrm{D}t} \right)
\end{equation} -
@mohui 非常感谢！原文已更新:cheeky:
有关
MULES的验证（暂时说不上研究，因为基金会已经植入了），我这里分为俩部分验证，第一部分验证已经完毕并且投稿，第二部分验证还在进行中，不过进度很慢毕竟还有别的课题在进行，我只能建议你参考Flux Correct Transport部分数值算法内容，MULES的思想来源于此，只不过基金会将此方法用于多相流有关求解一个方程或者两个方程的验证，应该不是很难，并且结果应该一致。在你之前发的相图上，我没有看出太大差别，是否有曲线图。另外，是否对比过1D
Riemann问题的解析解？ -
@李东岳 东岳老师，你可能还对我提的问题存在误区，我并不是对mules做验证，仅是分析OF中多相流求解器和两相流求解器对与两相流求解的区别。以我目前的测试检验，认为多相流求解器每一相求解并能够同时保证总的alpha之和等于的1的关键在于limitsum这个限制器，关于这个限制器的代码我也上传了，通过代码可以分析的出来，就是在mules求解之前对alphaphi进行了修正，以保证流入流出单元的alphaphi总和为0。其次，就如我前面所言，如果单纯的mules求解每一相是无法保证alpha1+alpha2=1的，虽然我给的图上没有很大区别，但是我输出alpha1+alpha2最大值为2，这是不对的。以上是我这几天的发现。
-
@李东岳 在 of中的多相流求解器 中说：
@mohui 非常感谢！原文已更新:cheeky:
有关
MULES的验证（暂时说不上研究，因为基金会已经植入了），我这里分为俩部分验证，第一部分验证已经完毕并且投稿，第二部分验证还在进行中，不过进度很慢毕竟还有别的课题在进行，我只能建议你参考Flux Correct Transport部分数值算法内容，MULES的思想来源于此，只不过基金会将此方法用于多相流有关求解一个方程或者两个方程的验证，应该不是很难，并且结果应该一致。在你之前发的相图上，我没有看出太大差别，是否有曲线图。另外，是否对比过1D
Riemann问题的解析解？岳哥，公式20更新了，公式19却没更新。建议更新一下，否则看着怪怪的。
-
确实！好眼力。
@李东岳
