of中的多相流求解器
-
@李东岳 ，首先我是看了您写的关于compressibleInterFoam求解器解析，然后推导完之后发现了个小小的错误。dgdt里面是没有alpha的。
然后这个是我推导的可压缩两相方程：
，
根据这个方程结合源码，这是修改compressibleInterFoam中的alpha方程的代码：for (int gCorr=0; gCorr<nAlphaCorr; gCorr++) { volScalarField::Internal Sp ( IOobject ( "Sp", runTime.timeName(), mesh ), mesh, dimensionedScalar("Sp", dgdt.dimensions(), 0.0) ); volScalarField::Internal Su ( IOobject ( "Su", runTime.timeName(), mesh ), divU*min(alpha1, scalar(1)) ); forAll(dgdt, celli) { if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0) { Sp[celli] -= dgdt[celli]*alpha1[celli]; Su[celli] += dgdt[celli]*alpha1[celli]; } else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0) { Sp[celli] += dgdt[celli]*(max(1.0 - alpha1[celli],scalar(0))); } } surfaceScalarField alphaPhi1 ( fvc::flux ( phi, alpha1, alphaScheme ) + fvc::flux ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); MULES::explicitSolve ( geometricOneField(), alpha1, phi, alphaPhi1, Sp, Su, 1, 0 ); //slove for alpha2 volScalarField::Internal Sp1 ( IOobject ( "Sp1", runTime.timeName(), mesh ), mesh, dimensionedScalar("Sp1", dgdt.dimensions(), 0.0) ); volScalarField::Internal Su1 ( IOobject ( "Su1", runTime.timeName(), mesh ), divU*min(alpha2, scalar(1)) ); forAll(dgdt, celli) { if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0) { Sp1[celli] -= dgdt[celli]*alpha1[celli]; } else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0) { Sp1[celli] += dgdt[celli]*alpha2[celli]; Su1[celli] -= dgdt[celli]*(max(1.0 - alpha1[celli],scalar(0))); } } surfaceScalarField alphaPhi2 ( fvc::flux ( phi, alpha2, alphaScheme ) + fvc::flux ( -fvc::flux(phir, alpha1, alpharScheme), alpha2, alpharScheme ) ); MULES::explicitSolve ( geometricOneField(), alpha2, phi, alphaPhi2, Sp1, Su1, 1, 0 ); -
-
@mohui Hello, 不好意思现在才回复你。OpenFOAM的MULES求解器我不大熟悉，我自己本身是做多面体网格上的PLIC算法的，最后评估算法优劣的参数中就有一个mass/volume conservation error，这个值越小越好，如果算法能够保证这个参数恒等于零，那么分别求解alpha1和alpha2没有什么问题，但现实中算法误差和数值误差的存在，这个值不可能恒等于零，所以一般次相alpha2直接用1-alpha1计算，如果这个时候分别计算alpha1和alpha2，我觉得会叠加误差，以至于出现你遇到的alpha1+alpha2>1的情况，更新密度粘度场后误差再传进动量方程和压力泊松方程，最后造成求解发散。这些是我的看法。
-
@队长别开枪 这个是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却没更新。建议更新一下，否则看着怪怪的。
-
确实！好眼力。
@李东岳
，我说的错误是这个，里面是没有alpha1，alpha2的。对，算的结果是alpha1+alpha2最大值为2.。
左边是原程序，右边是我改的程序。两者还是有点区别的。
