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}
，我说的错误是这个，里面是没有alpha1，alpha2的。对，算的结果是alpha1+alpha2最大值为2.。
左边是原程序，右边是我改的程序。两者还是有点区别的。