diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C index 2fa154d935..13d95ceba0 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C @@ -45,6 +45,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() { volScalarField& p(p_); + const phaseSystem::phaseModelPartialList& movingPhases + ( + fluid.movingPhases() + ); + // Face volume fractions PtrList alphafs(phases.size()); forAll(phases, phasei) @@ -57,21 +62,25 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } // Diagonal coefficients - rAUs.clear(); - rAUs.setSize(phases.size()); - - PtrList rAUfs(phases.size()); - + rAs.clear(); + if (fluid.implicitPhasePressure()) { - PtrList Kds(fluid.Kds()); + rAs.setSize(phases.size()); + } - forAll(fluid.movingPhases(), movingPhasei) + PtrList> invADs; + PtrList> invADfs; + { + PtrList As(movingPhases.size()); + PtrList Afs(movingPhases.size()); + forAll(movingPhases, movingPhasei) { - const phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const phaseModel& phase = movingPhases[movingPhasei]; const volScalarField& alpha = phase; - volScalarField AU + As.set ( + movingPhasei, UEqns[phase.index()].A() + byDt ( @@ -80,153 +89,189 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ) ); - if (Kds.set(phase.index())) - { - AU += Kds[phase.index()]; - } - - rAUs.set + Afs.set ( - phase.index(), - new volScalarField - ( - IOobject::groupName("rAU", phase.name()), - 1/AU - ) + movingPhasei, + fvc::interpolate(As[movingPhasei]) ); - rAUfs.set(phase.index(), 1/fvc::interpolate(AU)); + if (fluid.implicitPhasePressure()) + { + rAs.set + ( + phase.index(), + new volScalarField + ( + IOobject::groupName("rA", phase.name()), + 1/As[movingPhasei] + ) + ); + } } - fluid.fillFields("rAU", dimTime/dimDensity, rAUs); - fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs); + + fluid.invADs(As, invADs, invADfs); } - // Phase diagonal coefficients - PtrList alpharAUfs(phases.size()); - forAll(phases, phasei) - { - const phaseModel& phase = phases[phasei]; - const volScalarField& alpha = phase; - - alpharAUfs.set - ( - phasei, - ( - fvc::interpolate(max(alpha, phase.residualAlpha())) - *rAUfs[phasei] - ).ptr() - ); - } + volScalarField rho("rho", fluid.rho()); // Explicit force fluxes - PtrList Fs(fluid.Fs()); + PtrList alphaByADfs; + PtrList FgByADfs; + { + PtrList Ffs(fluid.Fs()); + + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + PtrList lalphafs(movingPhases.size()); + PtrList Fgfs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + lalphafs.set + ( + movingPhasei, + fvc::interpolate(max(alpha, phase.residualAlpha())) + ); + + Fgfs.set + ( + movingPhasei, + ( + Ffs[phase.index()] + + lalphafs[movingPhasei] + *( + ghSnGradRho + - (fvc::interpolate(phase.rho() - rho)) + *(buoyancy.g & mesh.Sf()) + - fluid.surfaceTension(phase)*mesh.magSf() + ) + ).ptr() + ); + } + + alphaByADfs = invADfs & lalphafs; + FgByADfs = invADfs & Fgfs; + } + // Mass transfer rates PtrList dmdts(fluid.dmdts()); PtrList d2mdtdps(fluid.d2mdtdps()); + // --- Optional momentum predictor + if (predictMomentum) + { + InfoInFunction << endl; + PtrList HbyADs; + { + PtrList Hs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + Hs.set + ( + movingPhasei, + UEqns[phase.index()].H() + + byDt + ( + max(phase.residualAlpha() - alpha, scalar(0)) + *phase.rho() + ) + *phase.U()().oldTime() + ); + } + + HbyADs = invADs & Hs; + } + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh); + } + + const surfaceScalarField mSfGradp(-mesh.magSf()*fvc::snGrad(p_rgh)); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = fluid_.movingPhases()[movingPhasei]; + + phase.URef() = + HbyADs[movingPhasei] + + fvc::reconstruct + ( + alphaByADfs[movingPhasei]*mSfGradp + - FgByADfs[movingPhasei] + ); + + phase.URef().correctBoundaryConditions(); + fvConstraints().constrain(phase.URef()); + } + } + + // --- Pressure corrector loop while (pimple.correct()) { - volScalarField rho("rho", fluid.rho()); - - // Correct p_rgh for consistency with p and the updated densities - p_rgh = p - rho*buoyancy.gh; - // Correct fixed-flux BCs to be consistent with the velocity BCs fluid_.correctBoundaryFlux(); - // Combined buoyancy and force fluxes - PtrList phigFs(phases.size()); + PtrList HbyADs; + PtrList phiHbyADs; { - const surfaceScalarField ghSnGradRho - ( - "ghSnGradRho", - buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() - ); + // Predicted velocities and fluxes for each phase + PtrList Hs(movingPhases.size()); + PtrList phiHs(movingPhases.size()); - forAll(phases, phasei) - { - const phaseModel& phase = phases[phasei]; - - phigFs.set - ( - phasei, - ( - ( - ghSnGradRho - - (fvc::interpolate(phase.rho() - rho)) - *(buoyancy.g & mesh.Sf()) - - fluid.surfaceTension(phase)*mesh.magSf() - ) - ).ptr() - ); - } - } - - // Predicted velocities and fluxes for each phase - PtrList HbyAs(phases.size()); - PtrList phiHbyAs(phases.size()); - { // Correction force fluxes PtrList ddtCorrs(fluid.ddtCorrs()); - forAll(fluid.movingPhases(), movingPhasei) + forAll(movingPhases, movingPhasei) { - const phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const phaseModel& phase = movingPhases[movingPhasei]; const volScalarField& alpha = phase; - const label phasei = phase.index(); - const volVectorField H + Hs.set ( - constrainH + movingPhasei, + UEqns[phase.index()].H() + + byDt ( - UEqns[phasei].H() - + byDt - ( - max(phase.residualAlpha() - alpha, scalar(0)) - *phase.rho() - ) - *phase.U()().oldTime(), - rAUs[phasei], - phase.U(), - p_rgh + max(phase.residualAlpha() - alpha, scalar(0)) + *phase.rho() ) + *phase.U()().oldTime() ); - HbyAs.set(phasei, rAUs[phasei]*H); - - phiHbyAs.set + phiHs.set ( - phasei, - new surfaceScalarField - ( - IOobject::groupName("phiHbyA", phase.name()), - rAUfs[phasei]*(fvc::flux(H) + ddtCorrs[phasei]) - - alpharAUfs[phasei]*phigFs[phasei] - - rAUfs[phasei]*Fs[phasei] - ) + movingPhasei, + fvc::flux(Hs[movingPhasei]) + ddtCorrs[phase.index()] ); } + + HbyADs = invADs & Hs; + phiHbyADs = invADfs & phiHs; } - fluid.fillFields("HbyA", dimVelocity, HbyAs); - fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs); - // Add explicit drag forces and fluxes - PtrList KdUs(fluid.KdUs()); - PtrList KdPhis(fluid.KdPhis()); - - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - if (KdUs.set(phasei)) - { - HbyAs[phasei] -= rAUs[phasei]*KdUs[phasei]; - } + const phaseModel& phase = movingPhases[movingPhasei]; - if (KdPhis.set(phasei)) - { - phiHbyAs[phasei] -= rAUfs[phasei]*KdPhis[phasei]; - } + constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh); + constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh); + + phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei]; } // Total predicted flux @@ -236,28 +281,28 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ( "phiHbyA", runTime.name(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE + mesh ), mesh, dimensionedScalar(dimFlux, 0) ); - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - phiHbyA += alphafs[phasei]*phiHbyAs[phasei]; + const phaseModel& phase = movingPhases[movingPhasei]; + + phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei]; } MRF.makeRelative(phiHbyA); - fvc::makeRelative(phiHbyA, fluid.movingPhases()[0].U()); + fvc::makeRelative(phiHbyA, movingPhases[0].U()); - // Construct pressure "diffusivity" - surfaceScalarField rAUf + // Pressure "diffusivity" + surfaceScalarField rAf ( IOobject ( - "rAUf", + "rAf", runTime.name(), mesh ), @@ -265,9 +310,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() dimensionedScalar(dimTime/dimDensity, 0) ); - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - rAUf += alphafs[phasei]*alpharAUfs[phasei]; + const phaseModel& phase = movingPhases[movingPhasei]; + + rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei]; } // Update the fixedFluxPressure BCs to ensure flux consistency @@ -279,11 +326,12 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ); phib = 0; - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - const phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid_.movingPhases()[movingPhasei]; + phib += - alphafs[phasei].boundaryField() + alphafs[phase.index()].boundaryField() *phase.phi()().boundaryField(); } @@ -292,7 +340,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - phib - )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + )/(mesh.magSf().boundaryField()*rAf.boundaryField()) ); } @@ -309,7 +357,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() fvScalarMatrix pEqnIncomp ( fvc::div(phiHbyA) - - fvm::laplacian(rAUf, p_rgh) + - fvm::laplacian(rAf, p_rgh) ); // Solve @@ -340,86 +388,76 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() { phi_ = phiHbyA + pEqnIncomp.flux(); - surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf); - forAll(fluid.movingPhases(), movingPhasei) + forAll(movingPhases, movingPhasei) { phaseModel& phase = fluid_.movingPhases()[movingPhasei]; phase.phiRef() = - phiHbyAs[phase.index()] - + alpharAUfs[phase.index()]*mSfGradp; + phiHbyADs[movingPhasei] + + alphaByADfs[movingPhasei]*mSfGradp; // Set the phase dilatation rate phase.divU(-pEqnComps[phase.index()] & p_rgh); + + MRF.makeRelative(phase.phiRef()); + fvc::makeRelative(phase.phiRef(), phase.U()); } // Optionally relax pressure for velocity correction p_rgh.relax(); - mSfGradp = pEqnIncomp.flux()/rAUf; + mSfGradp = pEqnIncomp.flux()/rAf; - if (!dragCorrection) + if (dragCorrection) { - forAll(fluid.movingPhases(), movingPhasei) + PtrList dragCorrs(movingPhases.size()); + PtrList dragCorrfs(movingPhases.size()); + fluid.dragCorrs(dragCorrs, dragCorrfs); + + PtrList dragCorrByADs + ( + invADs & dragCorrs + ); + + PtrList dragCorrByADfs + ( + invADfs & dragCorrfs + ); + + forAll(movingPhases, movingPhasei) { phaseModel& phase = fluid_.movingPhases()[movingPhasei]; - const label phasei = phase.index(); phase.URef() = - HbyAs[phasei] + HbyADs[movingPhasei] + fvc::reconstruct ( - alpharAUfs[phasei]*(mSfGradp - phigFs[phasei]) - - rAUfs[phasei]*Fs[phasei] + alphaByADfs[movingPhasei]*mSfGradp + - FgByADfs[movingPhasei] + + dragCorrByADfs[movingPhasei] + ) + - dragCorrByADs[movingPhasei]; + } + } + else + { + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = fluid_.movingPhases()[movingPhasei]; + + phase.URef() = + HbyADs[movingPhasei] + + fvc::reconstruct + ( + alphaByADfs[movingPhasei]*mSfGradp + - FgByADfs[movingPhasei] ); } } - else - { - PtrList dragCorrs(phases.size()); - PtrList dragCorrfs(phases.size()); - fluid.dragCorrs(dragCorrs, dragCorrfs); - forAll(fluid.movingPhases(), movingPhasei) - { - phaseModel& phase = fluid_.movingPhases()[movingPhasei]; - const label phasei = phase.index(); - - phase.URef() = - HbyAs[phasei] - + fvc::reconstruct - ( - alpharAUfs[phasei]*(mSfGradp - phigFs[phasei]) - + rAUfs[phasei]*(dragCorrfs[phasei] - Fs[phasei]) - ) - - rAUs[phasei]*dragCorrs[phasei]; - } - } - - if (partialElimination) - { - fluid_.partialElimination - ( - rAUs, - KdUs, - alphafs, - rAUfs, - KdPhis - ); - } - else - { - forAll(fluid.movingPhases(), movingPhasei) - { - phaseModel& phase = fluid_.movingPhases()[movingPhasei]; - - MRF.makeRelative(phase.phiRef()); - fvc::makeRelative(phase.phiRef(), phase.U()); - } - } - - forAll(fluid.movingPhases(), movingPhasei) + forAll(movingPhases, movingPhasei) { phaseModel& phase = fluid_.movingPhases()[movingPhasei]; @@ -472,11 +510,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() } UEqns.clear(); - - if (!fluid.implicitPhasePressure()) - { - rAUs.clear(); - } } diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C index 931481d5ba..b10d5dd223 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C @@ -45,9 +45,13 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() { volScalarField& p(p_); + const phaseSystem::phaseModelPartialList& movingPhases + ( + fluid.movingPhases() + ); + // Face volume fractions PtrList alphafs(phases.size()); - PtrList alphaRho0fs(phases.size()); forAll(phases, phasei) { const phaseModel& phase = phases[phasei]; @@ -55,65 +59,113 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() alphafs.set(phasei, fvc::interpolate(alpha).ptr()); alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); - - alphaRho0fs.set - ( - phasei, - ( - fvc::interpolate - ( - max(alpha.oldTime(), phase.residualAlpha()) - *phase.rho().oldTime() - ) - ).ptr() - ); } // Diagonal coefficients - rAUfs.clear(); - rAUfs.setSize(phases.size()); + rAs.clear(); + if (fluid.implicitPhasePressure()) { - PtrList KdVmfs(fluid.KdVmfs()); + rAs.setSize(phases.size()); + } + + PtrList> invADfs; + { + PtrList Vmfs(fluid.Vmfs()); + PtrList Afs(movingPhases.size()); forAll(fluid.movingPhases(), movingPhasei) { const phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const volScalarField& alpha = phase; - rAUfs.set + const volScalarField A ( - phase.index(), + byDt + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho().oldTime() + ) + + UEqns[phase.index()].A() + ); + + if (fluid.implicitPhasePressure()) + { + rAs.set + ( + phase.index(), + new volScalarField + ( + IOobject::groupName("rA", phase.name()), + 1/A + ) + ); + } + + Afs.set + ( + movingPhasei, new surfaceScalarField ( - IOobject::groupName("rAUf", phase.name()), - 1.0 - /( - byDt(alphaRho0fs[phase.index()]) - + fvc::interpolate(UEqns[phase.index()].A()) - + KdVmfs[phase.index()] - ) + IOobject::groupName("rAf", phase.name()), + fvc::interpolate(A) + ) + ); + + if (Vmfs.set(phase.index())) + { + Afs[movingPhasei] += Vmfs[phase.index()]; + } + } + + invADfs = fluid.invADfs(Afs); + } + + volScalarField rho("rho", fluid.rho()); + + // Phase diagonal coefficients + PtrList alphaByADfs; + PtrList FgByADfs; + { + // Explicit force fluxes + PtrList Ffs(fluid.Ffs()); + + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + PtrList lalphafs(movingPhases.size()); + PtrList Fgfs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + lalphafs.set + ( + movingPhasei, + max(alphafs[phase.index()], phase.residualAlpha()) + ); + + Fgfs.set + ( + movingPhasei, + Ffs[phase.index()] + + lalphafs[movingPhasei] + *( + ghSnGradRho + - (fvc::interpolate(phase.rho() - rho)) + *(buoyancy.g & mesh.Sf()) + - fluid.surfaceTension(phase)*mesh.magSf() ) ); } - } - fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs); - // Phase diagonal coefficients - PtrList alpharAUfs(phases.size()); - forAll(phases, phasei) - { - const phaseModel& phase = phases[phasei]; - alpharAUfs.set - ( - phase.index(), - ( - max(alphafs[phase.index()], phase.residualAlpha()) - *rAUfs[phase.index()] - ).ptr() - ); + alphaByADfs = invADfs & lalphafs; + FgByADfs = invADfs & Fgfs; } - // Explicit force fluxes - PtrList Ffs(fluid.Ffs()); // Mass transfer rates PtrList dmdts(fluid.dmdts()); @@ -122,87 +174,49 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() // --- Pressure corrector loop while (pimple.correct()) { - volScalarField rho("rho", fluid.rho()); - - // Correct p_rgh for consistency with p and the updated densities - p_rgh = p - rho*buoyancy.gh; - // Correct fixed-flux BCs to be consistent with the velocity BCs fluid_.correctBoundaryFlux(); - // Combined buoyancy and force fluxes - PtrList phigFs(phases.size()); - { - const surfaceScalarField ghSnGradRho - ( - "ghSnGradRho", - buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() - ); - - forAll(phases, phasei) - { - const phaseModel& phase = phases[phasei]; - - phigFs.set - ( - phasei, - ( - alpharAUfs[phasei] - *( - ghSnGradRho - - (fvc::interpolate(phase.rho() - rho)) - *(buoyancy.g & mesh.Sf()) - - fluid.surfaceTension(phase)*mesh.magSf() - ) - ).ptr() - ); - - if (Ffs.set(phasei)) - { - phigFs[phasei] += rAUfs[phasei]*Ffs[phasei]; - } - } - } - // Predicted fluxes for each phase - PtrList phiHbyAs(phases.size()); - forAll(fluid.movingPhases(), movingPhasei) + PtrList phiHbyADs; { - const phaseModel& phase = fluid.movingPhases()[movingPhasei]; + PtrList phiHs(movingPhases.size()); - phiHbyAs.set - ( - phase.index(), - constrainPhiHbyA + forAll(fluid.movingPhases(), movingPhasei) + { + const phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const volScalarField& alpha = phase; + + phiHs.set ( - rAUfs[phase.index()] - *( - fvc::flux(UEqns[phase.index()].H()) - + alphaRho0fs[phase.index()] + movingPhasei, + ( + fvc::interpolate + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho().oldTime() + ) *byDt ( phase.Uf().valid() ? (mesh.Sf() & phase.Uf()().oldTime()) : MRF.absolute(phase.phi()().oldTime()) ) + + fvc::flux(UEqns[phase.index()].H()) ) - - phigFs[phase.index()], - phase.U(), - p_rgh - ) - ); - } - fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs); - - // Add explicit drag forces and fluxes - PtrList KdPhifs(fluid.KdPhifs()); - - forAll(phases, phasei) - { - if (KdPhifs.set(phasei)) - { - phiHbyAs[phasei] -= rAUfs[phasei]*KdPhifs[phasei]; + ); } + + phiHbyADs = invADfs & phiHs; + } + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh); + + phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei]; } // Total predicted flux @@ -220,20 +234,22 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() dimensionedScalar(dimFlux, 0) ); - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - phiHbyA += alphafs[phasei]*phiHbyAs[phasei]; + const phaseModel& phase = movingPhases[movingPhasei]; + + phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei]; } MRF.makeRelative(phiHbyA); fvc::makeRelative(phiHbyA, phases[0].U()); // Construct pressure "diffusivity" - surfaceScalarField rAUf + surfaceScalarField rAf ( IOobject ( - "rAUf", + "rAf", runTime.name(), mesh ), @@ -241,12 +257,12 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), 0) ); - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - rAUf += alphafs[phasei]*alpharAUfs[phasei]; - } + const phaseModel& phase = movingPhases[movingPhasei]; - rAUf = mag(rAUf); + rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei]; + } // Update the fixedFluxPressure BCs to ensure flux consistency { @@ -257,11 +273,12 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() ); phib = 0; - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { - const phaseModel& phase = phases[phasei]; + phaseModel& phase = fluid_.movingPhases()[movingPhasei]; + phib += - alphafs[phasei].boundaryField() + alphafs[phase.index()].boundaryField() *phase.phi()().boundaryField(); } @@ -270,7 +287,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() p_rgh.boundaryFieldRef(), ( phiHbyA.boundaryField() - phib - )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + )/(mesh.magSf().boundaryField()*rAf.boundaryField()) ); } @@ -287,7 +304,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() fvScalarMatrix pEqnIncomp ( fvc::div(phiHbyA) - - fvm::laplacian(rAUf, p_rgh) + - fvm::laplacian(rAf, p_rgh) ); // Solve @@ -318,39 +335,28 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() { phi_ = phiHbyA + pEqnIncomp.flux(); - surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf); forAll(fluid.movingPhases(), movingPhasei) { phaseModel& phase = fluid_.movingPhases()[movingPhasei]; + const label phasei = phase.index(); phase.phiRef() = - phiHbyAs[phase.index()] - + alpharAUfs[phase.index()]*mSfGradp; + phiHbyADs[movingPhasei] + + alphaByADfs[movingPhasei]*mSfGradp; // Set the phase dilatation rate - phase.divU(-pEqnComps[phase.index()] & p_rgh); - } - - if (partialElimination) - { - fluid_.partialEliminationf(rAUfs, alphafs, KdPhifs); - } - else - { - forAll(fluid.movingPhases(), movingPhasei) - { - phaseModel& phase = fluid_.movingPhases()[movingPhasei]; - - MRF.makeRelative(phase.phiRef()); - fvc::makeRelative(phase.phiRef(), phase.U()); - } + phase.divU(-pEqnComps[phasei] & p_rgh); } forAll(fluid.movingPhases(), movingPhasei) { phaseModel& phase = fluid_.movingPhases()[movingPhasei]; + MRF.makeRelative(phase.phiRef()); + fvc::makeRelative(phase.phiRef(), phase.U()); + phase.URef() = fvc::reconstruct ( fvc::absolute(MRF.absolute(phase.phi()), phase.U()) @@ -405,11 +411,6 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() } UEqns.clear(); - - if (!fluid.implicitPhasePressure()) - { - rAUfs.clear(); - } } diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C b/applications/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C index c7deefbb46..40039a0186 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C @@ -40,8 +40,6 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor() phaseSystem::momentumTransferTable& momentumTransfer(momentumTransferPtr()); - const PtrList Kds(fluid.Kds()); - forAll(fluid.movingPhases(), movingPhasei) { phaseModel& phase = fluid.movingPhases()[movingPhasei]; @@ -59,7 +57,6 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor() == *momentumTransfer[phase.name()] + fvModels().source(alpha, rho, U) - // - fvm::Sp(Kds[phase.index()], U) ) ); diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C index 56cac2dbb9..9102336450 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C @@ -51,15 +51,15 @@ void Foam::solvers::multiphaseEuler::readControls(const bool construct) if (construct || mesh.solution().modified()) { + predictMomentum = + pimple.dict().lookupOrDefault("momentumPredictor", false); + faceMomentum = pimple.dict().lookupOrDefault("faceMomentum", false); dragCorrection = pimple.dict().lookupOrDefault("dragCorrection", false); - partialElimination = - pimple.dict().lookupOrDefault("partialElimination", false); - nEnergyCorrectors = pimple.dict().lookupOrDefault("nEnergyCorrectors", 1); } @@ -98,6 +98,11 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh) : fluidSolver(mesh), + predictMomentum + ( + pimple.dict().lookupOrDefault("momentumPredictor", false) + ), + faceMomentum ( pimple.dict().lookupOrDefault("faceMomentum", false) @@ -108,11 +113,6 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh) pimple.dict().lookupOrDefault("dragCorrection", false) ), - partialElimination - ( - pimple.dict().lookupOrDefault("partialElimination", false) - ), - nEnergyCorrectors ( pimple.dict().lookupOrDefault("nEnergyCorrectors", 1) @@ -247,7 +247,7 @@ void Foam::solvers::multiphaseEuler::prePredictor() { if (pimple.thermophysics() || pimple.flow()) { - fluid_.solve(rAUs, rAUfs); + fluid_.solve(rAs); fluid_.correct(); fluid_.correctContinuityError(); } diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H index bb0d19ad4f..4b4a692495 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H +++ b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H @@ -77,6 +77,10 @@ protected: // Controls + //- Momentum equation predictor switch + // Defaults to false + Switch predictMomentum; + //- Cell/face momentum equation switch // Defaults to false, i.e. uses the cell momentum equation Switch faceMomentum; @@ -85,10 +89,6 @@ protected: // Defaults to false Switch dragCorrection; - //- Partial elimination drag contribution optimisation - // Defaults to false - Switch partialElimination; - //- Number of energy correctors // Used to improve stability of phase-change sibulations // Defaults to 1 @@ -144,11 +144,7 @@ protected: //- Temporary storage for the reciprocal momentum equation diagonal // Used by the phase-fraction predictor and pressure corrector - PtrList rAUs; - - //- Temporary storage for the reciprocal momentum equation diagonal - // Used by the phase-fraction predictor and face pressure corrector - PtrList rAUfs; + PtrList rAs; //- Stored divU from the previous mesh so that it can be // mapped and used in correctPhi to ensure the corrected phi diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index 8450bd912b..2d1df14ea6 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -31,6 +31,8 @@ License #include "wallLubricationModel.H" #include "turbulentDispersionModel.H" +#include "scalarMatrices.H" + #include "fvmDdt.H" #include "fvmDiv.H" #include "fvmSup.H" @@ -474,36 +476,9 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() template Foam::PtrList -Foam::MomentumTransferPhaseSystem::KdVmfs() const +Foam::MomentumTransferPhaseSystem::Vmfs() const { - PtrList KdVmfs(this->phaseModels_.size()); - - // Add the implicit part of the drag force - forAllConstIter(KdfTable, Kdfs_, KdfIter) - { - const surfaceScalarField& Kf(*KdfIter()); - const phaseInterface interface(*this, KdfIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - const phaseModel& phase = iter(); - const phaseModel& otherPhase = iter.otherPhase(); - - const surfaceScalarField alphaf - ( - fvc::interpolate(otherPhase) - ); - - addField - ( - phase, - "KdVmf", - (alphaf/max(alphaf, otherPhase.residualAlpha())) - *Kf, - KdVmfs - ); - } - } + PtrList Vmfs(this->phaseModels_.size()); // Add the implicit part of the virtual mass force forAllConstIter(VmTable, Vms_, VmIter) @@ -522,13 +497,11 @@ Foam::MomentumTransferPhaseSystem::KdVmfs() const *Vm ); - addField(phase, "KdVmf", byDt(fvc::interpolate(VmPhase)), KdVmfs); + addField(phase, "Vmf", byDt(fvc::interpolate(VmPhase)), Vmfs); } } - this->fillFields("KdVmf", dimDensity/dimTime, KdVmfs); - - return KdVmfs; + return Vmfs; } @@ -650,11 +623,6 @@ Foam::MomentumTransferPhaseSystem::Fs() const addField(interface.phase2(), "F", Df*snGradAlpha2By12, Fs); } - if (this->fillFields_) - { - this->fillFields("F", dimArea*dimDensity*dimAcceleration, Fs); - } - return Fs; } @@ -813,22 +781,121 @@ Foam::MomentumTransferPhaseSystem::Ffs() const addField(interface.phase2(), "F", Df*snGradAlpha2By12, Ffs); } - if (this->fillFields_) - { - this->fillFields("Ff", dimArea*dimDensity*dimAcceleration, Ffs); - } - return Ffs; } template -Foam::PtrList -Foam::MomentumTransferPhaseSystem::KdPhis() const +void Foam::MomentumTransferPhaseSystem::invADs +( + List>& ADs +) const { - PtrList KdPhis(this->phaseModels_.size()); + const label n = ADs.size(); + + scalarSquareMatrix AD(n); + scalarField source(n); + labelList pivotIndices(n); + + forAll(ADs[0][0], ci) + { + for (label i=0; i +template class PatchField, class GeoMesh> +void Foam::MomentumTransferPhaseSystem::invADs +( + PtrList>>& ADs +) const +{ + const label n = ADs.size(); + + List> ADps(n); + + for (label i=0; i +void Foam::MomentumTransferPhaseSystem::invADs +( + const PtrList& As, + PtrList>& invADs, + PtrList>& invADfs +) const +{ + const label n = As.size(); + + invADs.setSize(n); + invADfs.setSize(n); + + forAll(invADs, i) + { + invADs.set(i, new PtrList(n)); + invADs[i].set(i, As[i].clone()); + + invADfs.set(i, new PtrList(n)); + invADfs[i].set(i, fvc::interpolate(As[i])); + } + + labelList movingPhases(this->phases().size(), -1); + + forAll(this->movingPhases(), movingPhasei) + { + movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei; + } - // Add the explicit part of the drag force forAllConstIter(KdTable, Kds_, KdIter) { const volScalarField& K(*KdIter()); @@ -839,46 +906,91 @@ Foam::MomentumTransferPhaseSystem::KdPhis() const const phaseModel& phase = iter(); const phaseModel& otherPhase = iter.otherPhase(); - addField - ( - phase, - "KdPhi", - fvc::interpolate + const label i = movingPhases[phase.index()]; + + if (i != -1) + { + const volScalarField Kij ( - -(otherPhase/max(otherPhase, otherPhase.residualAlpha())) - *K - ) - *fvc::absolute - ( - this->MRF().absolute(otherPhase.phi()), - otherPhase.U() - ), - KdPhis - ); + (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*K + ); + + const surfaceScalarField Kijf(fvc::interpolate(Kij)); + + invADs[i][i] += Kij; + invADfs[i][i] += Kijf; + + const label j = movingPhases[otherPhase.index()]; + + if (j != -1) + { + invADs[i].set(j, -Kij); + invADfs[i].set(j, -Kijf); + } + } } } - if (this->fillFields_) + for (label i=0; ifillFields - ( - "KdPhi", - dimArea*dimDensity*dimAcceleration, - KdPhis - ); + for (label j=0; jmesh(), + dimensionedScalar(As[0].dimensions(), 0) + ) + ); + + invADfs[i].set + ( + j, + surfaceScalarField::New + ( + "0", + this->mesh(), + dimensionedScalar(As[0].dimensions(), 0) + ) + ); + } + } } - return KdPhis; + MomentumTransferPhaseSystem::invADs(invADs); + MomentumTransferPhaseSystem::invADs(invADfs); } template -Foam::PtrList -Foam::MomentumTransferPhaseSystem::KdPhifs() const +Foam::PtrList> +Foam::MomentumTransferPhaseSystem::invADfs +( + const PtrList& As +) const { - PtrList KdPhifs(this->phaseModels_.size()); + const label n = As.size(); + + PtrList> invADfs(n); + + forAll(invADfs, i) + { + invADfs.set(i, new PtrList(n)); + invADfs[i].set(i, As[i].clone()); + } + + labelList movingPhases(this->phases().size(), -1); + + forAll(this->movingPhases(), movingPhasei) + { + movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei; + } - // Add the explicit part of the drag force forAllConstIter(KdfTable, Kdfs_, KdfIter) { const surfaceScalarField& Kf(*KdfIter()); @@ -889,111 +1001,55 @@ Foam::MomentumTransferPhaseSystem::KdPhifs() const const phaseModel& phase = iter(); const phaseModel& otherPhase = iter.otherPhase(); - const surfaceScalarField alphaf - ( - fvc::interpolate(otherPhase) - ); + const label i = movingPhases[phase.index()]; - addField - ( - phase, - "KdPhif", - -(alphaf/max(alphaf, otherPhase.residualAlpha())) - *Kf - *fvc::absolute + if (i != -1) + { + const surfaceScalarField alphaf ( - this->MRF().absolute(otherPhase.phi()), - otherPhase.U() - ), - KdPhifs - ); + fvc::interpolate(otherPhase) + ); + + const surfaceScalarField Kfij + ( + (alphaf/max(alphaf, otherPhase.residualAlpha()))*Kf + ); + + invADfs[i][i] += Kfij; + + const label j = movingPhases[otherPhase.index()]; + + if (j != -1) + { + invADfs[i].set(j, -Kfij); + } + } } } - if (this->fillFields_) + for (label i=0; ifillFields - ( - "KdPhif", - dimArea*dimDensity*dimAcceleration, - KdPhifs - ); - } - - return KdPhifs; -} - - -template -Foam::PtrList -Foam::MomentumTransferPhaseSystem::Kds() const -{ - PtrList Kds(this->phaseModels_.size()); - - forAllConstIter(KdTable, Kds_, KdIter) - { - const volScalarField& K(*KdIter()); - const phaseInterface interface(*this, KdIter.key()); - - forAllConstIter(phaseInterface, interface, iter) + for (label j=0; jmesh(), + dimensionedScalar(As[0].dimensions(), 0) + ) + ); + } } } - if (this->fillFields_) - { - this->fillFields("Kd", dimDensity/dimTime, Kds); - } + invADs(invADfs); - return Kds; -} - - -template -Foam::PtrList -Foam::MomentumTransferPhaseSystem::KdUs() const -{ - PtrList KdUs(this->phaseModels_.size()); - - // Add the explicit part of the drag force - forAllConstIter(KdTable, Kds_, KdIter) - { - const volScalarField& K(*KdIter()); - const phaseInterface interface(*this, KdIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - const phaseModel& phase = iter(); - const phaseModel& otherPhase = iter.otherPhase(); - - addField - ( - phase, - "KdU", - -(otherPhase/max(otherPhase, otherPhase.residualAlpha())) - *K*otherPhase.U(), - KdUs - ); - } - } - - if (this->fillFields_) - { - this->fillFields("KdU", dimDensity*dimAcceleration, KdUs); - } - - return KdUs; + return invADfs; } @@ -1032,8 +1088,7 @@ template Foam::tmp Foam::MomentumTransferPhaseSystem::alphaDByAf ( - const PtrList& rAUs, - const PtrList& rAUfs + const PtrList& rAs ) const { tmp alphaDByAf; @@ -1049,15 +1104,7 @@ Foam::MomentumTransferPhaseSystem::alphaDByAf ( alphaDByAf, fvc::interpolate(max(phase, scalar(0))) - *( - rAUfs.size() - - // Face-momentum form - ? rAUfs[phasei]*fvc::interpolate(pPrime()) - - // Cell-momentum form - : fvc::interpolate(rAUs[phasei]*pPrime(), pPrime().name()) - ) + *fvc::interpolate(rAs[phasei]*pPrime(), pPrime().name()) ); } @@ -1087,27 +1134,14 @@ Foam::MomentumTransferPhaseSystem::alphaDByAf alphaDByAf, alpha1f*alpha2f /max(alpha1f + alpha2f, interface.phase1().residualAlpha()) - *( - rAUfs.size() - - // Face-momentum form - ? max + *fvc::interpolate + ( + max ( - rAUfs[interface.phase1().index()], - rAUfs[interface.phase2().index()] - ) - *fvc::interpolate(turbulentDispersionModelIter()->D()) - - // Cell-momentum form - : fvc::interpolate - ( - max - ( - rAUs[interface.phase1().index()], - rAUs[interface.phase2().index()] - ) - *turbulentDispersionModelIter()->D() + rAs[interface.phase1().index()], + rAs[interface.phase2().index()] ) + *turbulentDispersionModelIter()->D() ) ); } @@ -1232,20 +1266,18 @@ void Foam::MomentumTransferPhaseSystem::dragCorrs PtrList& dragCorrfs ) const { - const phaseSystem::phaseModelList& phases = this->phaseModels_; + labelList movingPhases(this->phases().size(), -1); + PtrList Uphis(this->movingPhases().size()); - PtrList Uphis(phases.size()); - - forAll(phases, i) + forAll(this->movingPhases(), movingPhasei) { - if (!phases[i].stationary()) - { - Uphis.set - ( - i, - fvc::reconstruct(phases[i].phi()) - ); - } + movingPhases[this->movingPhases()[movingPhasei].index()] = movingPhasei; + + Uphis.set + ( + movingPhasei, + fvc::reconstruct(this->movingPhases()[movingPhasei].phi()) + ); } forAllConstIter(KdTable, Kds_, KdIter) @@ -1253,384 +1285,40 @@ void Foam::MomentumTransferPhaseSystem::dragCorrs const volScalarField& K(*KdIter()); const phaseInterface interface(*this, KdIter.key()); - const phaseModel& phase1 = interface.phase1(); - const phaseModel& phase2 = interface.phase2(); - - const label phase1i = phase1.index(); - const label phase2i = phase2.index(); - - if (!phase1.stationary()) + forAllConstIter(phaseInterface, interface, iter) { - const volScalarField K1 - ( - (phase2/max(phase2, phase2.residualAlpha()))*K - ); + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); - addField - ( - phase1, - "dragCorr", - K1 - *( - phase2.stationary() - ? -Uphis[phase1i] - : (Uphis[phase2i] - Uphis[phase1i]) - ), - dragCorrs - ); + const label i = movingPhases[phase.index()]; + const label j = movingPhases[otherPhase.index()]; - addField - ( - phase1, - "dragCorrf", - fvc::interpolate(K1) - *( - phase2.stationary() - ? -phase1.phi() - : (phase2.phi() - phase1.phi()) - ), - dragCorrfs - ); - } - - if (!phase2.stationary()) - { - const volScalarField K2 - ( - (phase1/max(phase1, phase1.residualAlpha()))*K - ); - - addField - ( - phase2, - "dragCorr", - K2 - *( - phase1.stationary() - ? -Uphis[phase2i] - : (Uphis[phase1i] - Uphis[phase2i]) - ), - dragCorrs - ); - - addField - ( - phase2, - "dragCorrf", - fvc::interpolate(K2) - *( - phase1.stationary() - ? -phase2.phi() - : (phase1.phi() - phase2.phi()) - ), - dragCorrfs - ); - } - } -} - - -template -void Foam::MomentumTransferPhaseSystem::partialElimination -( - const PtrList& rAUs, - const PtrList& KdUs, - const PtrList& alphafs, - const PtrList& rAUfs, - const PtrList& KdPhis -) -{ - Info<< "Inverting drag systems: "; - - phaseSystem::phaseModelList& phases = this->phaseModels_; - - // Remove the drag contributions from the velocity and flux of the phases - // in preparation for the partial elimination of these terms - forAll(phases, i) - { - if (!phases[i].stationary()) - { - phases[i].URef() += rAUs[i]*KdUs[i]; - phases[i].phiRef() += rAUfs[i]*KdPhis[i]; - } - } - - { - // Create drag coefficient matrices - PtrList> KdByAs(phases.size()); - PtrList> KdByAfs(phases.size()); - - forAll(phases, i) - { - KdByAs.set - ( - i, - new PtrList(phases.size()) - ); - - KdByAfs.set - ( - i, - new PtrList(phases.size()) - ); - } - - forAllConstIter(KdTable, Kds_, KdIter) - { - const volScalarField& K(*KdIter()); - const phaseInterface interface(*this, KdIter.key()); - - const label phase1i = interface.phase1().index(); - const label phase2i = interface.phase2().index(); - - const volScalarField K1 - ( - interface.phase2() - /max(interface.phase2(), interface.phase2().residualAlpha()) - *K - ); - - const volScalarField K2 - ( - interface.phase1() - /max(interface.phase1(), interface.phase1().residualAlpha()) - *K - ); - - addField - ( - interface.phase2(), - "KdByA", - -rAUs[phase1i]*K1, - KdByAs[phase1i] - ); - addField - ( - interface.phase1(), - "KdByA", - -rAUs[phase2i]*K2, - KdByAs[phase2i] - ); - - addField - ( - interface.phase2(), - "KdByAf", - -rAUfs[phase1i]*fvc::interpolate(K1), - KdByAfs[phase1i] - ); - addField - ( - interface.phase1(), - "KdByAf", - -rAUfs[phase2i]*fvc::interpolate(K2), - KdByAfs[phase2i] - ); - } - - forAll(phases, i) - { - this->fillFields("KdByAs", dimless, KdByAs[i]); - this->fillFields("KdByAfs", dimless, KdByAfs[i]); - - KdByAs[i][i] = 1; - KdByAfs[i][i] = 1; - } - - // Decompose - for (label i = 0; i < phases.size(); i++) - { - for (label j = i + 1; j < phases.size(); j++) + if (i != -1) { - KdByAs[j][i] /= KdByAs[i][i]; - KdByAfs[j][i] /= KdByAfs[i][i]; - for (label k = i + 1; k < phases.size(); ++ k) - { - KdByAs[j][k] -= KdByAs[j][i]*KdByAs[i][k]; - KdByAfs[j][k] -= KdByAfs[j][i]*KdByAfs[i][k]; - } - } - } + const volScalarField K1 + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha()))*K + ); - { - volScalarField detKdByAs(KdByAs[0][0]); - surfaceScalarField detPhiKdfs(KdByAfs[0][0]); + addField + ( + i, + "dragCorr", + K1*(j == -1 ? -Uphis[i] : (Uphis[j] - Uphis[i])), + dragCorrs + ); - for (label i = 1; i < phases.size(); i++) - { - detKdByAs *= KdByAs[i][i]; - detPhiKdfs *= KdByAfs[i][i]; - } - - Info<< "Min cell/face det = " << gMin(detKdByAs.primitiveField()) - << "/" << gMin(detPhiKdfs.primitiveField()) << endl; - } - - // Solve for the velocities and fluxes - for (label i = 1; i < phases.size(); i++) - { - if (!phases[i].stationary()) - { - for (label j = 0; j < i; j ++) - { - phases[i].URef() -= KdByAs[i][j]*phases[j].U(); - phases[i].phiRef() -= KdByAfs[i][j]*phases[j].phi(); - } - } - } - for (label i = phases.size() - 1; i >= 0; i--) - { - if (!phases[i].stationary()) - { - for (label j = phases.size() - 1; j > i; j--) - { - phases[i].URef() -= KdByAs[i][j]*phases[j].U(); - phases[i].phiRef() -= KdByAfs[i][j]*phases[j].phi(); - } - phases[i].URef() /= KdByAs[i][i]; - phases[i].phiRef() /= KdByAfs[i][i]; + addField + ( + i, + "dragCorrf", + fvc::interpolate(K1) + *(j == -1 ? -phase.phi() : (otherPhase.phi() - phase.phi())), + dragCorrfs + ); } } } - - // Adjust the phase-fluxes such that the mean flux - // obtained from the pressure solution is retained - this->setMixturePhi(alphafs, this->phi()); -} - - -template -void Foam::MomentumTransferPhaseSystem::partialEliminationf -( - const PtrList& rAUfs, - const PtrList& alphafs, - const PtrList& KdPhifs -) -{ - Info<< "Inverting drag system: "; - - phaseSystem::phaseModelList& phases = this->phaseModels_; - - // Remove the drag contributions from the flux of the phases - // in preparation for the partial elimination of these terms - forAll(phases, i) - { - if (!phases[i].stationary()) - { - phases[i].phiRef() += rAUfs[i]*KdPhifs[i]; - } - } - - { - // Create drag coefficient matrix - PtrList> phiKdfs(phases.size()); - - forAll(phases, phasei) - { - phiKdfs.set - ( - phasei, - new PtrList(phases.size()) - ); - } - - forAllConstIter(KdfTable, Kdfs_, KdfIter) - { - const surfaceScalarField& Kf(*KdfIter()); - const phaseInterface interface(*this, KdfIter.key()); - - const label phase1i = interface.phase1().index(); - const label phase2i = interface.phase2().index(); - - const surfaceScalarField alpha1f - ( - fvc::interpolate(interface.phase1()) - ); - - const surfaceScalarField alpha2f - ( - fvc::interpolate(interface.phase2()) - ); - - addField - ( - interface.phase2(), - "phiKdf", - -rAUfs[phase1i] - *alpha2f/max(alpha2f, interface.phase2().residualAlpha()) - *Kf, - phiKdfs[phase1i] - ); - addField - ( - interface.phase1(), - "phiKdf", - -rAUfs[phase2i] - *alpha1f/max(alpha1f, interface.phase1().residualAlpha()) - *Kf, - phiKdfs[phase2i] - ); - } - - forAll(phases, phasei) - { - this->fillFields("phiKdf", dimless, phiKdfs[phasei]); - - phiKdfs[phasei][phasei] = 1; - } - - // Decompose - for (label i = 0; i < phases.size(); i++) - { - for (label j = i + 1; j < phases.size(); j++) - { - phiKdfs[j][i] /= phiKdfs[i][i]; - for (label k = i + 1; k < phases.size(); ++ k) - { - phiKdfs[j][k] -= phiKdfs[j][i]*phiKdfs[i][k]; - } - } - } - - { - surfaceScalarField detPhiKdfs(phiKdfs[0][0]); - - for (label i = 1; i < phases.size(); i++) - { - detPhiKdfs *= phiKdfs[i][i]; - } - - Info<< "Min face det = " - << gMin(detPhiKdfs.primitiveField()) << endl; - } - - // Solve for the fluxes - for (label i = 1; i < phases.size(); i++) - { - if (!phases[i].stationary()) - { - for (label j = 0; j < i; j ++) - { - phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi(); - } - } - } - for (label i = phases.size() - 1; i >= 0; i--) - { - if (!phases[i].stationary()) - { - for (label j = phases.size() - 1; j > i; j--) - { - phases[i].phiRef() -= phiKdfs[i][j]*phases[j].phi(); - } - phases[i].phiRef() /= phiKdfs[i][i]; - } - } - } - - // Adjust the phase-fluxes such that the mean flux - // obtained from the pressure solution is retained - this->setMixturePhi(alphafs, this->phi()); } diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H index e33b5854d6..8c35852810 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H @@ -169,6 +169,16 @@ protected: const tmp& field ) const; + //- Invert the ADs matrix inplace + void invADs(List>& ADs) const; + + //- Invert the ADs matrix inplace + template class PatchField, class GeoMesh> + void invADs + ( + PtrList>>& ADs + ) const; + public: @@ -194,7 +204,7 @@ public: //- Return implicit force coefficients on the faces, for the face-based // algorithm. - virtual PtrList KdVmfs() const; + virtual PtrList Vmfs() const; //- Return the explicit force fluxes for the cell-based algorithm, that // do not depend on phase mass/volume fluxes, and can therefore be @@ -205,22 +215,19 @@ public: //- As Fs, but for the face-based algorithm virtual PtrList Ffs() const; - //- Return the explicit drag force fluxes for the cell-based algorithm. - // These depend on phase mass/volume fluxes, and must therefore be - // evaluated inside the corrector loop. - virtual PtrList KdPhis() const; + //- Return the inverse of (As + Ds) + virtual void invADs + ( + const PtrList& As, + PtrList>& invADs, + PtrList>& invADfs + ) const; - //- As KdPhis, but for the face-based algorithm - virtual PtrList KdPhifs() const; - - //- Return the implicit part of the drag force - virtual PtrList Kds() const; - - //- Return the explicit part of the drag force for the cell-based - // algorithm. This is the cell-equivalent of KdPhis. These depend on - // phase velocities, and must therefore be evaluated inside the - // corrector loop. - virtual PtrList KdUs() const; + //- Return the inverse of (Afs + Dfs) + virtual PtrList> invADfs + ( + const PtrList& Afs + ) const; //- Returns true if the phase pressure is treated implicitly // in the phase fraction equation @@ -234,8 +241,7 @@ public: // divided by the momentum central coefficient virtual tmp alphaDByAf ( - const PtrList& rAUs, - const PtrList& rAUfs + const PtrList& rAs ) const; //- Return the flux corrections for the cell-based algorithm. These @@ -250,25 +256,6 @@ public: PtrList& dragCorrf ) const; - //- Solve the drag system for the velocities and fluxes - virtual void partialElimination - ( - const PtrList& rAUs, - const PtrList& KdUs, - const PtrList& alphafs, - const PtrList& rAUfs, - const PtrList& KdPhis - ); - - //- As partialElimination, but for the face-based algorithm. Only solves - // for the fluxes. - virtual void partialEliminationf - ( - const PtrList& rAUfs, - const PtrList& alphafs, - const PtrList& KdPhifs - ); - //- Read base phaseProperties dictionary virtual bool read(); }; diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C index 29aa0bd9e9..1700346241 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C @@ -175,11 +175,10 @@ Foam::PopulationBalancePhaseSystem::specieTransfer() const template void Foam::PopulationBalancePhaseSystem::solve ( - const PtrList& rAUs, - const PtrList& rAUfs + const PtrList& rAs ) { - BasePhaseSystem::solve(rAUs, rAUfs); + BasePhaseSystem::solve(rAs); forAll(populationBalances_, i) { diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H index c550e6d5b5..c18bd65659 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H @@ -97,11 +97,7 @@ public: specieTransfer() const; //- Solve all population balance equations - virtual void solve - ( - const PtrList& rAUs, - const PtrList& rAUfs - ); + virtual void solve(const PtrList& rAs); //- Correct derived properties virtual void correct(); diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H index dc3ac0bf39..4a87aaa0a4 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H @@ -523,7 +523,7 @@ public: //- Return the implicit force coefficients for the face-based // algorithm - virtual PtrList KdVmfs() const = 0; + virtual PtrList Vmfs() const = 0; //- Return the force fluxes for the cell-based algorithm virtual PtrList Fs() const = 0; @@ -531,17 +531,19 @@ public: //- Return the force fluxes for the face-based algorithm virtual PtrList Ffs() const = 0; - //- Return the force fluxes for the cell-based algorithm - virtual PtrList KdPhis() const = 0; + //- Return the inverse of (As + Ds) + virtual void invADs + ( + const PtrList& As, + PtrList>& invADs, + PtrList>& invADfs + ) const = 0; - //- Return the force fluxes for the face-based algorithm - virtual PtrList KdPhifs() const = 0; - - //- Return the implicit part of the drag force - virtual PtrList Kds() const = 0; - - //- Return the explicit part of the drag force - virtual PtrList KdUs() const = 0; + //- Return the inverse of (Afs + Dfs) + virtual PtrList> invADfs + ( + const PtrList& Afs + ) const = 0; //- Returns true if the phase pressure is treated implicitly // in the phase fraction equation @@ -555,8 +557,7 @@ public: // divided by the momentum central coefficient virtual tmp alphaDByAf ( - const PtrList& rAUs, - const PtrList& rAUfs + const PtrList& rAs ) const = 0; //- Return the flux corrections for the cell-based algorithm @@ -569,24 +570,6 @@ public: PtrList& dragCorrf ) const = 0; - //- Solve the drag system for the new velocities and fluxes - virtual void partialElimination - ( - const PtrList& rAUs, - const PtrList& KdUs, - const PtrList& alphafs, - const PtrList& rAUfs, - const PtrList& KdPhis - ) = 0; - - //- Solve the drag system for the new fluxes - virtual void partialEliminationf - ( - const PtrList& rAUfs, - const PtrList& alphafs, - const PtrList& KdPhifs - ) = 0; - //- Re-normalise the flux of the phases // around the specified mixture mean void setMixturePhi @@ -611,11 +594,7 @@ public: // Evolution //- Solve for the phase fractions - virtual void solve - ( - const PtrList& rAUs, - const PtrList& rAUfs - ); + virtual void solve(const PtrList& rAs); //- Correct the fluid properties other than those listed below virtual void correct(); diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H index bb8966b07e..2b986358da 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H @@ -23,6 +23,67 @@ License \*---------------------------------------------------------------------------*/ +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +class wordListAndType +{ + public: + + wordList wl; + Type t; + + wordListAndType() + {} + + wordListAndType(Istream& is) + : + wl(is), + t(is) + {} +}; + +template +inline Istream& operator>>(Istream& is, wordListAndType& wlat) +{ + return is >> wlat.wl >> wlat.t; +} + +template +inline Ostream& operator<<(Ostream& os, const wordListAndType& wlat) +{ + return os << wlat.wl << wlat.t; +} + +template +inline bool operator== +( + const wordListAndType& a, + const wordListAndType& b +) +{ + return a.wl == b.wl && a.t == b.t; +} + +template +inline bool operator!= +( + const wordListAndType& a, + const wordListAndType& b +) +{ + return !(a == b); +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // inline const Foam::fvMesh& Foam::phaseSystem::mesh() const @@ -190,4 +251,151 @@ inline const Foam::dimensionedScalar& Foam::phaseSystem::deltaN() const } +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +inline void addField +( + const label phasei, + const word& name, + tmp field, + PtrList& fieldList +) +{ + if (fieldList.set(phasei)) + { + fieldList[phasei] += field; + } + else + { + fieldList.set + ( + phasei, + new GeoField + ( + name, // IOobject::groupName(name, group.name()), + field + ) + ); + } +} + + +template +inline void addField +( + const Group& group, + const word& name, + tmp field, + PtrList& fieldList +) +{ + if (fieldList.set(group.index())) + { + fieldList[group.index()] += field; + } + else + { + fieldList.set + ( + group.index(), + new GeoField + ( + IOobject::groupName(name, group.name()), + field + ) + ); + } +} + + +template +inline void addField +( + const Group& group, + const word& name, + const GeoField& field, + PtrList& fieldList +) +{ + addField(group, name, tmp(field), fieldList); +} + + +template +inline void addField +( + const Group& group, + const word& name, + tmp field, + HashPtrTable& fieldTable +) +{ + if (fieldTable.found(group.name())) + { + *fieldTable[group.name()] += field; + } + else + { + fieldTable.set + ( + group.name(), + new GeoField + ( + IOobject::groupName(name, group.name()), + field + ) + ); + } +} + + +template +inline void addField +( + const Group& group, + const word& name, + const GeoField& field, + HashPtrTable& fieldTable +) +{ + addField(group, name, tmp(field), fieldTable); +} + + +template class PatchField, class GeoMesh> +PtrList> operator& +( + const PtrList>>& As, + const PtrList>& fs +) +{ + PtrList> Afs(fs.size()); + + forAll(Afs, i) + { + Afs.set(i, As[i][i]*fs[i]); + + forAll(Afs, j) + { + if (j != i) + { + Afs[i] += As[i][j]*fs[j]; + } + } + } + + return Afs; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + // ************************************************************************* // diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C index 54826f13a7..c266eaa0e8 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C @@ -41,11 +41,7 @@ License // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -void Foam::phaseSystem::solve -( - const PtrList& rAUs, - const PtrList& rAUfs -) +void Foam::phaseSystem::solve(const PtrList& rAs) { const dictionary& alphaControls = mesh_.solution().solverDict("alpha"); @@ -251,9 +247,9 @@ void Foam::phaseSystem::solve PtrList alphaPhis(phases().size()); tmp alphaDByAf; - if (implicitPhasePressure() && (rAUs.size() || rAUfs.size())) + if (implicitPhasePressure() && (rAs.size())) { - alphaDByAf = this->alphaDByAf(rAUs, rAUfs); + alphaDByAf = this->alphaDByAf(rAs); } forAll(movingPhases(), movingPhasei) @@ -454,7 +450,7 @@ void Foam::phaseSystem::solve if (alphaDByAf.valid()) { // Update alphaDByAf due to changes in alpha - alphaDByAf = this->alphaDByAf(rAUs, rAUfs); + alphaDByAf = this->alphaDByAf(rAs); forAll(solvePhases, solvePhasei) { diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C index a89722ad23..b08a727f58 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C @@ -28,62 +28,6 @@ License // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // -namespace Foam -{ - -template -class wordListAndType -{ - public: - - wordList wl; - Type t; - - wordListAndType() - {} - - wordListAndType(Istream& is) - : - wl(is), - t(is) - {} -}; - -template -inline Istream& operator>>(Istream& is, wordListAndType& wlat) -{ - return is >> wlat.wl >> wlat.t; -} - -template -inline Ostream& operator<<(Ostream& os, const wordListAndType& wlat) -{ - return os << wlat.wl << wlat.t; -} - -template -inline bool operator== -( - const wordListAndType& a, - const wordListAndType& b -) -{ - return a.wl == b.wl && a.t == b.t; -} - -template -inline bool operator!= -( - const wordListAndType& a, - const wordListAndType& b -) -{ - return !(a == b); -} - -} - - template Foam::dictionary Foam::phaseSystem::interfacialDict(const word& name) const { @@ -362,7 +306,6 @@ void Foam::phaseSystem::generateInterfacialModels } - template void Foam::phaseSystem::generateInterfacialModels ( @@ -483,96 +426,4 @@ const ModelType& Foam::phaseSystem::lookupInterfacialModel } -// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -template -inline void addField -( - const Group& group, - const word& name, - tmp field, - PtrList& fieldList -) -{ - if (fieldList.set(group.index())) - { - fieldList[group.index()] += field; - } - else - { - fieldList.set - ( - group.index(), - new GeoField - ( - IOobject::groupName(name, group.name()), - field - ) - ); - } -} - - -template -inline void addField -( - const Group& group, - const word& name, - const GeoField& field, - PtrList& fieldList -) -{ - addField(group, name, tmp(field), fieldList); -} - - -template -inline void addField -( - const Group& group, - const word& name, - tmp field, - HashPtrTable& fieldTable -) -{ - if (fieldTable.found(group.name())) - { - *fieldTable[group.name()] += field; - } - else - { - fieldTable.set - ( - group.name(), - new GeoField - ( - IOobject::groupName(name, group.name()), - field - ) - ); - } -} - - -template -inline void addField -( - const Group& group, - const word& name, - const GeoField& field, - HashPtrTable& fieldTable -) -{ - addField(group, name, tmp(field), fieldTable); -} - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H index 48be7b3912..9b5e54e999 100644 --- a/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H +++ b/src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -90,6 +90,11 @@ public: return word("SymmetricSquareMatrix<") + pTraits::typeName + '>'; } + + + // Member Operators + + using Matrix, Type>::operator=; }; diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleSingleRegionControl/pimpleSingleRegionControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleSingleRegionControl/pimpleSingleRegionControl.C index 2e28cf6907..ed80807ac8 100644 --- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleSingleRegionControl/pimpleSingleRegionControl.C +++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleSingleRegionControl/pimpleSingleRegionControl.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -45,7 +45,7 @@ Foam::pimpleSingleRegionControl::pimpleSingleRegionControl { pimple_.pimpleLoopPtr_ = this; - read(); + pimpleLoop::read(); pimple_.printResidualControls(); diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C index 77a6fb7952..a84393f121 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C @@ -457,7 +457,7 @@ EulerDdtScheme::fvcDdtUfCorr return fluxFieldType::New ( "ddtCorr(" + U.name() + ',' + Uf.name() + ')', - this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr) *rDeltaT*phiCorr + this->fvcDdtPhiCoeff(U.oldTime(), phiUf0, phiCorr)*rDeltaT*phiCorr ); } diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H index 19d044783a..a6e6f0ac53 100644 --- a/src/finiteVolume/fvMesh/fvMesh.H +++ b/src/finiteVolume/fvMesh/fvMesh.H @@ -536,7 +536,7 @@ public: //- Return the fvSchemes const fvSchemes& schemes() const; - //- Return the fvSchemes + //- Return the fvSolution const fvSolution& solution() const; diff --git a/tutorials/multiRegion/CHT/wallBoiling/system/fluid/fvSolution b/tutorials/multiRegion/CHT/wallBoiling/system/fluid/fvSolution index 6eb378643d..ad6294a529 100644 --- a/tutorials/multiRegion/CHT/wallBoiling/system/fluid/fvSolution +++ b/tutorials/multiRegion/CHT/wallBoiling/system/fluid/fvSolution @@ -68,7 +68,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/Grossetete/system/fvSolution b/tutorials/multiphaseEuler/Grossetete/system/fvSolution index d4490ada20..207db178e2 100644 --- a/tutorials/multiphaseEuler/Grossetete/system/fvSolution +++ b/tutorials/multiphaseEuler/Grossetete/system/fvSolution @@ -77,7 +77,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/LBend/system/fvSolution b/tutorials/multiphaseEuler/LBend/system/fvSolution index 74e86e10d0..ef58bdcc9c 100644 --- a/tutorials/multiphaseEuler/LBend/system/fvSolution +++ b/tutorials/multiphaseEuler/LBend/system/fvSolution @@ -92,7 +92,6 @@ PIMPLE faceMomentum yes; VmDdtCorrection no; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bed/system/fvSolution b/tutorials/multiphaseEuler/bed/system/fvSolution index 9523801f62..59ceabec26 100644 --- a/tutorials/multiphaseEuler/bed/system/fvSolution +++ b/tutorials/multiphaseEuler/bed/system/fvSolution @@ -57,10 +57,11 @@ solvers PIMPLE { - nOuterCorrectors 3; - nCorrectors 1; + nOuterCorrectors 3; + nCorrectors 1; nNonOrthogonalCorrectors 0; + dragCorrection yes; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumn/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumn/system/fvSolution index c218156b43..4b726bdc2f 100644 --- a/tutorials/multiphaseEuler/bubbleColumn/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumn/system/fvSolution @@ -73,7 +73,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumnEvaporating/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumnEvaporating/system/fvSolution index f7ac899b18..d62ca7b3eb 100644 --- a/tutorials/multiphaseEuler/bubbleColumnEvaporating/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumnEvaporating/system/fvSolution @@ -73,7 +73,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumnEvaporatingDissolving/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumnEvaporatingDissolving/system/fvSolution index b8b0759221..6485458ab1 100644 --- a/tutorials/multiphaseEuler/bubbleColumnEvaporatingDissolving/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumnEvaporatingDissolving/system/fvSolution @@ -70,7 +70,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumnEvaporatingReacting/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumnEvaporatingReacting/system/fvSolution index 73e00700e1..c9643cc241 100644 --- a/tutorials/multiphaseEuler/bubbleColumnEvaporatingReacting/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumnEvaporatingReacting/system/fvSolution @@ -82,7 +82,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumnIATE/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumnIATE/system/fvSolution index 8a18a29ee7..18370ab628 100644 --- a/tutorials/multiphaseEuler/bubbleColumnIATE/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumnIATE/system/fvSolution @@ -63,7 +63,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumnLES/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumnLES/system/fvSolution index cc87db0b46..329c570631 100644 --- a/tutorials/multiphaseEuler/bubbleColumnLES/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumnLES/system/fvSolution @@ -73,7 +73,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubbleColumnLaminar/system/fvSolution b/tutorials/multiphaseEuler/bubbleColumnLaminar/system/fvSolution index 3c829cb38d..2ffe12cdb5 100644 --- a/tutorials/multiphaseEuler/bubbleColumnLaminar/system/fvSolution +++ b/tutorials/multiphaseEuler/bubbleColumnLaminar/system/fvSolution @@ -64,7 +64,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/bubblePipe/system/fvSolution b/tutorials/multiphaseEuler/bubblePipe/system/fvSolution index ac21994f9d..ee645ce655 100644 --- a/tutorials/multiphaseEuler/bubblePipe/system/fvSolution +++ b/tutorials/multiphaseEuler/bubblePipe/system/fvSolution @@ -73,7 +73,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/damBreak4phase/system/fvSolution b/tutorials/multiphaseEuler/damBreak4phase/system/fvSolution index fa0023dd27..735e804492 100644 --- a/tutorials/multiphaseEuler/damBreak4phase/system/fvSolution +++ b/tutorials/multiphaseEuler/damBreak4phase/system/fvSolution @@ -79,8 +79,7 @@ PIMPLE faceMomentum no; VmDdtCorrection no; - dragCorrection no; - partialElimination no; + dragCorrection yes; } relaxationFactors diff --git a/tutorials/multiphaseEuler/fluidisedBed/system/fvSolution b/tutorials/multiphaseEuler/fluidisedBed/system/fvSolution index 13f8b985d3..824ac18808 100644 --- a/tutorials/multiphaseEuler/fluidisedBed/system/fvSolution +++ b/tutorials/multiphaseEuler/fluidisedBed/system/fvSolution @@ -90,7 +90,6 @@ PIMPLE faceMomentum no; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/fluidisedBedLaminar/system/fvSolution b/tutorials/multiphaseEuler/fluidisedBedLaminar/system/fvSolution index 4b39a15a9e..6e70342d11 100644 --- a/tutorials/multiphaseEuler/fluidisedBedLaminar/system/fvSolution +++ b/tutorials/multiphaseEuler/fluidisedBedLaminar/system/fvSolution @@ -89,8 +89,7 @@ PIMPLE nNonOrthogonalCorrectors 0; faceMomentum no; - dragCorrection yes; - partialElimination no; + dragCorrection no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/hydrofoil/system/fvSolution b/tutorials/multiphaseEuler/hydrofoil/system/fvSolution index 56d65e01ff..94878ffe33 100644 --- a/tutorials/multiphaseEuler/hydrofoil/system/fvSolution +++ b/tutorials/multiphaseEuler/hydrofoil/system/fvSolution @@ -83,7 +83,6 @@ PIMPLE faceMomentum no; VmDdtCorrection no; dragCorrection yes; - partialElimination no; } cache diff --git a/tutorials/multiphaseEuler/injection/system/fvSolution b/tutorials/multiphaseEuler/injection/system/fvSolution index 3c829cb38d..2ffe12cdb5 100644 --- a/tutorials/multiphaseEuler/injection/system/fvSolution +++ b/tutorials/multiphaseEuler/injection/system/fvSolution @@ -64,7 +64,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/mixerVessel2D/0/alpha.map b/tutorials/multiphaseEuler/mixerVessel2D/0/alpha.map index d43b824229..e5be25d94b 100644 --- a/tutorials/multiphaseEuler/mixerVessel2D/0/alpha.map +++ b/tutorials/multiphaseEuler/mixerVessel2D/0/alpha.map @@ -14,7 +14,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions []; +dimensions [0 0 0 0 0 0 0]; internalField uniform 0.6; diff --git a/tutorials/multiphaseEuler/mixerVessel2D/system/fvSolution b/tutorials/multiphaseEuler/mixerVessel2D/system/fvSolution index 7112f44639..91decb1250 100644 --- a/tutorials/multiphaseEuler/mixerVessel2D/system/fvSolution +++ b/tutorials/multiphaseEuler/mixerVessel2D/system/fvSolution @@ -98,7 +98,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/mixerVessel2DMRF/system/fvSolution b/tutorials/multiphaseEuler/mixerVessel2DMRF/system/fvSolution index fe4010e9ce..e2a78fbe87 100644 --- a/tutorials/multiphaseEuler/mixerVessel2DMRF/system/fvSolution +++ b/tutorials/multiphaseEuler/mixerVessel2DMRF/system/fvSolution @@ -85,7 +85,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/pipeBend/system/fvSolution b/tutorials/multiphaseEuler/pipeBend/system/fvSolution index 1fc5fc3869..e5d4be527e 100644 --- a/tutorials/multiphaseEuler/pipeBend/system/fvSolution +++ b/tutorials/multiphaseEuler/pipeBend/system/fvSolution @@ -89,7 +89,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/steamInjection/system/fvSolution b/tutorials/multiphaseEuler/steamInjection/system/fvSolution index f076db9589..0a9e82af4e 100644 --- a/tutorials/multiphaseEuler/steamInjection/system/fvSolution +++ b/tutorials/multiphaseEuler/steamInjection/system/fvSolution @@ -86,7 +86,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination yes; } relaxationFactors diff --git a/tutorials/multiphaseEuler/titaniaSynthesis/system/fvSolution b/tutorials/multiphaseEuler/titaniaSynthesis/system/fvSolution index be55f8e5de..3d80ae9683 100644 --- a/tutorials/multiphaseEuler/titaniaSynthesis/system/fvSolution +++ b/tutorials/multiphaseEuler/titaniaSynthesis/system/fvSolution @@ -86,7 +86,6 @@ PIMPLE faceMomentum no; VmDdtCorrection no; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/titaniaSynthesisSurface/system/fvSolution b/tutorials/multiphaseEuler/titaniaSynthesisSurface/system/fvSolution index be55f8e5de..3d80ae9683 100644 --- a/tutorials/multiphaseEuler/titaniaSynthesisSurface/system/fvSolution +++ b/tutorials/multiphaseEuler/titaniaSynthesisSurface/system/fvSolution @@ -86,7 +86,6 @@ PIMPLE faceMomentum no; VmDdtCorrection no; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/wallBoilingIATE/system/fvSolution b/tutorials/multiphaseEuler/wallBoilingIATE/system/fvSolution index 5aa01332e3..7ec84f46aa 100644 --- a/tutorials/multiphaseEuler/wallBoilingIATE/system/fvSolution +++ b/tutorials/multiphaseEuler/wallBoilingIATE/system/fvSolution @@ -68,6 +68,7 @@ solvers PIMPLE { + momentumPredictor no; nOuterCorrectors 1; nCorrectors 1; nNonOrthogonalCorrectors 0; @@ -76,19 +77,19 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors { fields { - thermalPhaseChange:dmdtf 1.0; + thermalPhaseChange:dmdtf 1; } equations { ".*" 1; + "U\..*" 1; "h\..*" 1; } } diff --git a/tutorials/multiphaseEuler/wallBoilingPolydisperse/system/fvSolution b/tutorials/multiphaseEuler/wallBoilingPolydisperse/system/fvSolution index ca7c37f7b8..c0747a5bf0 100644 --- a/tutorials/multiphaseEuler/wallBoilingPolydisperse/system/fvSolution +++ b/tutorials/multiphaseEuler/wallBoilingPolydisperse/system/fvSolution @@ -85,7 +85,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors diff --git a/tutorials/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/fvSolution b/tutorials/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/fvSolution index ca7c37f7b8..c0747a5bf0 100644 --- a/tutorials/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/fvSolution +++ b/tutorials/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/fvSolution @@ -85,7 +85,6 @@ PIMPLE faceMomentum no; VmDdtCorrection yes; dragCorrection yes; - partialElimination no; } relaxationFactors