From e98dcc5aa8449f79d6dfe236ebf23123387d89ec Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Fri, 24 Mar 2023 17:23:14 +0000 Subject: [PATCH] solvers: Added ddtCorr support in MRF regions by extending the use of Uf and rhoUf to provide the old-time absolute flux. This avoids possible pressure-velocity-flux decoupling (staggering) within the MRF region using ddtCorr to better couple the velocity and flux fields. --- .../solvers/modules/VoFSolver/VoFSolver.C | 5 +- .../pressureCorrector.C | 2 +- .../compressibleVoF/pressureCorrector.C | 2 +- .../incompressibleFluid/correctPressure.C | 2 +- .../incompressibleFluid/incompressibleFluid.C | 5 +- .../pressureCorrector.C | 2 +- .../incompressibleVoF/pressureCorrector.C | 2 +- .../isothermalFluid/correctBuoyantPressure.C | 16 +++-- .../modules/isothermalFluid/correctPressure.C | 16 +++-- .../modules/isothermalFluid/isothermalFluid.C | 8 ++- .../MomentumTransferPhaseSystem.C | 13 ++-- .../MovingPhaseModel/MovingPhaseModel.C | 19 +++--- .../MovingPhaseModel/MovingPhaseModel.H | 6 +- .../StationaryPhaseModel.C | 7 +- .../StationaryPhaseModel.H | 4 +- .../phaseModel/phaseModel/phaseModel.H | 4 +- .../phaseSystems/phaseSystem/phaseSystem.C | 6 +- .../phaseSystems/phaseSystem/phaseSystem.H | 8 +-- .../multiphase/driftFluxFoam/createFields.H | 25 ++++++++ .../solvers/multiphase/driftFluxFoam/pEqn.H | 5 +- .../cfdTools/general/MRF/IOMRFZoneList.H | 44 ++++++++++++- .../cfdTools/general/MRF/MRFZoneList.C | 34 +++++++++- .../cfdTools/general/MRF/MRFZoneList.H | 19 ++---- .../general/MRF/MRFZoneListTemplates.C | 64 ------------------- .../EulerDdtScheme/EulerDdtScheme.C | 28 ++++---- src/finiteVolume/finiteVolume/fvc/fvcDdt.C | 9 ++- src/finiteVolume/finiteVolume/fvc/fvcDdt.H | 2 +- .../finiteVolume/fvc/fvcMeshPhiTemplates.C | 12 ++-- 28 files changed, 203 insertions(+), 166 deletions(-) delete mode 100644 src/finiteVolume/cfdTools/general/MRF/MRFZoneListTemplates.C diff --git a/applications/solvers/modules/VoFSolver/VoFSolver.C b/applications/solvers/modules/VoFSolver/VoFSolver.C index 38c73e1bcb..9769c90aa7 100644 --- a/applications/solvers/modules/VoFSolver/VoFSolver.C +++ b/applications/solvers/modules/VoFSolver/VoFSolver.C @@ -121,10 +121,13 @@ Foam::solvers::VoFSolver::VoFSolver { mesh.schemes().setFluxRequired(p_rgh.name()); - if (mesh.dynamic()) + if (mesh.dynamic() || MRF.size()) { Info<< "Constructing face momentum Uf" << endl; + // Ensure the U BCs are up-to-date before constructing Uf + U_.correctBoundaryConditions(); + Uf = new surfaceVectorField ( IOobject diff --git a/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C b/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C index 46d530c18e..cbc9005f2c 100644 --- a/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C +++ b/applications/solvers/modules/compressibleMultiphaseVoF/pressureCorrector.C @@ -63,7 +63,7 @@ void Foam::solvers::compressibleMultiphaseVoF::pressureCorrector() ( "phiHbyA", fvc::flux(HbyA) - + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)) + + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf) ); MRF.makeRelative(phiHbyA); diff --git a/applications/solvers/modules/compressibleVoF/pressureCorrector.C b/applications/solvers/modules/compressibleVoF/pressureCorrector.C index b39a547835..ba7e95b407 100644 --- a/applications/solvers/modules/compressibleVoF/pressureCorrector.C +++ b/applications/solvers/modules/compressibleVoF/pressureCorrector.C @@ -70,7 +70,7 @@ void Foam::solvers::compressibleVoF::pressureCorrector() ( "phiHbyA", fvc::flux(HbyA) - + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)) + + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf) ); MRF.makeRelative(phiHbyA); diff --git a/applications/solvers/modules/incompressibleFluid/correctPressure.C b/applications/solvers/modules/incompressibleFluid/correctPressure.C index 5edb3ce8e9..d667c31572 100644 --- a/applications/solvers/modules/incompressibleFluid/correctPressure.C +++ b/applications/solvers/modules/incompressibleFluid/correctPressure.C @@ -46,7 +46,7 @@ void Foam::solvers::incompressibleFluid::correctPressure() ( "phiHbyA", fvc::flux(HbyA) - + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf)) + + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf) ); MRF.makeRelative(phiHbyA); diff --git a/applications/solvers/modules/incompressibleFluid/incompressibleFluid.C b/applications/solvers/modules/incompressibleFluid/incompressibleFluid.C index b2d49dbeee..38673f4851 100644 --- a/applications/solvers/modules/incompressibleFluid/incompressibleFluid.C +++ b/applications/solvers/modules/incompressibleFluid/incompressibleFluid.C @@ -119,10 +119,13 @@ Foam::solvers::incompressibleFluid::incompressibleFluid(fvMesh& mesh) momentumTransport->validate(); - if (mesh.dynamic()) + if (mesh.dynamic() || MRF.size()) { Info<< "Constructing face momentum Uf" << endl; + // Ensure the U BCs are up-to-date before constructing Uf + U.correctBoundaryConditions(); + Uf = new surfaceVectorField ( IOobject diff --git a/applications/solvers/modules/incompressibleMultiphaseVoF/pressureCorrector.C b/applications/solvers/modules/incompressibleMultiphaseVoF/pressureCorrector.C index ee75c7f840..e2b9d4937d 100644 --- a/applications/solvers/modules/incompressibleMultiphaseVoF/pressureCorrector.C +++ b/applications/solvers/modules/incompressibleMultiphaseVoF/pressureCorrector.C @@ -62,7 +62,7 @@ void Foam::solvers::incompressibleMultiphaseVoF::pressureCorrector() ( "phiHbyA", fvc::flux(HbyA) - + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)) + + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf) ); MRF.makeRelative(phiHbyA); diff --git a/applications/solvers/modules/incompressibleVoF/pressureCorrector.C b/applications/solvers/modules/incompressibleVoF/pressureCorrector.C index ea7ceb7073..7058a13345 100644 --- a/applications/solvers/modules/incompressibleVoF/pressureCorrector.C +++ b/applications/solvers/modules/incompressibleVoF/pressureCorrector.C @@ -62,7 +62,7 @@ void Foam::solvers::incompressibleVoF::pressureCorrector() ( "phiHbyA", fvc::flux(HbyA) - + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf)) + + fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf) ); MRF.makeRelative(phiHbyA); diff --git a/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C b/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C index 6c60aa0e75..03b9c1dc7a 100644 --- a/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C +++ b/applications/solvers/modules/isothermalFluid/correctBuoyantPressure.C @@ -55,6 +55,8 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() // pressure solution const volScalarField psip0(psi*p); + const surfaceScalarField rhof(fvc::interpolate(rho)); + const volScalarField rAU("rAU", 1.0/UEqn.A()); const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); @@ -86,11 +88,11 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() surfaceScalarField phiHbyA ( "phiHbyA", - fvc::interpolate(rho)*fvc::flux(HbyA) - + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf)) + rhof*fvc::flux(HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf) ); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + MRF.makeRelative(rhof, phiHbyA); const bool adjustMass = mesh.schemes().steady() && adjustPhi(phiHbyA, U, p_rgh); @@ -107,7 +109,7 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() ( constrainPhid ( - fvc::relative(phiHbyA, rho, U)/fvc::interpolate(rho), + fvc::relative(phiHbyA, rho, U)/rhof, p_rgh ) ); @@ -259,14 +261,14 @@ void Foam::solvers::isothermalFluid::correctBuoyantPressure() rho = thermo.rho(); } - // Correct rhoUf if the mesh is moving - fvc::correctRhoUf(rhoUf, rho, U, phi, MRF); - if (mesh.schemes().steady() || pimple.simpleRho()) { rho.relax(); } + // Correct rhoUf if the mesh is moving + fvc::correctRhoUf(rhoUf, rho, U, phi, MRF); + if (thermo.dpdt()) { dpdt = fvc::ddt(p); diff --git a/applications/solvers/modules/isothermalFluid/correctPressure.C b/applications/solvers/modules/isothermalFluid/correctPressure.C index e5c24e0a6b..4e648eb16b 100644 --- a/applications/solvers/modules/isothermalFluid/correctPressure.C +++ b/applications/solvers/modules/isothermalFluid/correctPressure.C @@ -51,6 +51,8 @@ void Foam::solvers::isothermalFluid::correctPressure() // pressure solution const volScalarField psip0(psi*p); + const surfaceScalarField rhof(fvc::interpolate(rho)); + const volScalarField rAU("rAU", 1.0/UEqn.A()); const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); @@ -82,11 +84,11 @@ void Foam::solvers::isothermalFluid::correctPressure() surfaceScalarField phiHbyA ( "phiHbyA", - fvc::interpolate(rho)*fvc::flux(HbyA) - + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf)) + rhof*fvc::flux(HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf) ); - MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + MRF.makeRelative(rhof, phiHbyA); bool adjustMass = false; @@ -96,7 +98,7 @@ void Foam::solvers::isothermalFluid::correctPressure() ( constrainPhid ( - fvc::relative(phiHbyA, rho, U)/fvc::interpolate(rho), + fvc::relative(phiHbyA, rho, U)/rhof, p ) ); @@ -236,14 +238,14 @@ void Foam::solvers::isothermalFluid::correctPressure() rho = thermo.rho(); } - // Correct rhoUf if the mesh is moving - fvc::correctRhoUf(rhoUf, rho, U, phi, MRF); - if (mesh.schemes().steady() || pimple.simpleRho()) { rho.relax(); } + // Correct rhoUf if the mesh is moving + fvc::correctRhoUf(rhoUf, rho, U, phi, MRF); + if (thermo.dpdt()) { dpdt = fvc::ddt(p); diff --git a/applications/solvers/modules/isothermalFluid/isothermalFluid.C b/applications/solvers/modules/isothermalFluid/isothermalFluid.C index 78c23b8ac6..dcd921a99e 100644 --- a/applications/solvers/modules/isothermalFluid/isothermalFluid.C +++ b/applications/solvers/modules/isothermalFluid/isothermalFluid.C @@ -217,10 +217,13 @@ Foam::solvers::isothermalFluid::isothermalFluid ); } - if (mesh.dynamic()) + if (mesh.dynamic() || MRF.size()) { Info<< "Constructing face momentum rhoUf" << endl; + // Ensure the U BCs are up-to-date before constructing Uf + U.correctBoundaryConditions(); + rhoUf = new surfaceVectorField ( IOobject @@ -362,6 +365,9 @@ void Foam::solvers::isothermalFluid::postSolve() if (!mesh.schemes().steady()) { rho = thermo.rho(); + + // Correct rhoUf with the updated density if the mesh is moving + fvc::correctRhoUf(rhoUf, rho, U, phi, MRF); } } diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index 6434229be9..98f14ceaec 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -1067,15 +1067,12 @@ Foam::MomentumTransferPhaseSystem::ddtCorrByAs phiCorrs.set ( phasei, - this->MRF().zeroFilter ( - ( - phase.Uf().valid() - ? (this->mesh_.Sf() & phase.Uf()().oldTime())() - : phase.phi()().oldTime() - ) - - fvc::flux(phase.U()().oldTime()) - )() + phase.Uf().valid() + ? (this->mesh_.Sf() & phase.Uf()().oldTime())() + : phase.phi()().oldTime() + ) + - fvc::flux(phase.U()().oldTime()) ); } diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C index 97097b2034..d866bdedde 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -196,7 +196,7 @@ Foam::MovingPhaseModel::MovingPhaseModel { phi_.writeOpt() = IOobject::AUTO_WRITE; - if (fluid.mesh().dynamic()) + if (fluid.mesh().dynamic() || this->fluid().MRF().size()) { Uf_ = new surfaceVectorField ( @@ -305,11 +305,11 @@ void Foam::MovingPhaseModel::correctUf() { const fvMesh& mesh = this->fluid().mesh(); - if (mesh.dynamic()) + if (Uf_.valid()) { - Uf_.ref() = fvc::interpolate(U_); + Uf_() = fvc::interpolate(U_); surfaceVectorField n(mesh.Sf()/mesh.magSf()); - Uf_.ref() += + Uf_() += n*( this->fluid().MRF().absolute(fvc::absolute(phi_, U_)) /mesh.magSf() @@ -396,13 +396,10 @@ Foam::MovingPhaseModel::phiRef() template -Foam::tmp +const Foam::autoPtr& Foam::MovingPhaseModel::Uf() const { - return - Uf_.valid() - ? tmp(Uf_()) - : tmp(); + return Uf_; } @@ -412,7 +409,7 @@ Foam::MovingPhaseModel::UfRef() { if (Uf_.valid()) { - return Uf_.ref(); + return Uf_(); } else { diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H index 076e106c80..25d363b812 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -110,7 +110,7 @@ protected: surfaceScalarField alphaRhoPhi_; //- Face velocity field - tmp Uf_; + autoPtr Uf_; //- Lagrangian acceleration field (needed for virtual-mass) mutable tmp DUDt_; @@ -219,7 +219,7 @@ public: //- Return the face velocity // Required for moving mesh cases - virtual tmp Uf() const; + virtual const autoPtr& Uf() const; //- Access the face velocity // Required for moving mesh cases diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C index fbef4be716..6c81294ca4 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -132,14 +132,15 @@ Foam::StationaryPhaseModel::phiRef() template -Foam::tmp +const Foam::autoPtr& Foam::StationaryPhaseModel::Uf() const { FatalErrorInFunction << "Cannot access the face velocity of a stationary phase" << abort(FatalError); - return tmp(); + static autoPtr Uf_; + return Uf_; } diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H index e49eb0b5eb..46f1f750ae 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -102,7 +102,7 @@ public: //- Return the face velocity // Required for moving mesh cases - virtual tmp Uf() const; + virtual const autoPtr& Uf() const; //- Access the face velocity // Required for moving mesh cases diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H index 9d9838695b..cf7c6a6b0c 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/phaseModel/phaseModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -322,7 +322,7 @@ public: //- Return the face velocity // Required for moving mesh cases - virtual tmp Uf() const = 0; + virtual const autoPtr& Uf() const = 0; //- Access the face velocity // Required for moving mesh cases diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C index fce607a0a2..39e2859fc2 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C @@ -225,6 +225,8 @@ Foam::phaseSystem::phaseSystem mesh_(mesh), + MRF_(mesh_), + referencePhaseName_(lookupOrDefault("referencePhase", word::null)), phaseModels_ @@ -247,8 +249,6 @@ Foam::phaseSystem::phaseSystem dimensionedScalar(dimPressure/dimTime, 0) ), - MRF_(mesh_), - deltaN_ ( "deltaN", @@ -773,7 +773,7 @@ void Foam::phaseSystem::correctPhi // Calculate absolute flux // from the mapped surface velocity - phi_ += alphafs[phasei]*(mesh_.Sf() & phase.Uf()); + phi_ += alphafs[phasei]*(mesh_.Sf() & phase.UfRef()); } if (correctPhi) diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H index bc4b9bf3d7..035c52050d 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -126,6 +126,9 @@ protected: //- Reference to the mesh const fvMesh& mesh_; + //- Optional MRF zones + IOMRFZoneList MRF_; + //- Name of optional reference phase which is not solved for // but obtained from the sum of the other phases word referencePhaseName_; @@ -151,9 +154,6 @@ protected: //- Rate of change of pressure volScalarField dpdt_; - //- Optional MRF zones - IOMRFZoneList MRF_; - //- Interface compression coefficients cAlphaTable cAlphas_; diff --git a/applications/solvers/multiphase/driftFluxFoam/createFields.H b/applications/solvers/multiphase/driftFluxFoam/createFields.H index 44627327b2..4573153cdf 100644 --- a/applications/solvers/multiphase/driftFluxFoam/createFields.H +++ b/applications/solvers/multiphase/driftFluxFoam/createFields.H @@ -97,3 +97,28 @@ tmp talphaPhiCorr0; #include "createFvModels.H" #include "createFvConstraints.H" + +//- Pointer to the surface momentum field +// used for ddtCorr with MRF +autoPtr Uf; + +if (mixture.MRF().size()) +{ + Info<< "Constructing face momentum Uf" << endl; + + // Ensure the U BCs are up-to-date before constructing Uf + U.correctBoundaryConditions(); + + Uf = new surfaceVectorField + ( + IOobject + ( + "Uf", + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + fvc::interpolate(U) + ); +} diff --git a/applications/solvers/multiphase/driftFluxFoam/pEqn.H b/applications/solvers/multiphase/driftFluxFoam/pEqn.H index 24c244cb94..f75940186e 100644 --- a/applications/solvers/multiphase/driftFluxFoam/pEqn.H +++ b/applications/solvers/multiphase/driftFluxFoam/pEqn.H @@ -6,7 +6,7 @@ ( "phiHbyA", fvc::flux(HbyA) - + mixture.MRF().zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi, Uf) ); mixture.MRF().makeRelative(phiHbyA); adjustPhi(phiHbyA, U, p_rgh); @@ -65,4 +65,7 @@ ); p_rgh = p - rho*gh; } + + // Correct Uf for MRF + fvc::correctUf(Uf, U, phi, mixture.MRF()); } diff --git a/src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.H b/src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.H index 59efd0856d..43df79c848 100644 --- a/src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.H +++ b/src/finiteVolume/cfdTools/general/MRF/IOMRFZoneList.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2012-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -94,6 +94,9 @@ public: //- Read dictionary virtual bool read(); + //- Use MRFZoneList::size for looping + using MRFZoneList::size; + // Member Operators @@ -146,6 +149,45 @@ public: { return U; } + + + //- Return the argument unchanged + template + inline const Type& relative(const Type& U, const RhoType& rho) const + { + return U; + } + + //- Return the argument unchanged + template + inline const Type& relative + ( + const Type& U, + const RhoType& rho, + const label patchi + ) const + { + return U; + } + + //- Return the argument unchanged + template + inline const Type& absolute(const Type& U, const RhoType& rho) const + { + return U; + } + + //- Return the argument unchanged + template + inline const Type& absolute + ( + const Type& U, + const RhoType& rho, + const label patchi + ) const + { + return U; + } }; diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C index f69230ecc4..d4b3a2d22c 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2012-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -331,6 +331,38 @@ void Foam::MRFZoneList::makeAbsolute } +Foam::tmp Foam::MRFZoneList::absolute +( + const tmp& tphi, + const volScalarField& rho +) const +{ + if (size()) + { + tmp rphi + ( + New + ( + tphi, + "absolute(" + tphi().name() + ')', + tphi().dimensions(), + true + ) + ); + + makeAbsolute(fvc::interpolate(rho), rphi.ref()); + + tphi.clear(); + + return rphi; + } + else + { + return tmp(tphi, true); + } +} + + void Foam::MRFZoneList::update() { if (mesh_.topoChanged()) diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H index 7dee8f4dc7..c3ebc942ed 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneList.H @@ -153,16 +153,15 @@ public: surfaceScalarField& phi ) const; - void makeAbsolute(Field& Up, const label patchi) const; - - //- Filter-out the MRF region contribution from the given field - // setting the corresponding values to zero - template - tmp> zeroFilter + //- Return the given relative flux absolute within the MRF region + tmp absolute ( - const tmp>& tphi + const tmp& phi, + const volScalarField& rho ) const; + void makeAbsolute(Field& Up, const label patchi) const; + //- Update MRFZone faces if the mesh topology changes void update(); @@ -186,12 +185,6 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -#ifdef NoRepository - #include "MRFZoneListTemplates.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - #endif // ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneListTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneListTemplates.C deleted file mode 100644 index 962466c84d..0000000000 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneListTemplates.C +++ /dev/null @@ -1,64 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2018-2022 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "MRFZoneList.H" - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -Foam::tmp> -Foam::MRFZoneList::zeroFilter -( - const tmp>& tphi -) const -{ - if (size()) - { - tmp zphi - ( - New - ( - tphi, - "zeroFilter(" + tphi().name() + ')', - tphi().dimensions(), - true - ) - ); - - forAll(*this, i) - { - operator[](i).zero(zphi.ref()); - } - - return zphi; - } - else - { - return tmp(tphi, true); - } -} - - -// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C index e6f725e1d0..791200e1e4 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -492,7 +492,7 @@ EulerDdtScheme::fvcDdtUfCorr ( const volScalarField& rho, const VolField& U, - const SurfaceField& Uf + const SurfaceField& rhoUf ) { const dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); @@ -500,31 +500,27 @@ EulerDdtScheme::fvcDdtUfCorr if ( U.dimensions() == dimVelocity - && Uf.dimensions() == rho.dimensions()*dimVelocity + && rhoUf.dimensions() == rho.dimensions()*dimVelocity ) { - VolField rhoU0 - ( - rho.oldTime()*U.oldTime() - ); - - fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime()); + VolField rhoU0(rho.oldTime()*U.oldTime()); + fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime()); fluxFieldType phiCorr(phiUf0 - fvc::dotInterpolate(mesh().Sf(), rhoU0)); return fluxFieldType::New ( - "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', - this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime()) - *rDeltaT*phiCorr + "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')', + this->fvcDdtPhiCoeff(rhoU0, phiUf0, phiCorr, rho.oldTime()) + *rDeltaT*phiCorr ); } else if ( U.dimensions() == rho.dimensions()*dimVelocity - && Uf.dimensions() == rho.dimensions()*dimVelocity + && rhoUf.dimensions() == rho.dimensions()*dimVelocity ) { - fluxFieldType phiUf0(mesh().Sf() & Uf.oldTime()); + fluxFieldType phiUf0(mesh().Sf() & rhoUf.oldTime()); fluxFieldType phiCorr ( phiUf0 - fvc::dotInterpolate(mesh().Sf(), U.oldTime()) @@ -532,7 +528,7 @@ EulerDdtScheme::fvcDdtUfCorr return fluxFieldType::New ( - "ddtCorr(" + rho.name() + ',' + U.name() + ',' + Uf.name() + ')', + "ddtCorr(" + rho.name() + ',' + U.name() + ',' + rhoUf.name() + ')', this->fvcDdtPhiCoeff ( U.oldTime(), @@ -545,7 +541,7 @@ EulerDdtScheme::fvcDdtUfCorr else { FatalErrorInFunction - << "dimensions of Uf are not correct" + << "dimensions of rhoUf are not correct" << abort(FatalError); return fluxFieldType::null(); diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C index 01e15a4675..b5bc966b36 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C +++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C @@ -229,7 +229,7 @@ tmp::type>> ddtCorr const autoPtr>& Uf ) { - if (U.mesh().dynamic()) + if (Uf.valid()) { return ddtCorr(U, Uf()); } @@ -271,19 +271,18 @@ tmp::type>> ddtCorr ).ref().fvcDdtPhiCorr(rho, U, phi); } - template tmp::type>> ddtCorr ( const volScalarField& rho, const VolField& U, const SurfaceField::type>& phi, - const autoPtr>& Uf + const autoPtr>& rhoUf ) { - if (U.mesh().dynamic()) + if (rhoUf.valid()) { - return ddtCorr(rho, U, Uf()); + return ddtCorr(rho, U, rhoUf()); } else { diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.H b/src/finiteVolume/finiteVolume/fvc/fvcDdt.H index 9e525d7021..64d6455018 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.H +++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhiTemplates.C b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhiTemplates.C index c81e7f5766..637920b0e2 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcMeshPhiTemplates.C +++ b/src/finiteVolume/finiteVolume/fvc/fvcMeshPhiTemplates.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 @@ -34,10 +34,10 @@ void Foam::fvc::correctUf { const fvMesh& mesh = U.mesh(); - if (mesh.dynamic()) + if (Uf.valid()) { Uf() = fvc::interpolate(U); - surfaceVectorField n(mesh.Sf()/mesh.magSf()); + const surfaceVectorField n(mesh.Sf()/mesh.magSf()); Uf() += n*(MRF.absolute(phi)/mesh.magSf() - (n & Uf())); } } @@ -55,13 +55,13 @@ void Foam::fvc::correctRhoUf { const fvMesh& mesh = U.mesh(); - if (mesh.dynamic()) + if (rhoUf.valid()) { rhoUf() = fvc::interpolate(rho*U); - surfaceVectorField n(mesh.Sf()/mesh.magSf()); + const surfaceVectorField n(mesh.Sf()/mesh.magSf()); rhoUf() += n*( - MRF.absolute(fvc::absolute(phi, rho, U))/mesh.magSf() + MRF.absolute(fvc::absolute(phi, rho, U), rho)/mesh.magSf() - (n & rhoUf()) ); }