From 5fd30443f30aa9835527f4644c9c147a82b1c8b6 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Sat, 26 Aug 2023 10:09:38 +0100 Subject: [PATCH] multiphaseEuler: New implicit drag algorithm replaces partial-elimination corrector The momentum equation central coefficient and drag matrix is formulated, inverted and used to eliminate the drag terms from each of the phase momentum equations which are combined for formulate a drag-implicit pressure equation. This eliminates the lagged drag terms from the previous formulation which significantly improves convergence for small particle and Euler-VoF high-drag cases. It would also be possible to refactor the virtual-mass terms and include the central coefficients of the phase acceleration terms in the drag matrix before inversion to further improve the implicitness of the phase momentum-pressure coupling for bubbly flows. This work is pending funding. --- .../multiphaseEuler/cellPressureCorrector.C | 419 +++++---- .../multiphaseEuler/facePressureCorrector.C | 291 +++--- .../multiphaseEuler/momentumPredictor.C | 3 - .../multiphaseEuler/multiphaseEuler.C | 18 +- .../multiphaseEuler/multiphaseEuler.H | 14 +- .../MomentumTransferPhaseSystem.C | 846 ++++++------------ .../MomentumTransferPhaseSystem.H | 61 +- .../PopulationBalancePhaseSystem.C | 5 +- .../PopulationBalancePhaseSystem.H | 6 +- .../phaseSystems/phaseSystem/phaseSystem.H | 51 +- .../phaseSystems/phaseSystem/phaseSystemI.H | 208 +++++ .../phaseSystem/phaseSystemSolve.C | 12 +- .../phaseSystem/phaseSystemTemplates.C | 149 --- .../SymmetricSquareMatrix.H | 7 +- .../pimpleSingleRegionControl.C | 4 +- .../EulerDdtScheme/EulerDdtScheme.C | 2 +- src/finiteVolume/fvMesh/fvMesh.H | 2 +- .../CHT/wallBoiling/system/fluid/fvSolution | 1 - .../Grossetete/system/fvSolution | 1 - .../multiphaseEuler/LBend/system/fvSolution | 1 - .../multiphaseEuler/bed/system/fvSolution | 5 +- .../bubbleColumn/system/fvSolution | 1 - .../bubbleColumnEvaporating/system/fvSolution | 1 - .../system/fvSolution | 1 - .../system/fvSolution | 1 - .../bubbleColumnIATE/system/fvSolution | 1 - .../bubbleColumnLES/system/fvSolution | 1 - .../bubbleColumnLaminar/system/fvSolution | 1 - .../bubblePipe/system/fvSolution | 1 - .../damBreak4phase/system/fvSolution | 3 +- .../fluidisedBed/system/fvSolution | 1 - .../fluidisedBedLaminar/system/fvSolution | 3 +- .../hydrofoil/system/fvSolution | 1 - .../injection/system/fvSolution | 1 - .../multiphaseEuler/mixerVessel2D/0/alpha.map | 2 +- .../mixerVessel2D/system/fvSolution | 1 - .../mixerVessel2DMRF/system/fvSolution | 1 - .../pipeBend/system/fvSolution | 1 - .../steamInjection/system/fvSolution | 1 - .../titaniaSynthesis/system/fvSolution | 1 - .../titaniaSynthesisSurface/system/fvSolution | 1 - .../wallBoilingIATE/system/fvSolution | 5 +- .../wallBoilingPolydisperse/system/fvSolution | 1 - .../system/fvSolution | 1 - 44 files changed, 926 insertions(+), 1212 deletions(-) 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