From 3c18db77ede67cb4d83e89ceea951f693685fc9b Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Sun, 3 Sep 2023 10:21:31 +0100 Subject: [PATCH] multiphaseEuler::StationaryPhaseModel: Replaced functions returning 0 with notImplemented requiring all parts of the moving phase solution algorithm to loop and operate on the moving phases only making the code easier to understand and maintain. --- .../multiphaseEuler/cellPressureCorrector.C | 5 - .../multiphaseEuler/facePressureCorrector.C | 5 - .../multiphaseEuler/multiphaseEuler.C | 7 +- .../multiphaseEuler/multiphaseEuler.H | 5 + .../MomentumTransferPhaseSystem.C | 9 +- .../StationaryPhaseModel.C | 119 +++++++++--------- .../phaseSystems/phaseSystem/phaseSystem.C | 47 +++---- .../phaseSystems/phaseSystem/phaseSystem.H | 44 ------- .../phaseSystem/phaseSystemSolve.C | 4 +- .../phaseSystem/phaseSystemTemplates.C | 72 ----------- 10 files changed, 94 insertions(+), 223 deletions(-) diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C index 8b84be088d..c264ebc8a7 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C @@ -45,11 +45,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() { volScalarField& p(p_); - const phaseSystem::phaseModelPartialList& movingPhases - ( - fluid.movingPhases() - ); - // Face volume fractions PtrList alphafs(phases.size()); forAll(phases, phasei) diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C index 68c4ca1634..fcade0cefb 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C @@ -45,11 +45,6 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() { volScalarField& p(p_); - const phaseSystem::phaseModelPartialList& movingPhases - ( - fluid.movingPhases() - ); - // Face volume fractions PtrList alphafs(phases.size()); forAll(phases, phasei) diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C index baa70283cd..af4df8acff 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C +++ b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C @@ -72,12 +72,12 @@ void Foam::solvers::multiphaseEuler::correctCoNum() fvc::surfaceSum(mag(phi))().primitiveField() ); - forAll(phases, phasei) + forAll(movingPhases, movingPhasei) { sumPhi = max ( sumPhi, - fvc::surfaceSum(mag(phases[phasei].phi()))().primitiveField() + fvc::surfaceSum(mag(phases[movingPhasei].phi()))().primitiveField() ); } @@ -164,6 +164,8 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh) phases_(fluid_.phases()), + movingPhases_(fluid_.movingPhases()), + phi_(fluid_.phi()), p_(phases_[0].thermo().p()), @@ -182,6 +184,7 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh) fluid(fluid_), phases(phases_), + movingPhases(movingPhases_), p(p_), phi(phi_) { diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H index 34bf2c36ac..22b5d75264 100644 --- a/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H +++ b/applications/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H @@ -113,6 +113,8 @@ protected: phaseSystem::phaseModelList& phases_; + phaseSystem::phaseModelPartialList& movingPhases_; + surfaceScalarField& phi_; @@ -202,6 +204,9 @@ public: //- Reference to the phases const phaseSystem::phaseModelList& phases; + //- Reference to the moving phases + const phaseSystem::phaseModelPartialList& movingPhases; + //- Reference to the pressure field const volScalarField& p; diff --git a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index f982c16fdc..0e92dac3c4 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -253,10 +253,9 @@ Foam::MomentumTransferPhaseSystem::Fs() const } // Add the phase pressure - forAll(this->phaseModels_, phasei) + forAll(this->movingPhases(), movingPhasei) { - const phaseModel& phase = this->phaseModels_[phasei]; - + const phaseModel& phase = this->movingPhases()[movingPhasei]; const tmp pPrime(phase.pPrime()); addField @@ -375,9 +374,9 @@ Foam::MomentumTransferPhaseSystem::Ffs() const } // Add the phase pressure - forAll(this->phaseModels_, phasei) + forAll(this->movingPhases(), movingPhasei) { - const phaseModel& phase = this->phaseModels_[phasei]; + const phaseModel& phase = this->movingPhases()[movingPhasei]; addField ( diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C index fddd2054aa..5d84e7ff53 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C @@ -63,7 +63,7 @@ Foam::StationaryPhaseModel::UEqn() { FatalErrorInFunction << "Cannot construct a momentum equation for a stationary phase" - << exit(FatalError); + << abort(FatalError); return tmp(); } @@ -75,7 +75,7 @@ Foam::StationaryPhaseModel::UfEqn() { FatalErrorInFunction << "Cannot construct a momentum equation for a stationary phase" - << exit(FatalError); + << abort(FatalError); return tmp(); } @@ -85,12 +85,11 @@ template Foam::tmp Foam::StationaryPhaseModel::U() const { - return volVectorField::New - ( - IOobject::groupName("U", this->name()), - this->mesh(), - dimensionedVector(dimVelocity, Zero) - ); + FatalErrorInFunction + << "Cannot access the velocity of a stationary phase" + << abort(FatalError); + + return volVectorField::null(); } @@ -100,7 +99,7 @@ Foam::StationaryPhaseModel::URef() { FatalErrorInFunction << "Cannot access the velocity of a stationary phase" - << exit(FatalError); + << abort(FatalError); return const_cast(volVectorField::null()); } @@ -112,7 +111,7 @@ Foam::StationaryPhaseModel::URef() const { FatalErrorInFunction << "Cannot access the velocity of a stationary phase" - << exit(FatalError); + << abort(FatalError); return volVectorField::null(); } @@ -122,12 +121,11 @@ template Foam::tmp Foam::StationaryPhaseModel::phi() const { - return surfaceScalarField::New - ( - IOobject::groupName("phi", this->name()), - this->mesh(), - dimensionedScalar(dimVolume/dimTime, 0) - ); + FatalErrorInFunction + << "Cannot access the flux of a stationary phase" + << abort(FatalError); + + return surfaceScalarField::null(); } @@ -137,7 +135,7 @@ Foam::StationaryPhaseModel::phiRef() { FatalErrorInFunction << "Cannot access the flux of a stationary phase" - << exit(FatalError); + << abort(FatalError); return const_cast(surfaceScalarField::null()); } @@ -149,7 +147,7 @@ Foam::StationaryPhaseModel::phiRef() const { FatalErrorInFunction << "Cannot access the flux of a stationary phase" - << exit(FatalError); + << abort(FatalError); return surfaceScalarField::null(); } @@ -174,7 +172,7 @@ Foam::StationaryPhaseModel::UfRef() { FatalErrorInFunction << "Cannot access the face velocity of a stationary phase" - << exit(FatalError); + << abort(FatalError); return const_cast(surfaceVectorField::null()); } @@ -186,7 +184,7 @@ Foam::StationaryPhaseModel::UfRef() const { FatalErrorInFunction << "Cannot access the face velocity of a stationary phase" - << exit(FatalError); + << abort(FatalError); return surfaceVectorField::null(); } @@ -196,12 +194,11 @@ template Foam::tmp Foam::StationaryPhaseModel::alphaPhi() const { - return surfaceScalarField::New - ( - IOobject::groupName("alphaPhi", this->name()), - this->mesh(), - dimensionedScalar(dimVolume/dimTime, 0) - ); + FatalErrorInFunction + << "Cannot access the flux of a stationary phase" + << abort(FatalError); + + return surfaceScalarField::null(); } @@ -211,7 +208,7 @@ Foam::StationaryPhaseModel::alphaPhiRef() { FatalErrorInFunction << "Cannot access the volumetric flux of a stationary phase" - << exit(FatalError); + << abort(FatalError); return const_cast(surfaceScalarField::null()); } @@ -222,8 +219,8 @@ const Foam::surfaceScalarField& Foam::StationaryPhaseModel::alphaPhiRef() const { FatalErrorInFunction - << "Cannot access the volumetric flux of a stationary phase" - << exit(FatalError); + << "Cannot access the flux of a stationary phase" + << abort(FatalError); return surfaceScalarField::null(); } @@ -233,12 +230,11 @@ template Foam::tmp Foam::StationaryPhaseModel::alphaRhoPhi() const { - return surfaceScalarField::New - ( - IOobject::groupName("alphaRhoPhi", this->name()), - this->mesh(), - dimensionedScalar(dimMass/dimTime, 0) - ); + FatalErrorInFunction + << "Cannot access the flux of a stationary phase" + << abort(FatalError); + + return surfaceScalarField::null(); } @@ -247,8 +243,8 @@ Foam::surfaceScalarField& Foam::StationaryPhaseModel::alphaRhoPhiRef() { FatalErrorInFunction - << "Cannot access the mass flux of a stationary phase" - << exit(FatalError); + << "Cannot access the flux of a stationary phase" + << abort(FatalError); return const_cast(surfaceScalarField::null()); } @@ -259,8 +255,8 @@ const Foam::surfaceScalarField& Foam::StationaryPhaseModel::alphaRhoPhiRef() const { FatalErrorInFunction - << "Cannot access the mass flux of a stationary phase" - << exit(FatalError); + << "Cannot access the flux of a stationary phase" + << abort(FatalError); return surfaceScalarField::null(); } @@ -271,8 +267,8 @@ Foam::tmp Foam::StationaryPhaseModel::UgradU() const { FatalErrorInFunction - << "Cannot calculate DUDt of a stationary phase" - << exit(FatalError); + << "Cannot calculate UgradU of a stationary phase" + << abort(FatalError); return tmp(nullptr); } @@ -284,7 +280,7 @@ Foam::StationaryPhaseModel::DUDt() const { FatalErrorInFunction << "Cannot calculate DUDt of a stationary phase" - << exit(FatalError); + << abort(FatalError); return tmp(nullptr); } @@ -294,12 +290,11 @@ template Foam::tmp Foam::StationaryPhaseModel::continuityError() const { - return volScalarField::New - ( - IOobject::groupName("continuityError", this->name()), - this->mesh(), - dimensionedScalar(dimDensity/dimTime, 0) - ); + FatalErrorInFunction + << "Cannot access the continuityError of a stationary phase" + << abort(FatalError); + + return volScalarField::null(); } @@ -307,12 +302,11 @@ template Foam::tmp Foam::StationaryPhaseModel::K() const { - return volScalarField::New - ( - IOobject::groupName("K", this->name()), - this->mesh(), - dimensionedScalar(sqr(dimVelocity), 0) - ); + FatalErrorInFunction + << "Cannot access the kinetic energy of a stationary phase" + << abort(FatalError); + + return volScalarField::null(); } @@ -320,6 +314,10 @@ template const Foam::autoPtr& Foam::StationaryPhaseModel::divU() const { + FatalErrorInFunction + << "Cannot access the dilatation rate of a stationary phase" + << abort(FatalError); + static autoPtr divU_; return divU_; } @@ -333,7 +331,7 @@ void Foam::StationaryPhaseModel::divU { FatalErrorInFunction << "Cannot set the dilatation rate of a stationary phase" - << exit(FatalError); + << abort(FatalError); } @@ -362,12 +360,11 @@ template Foam::tmp Foam::StationaryPhaseModel::pPrime() const { - return volScalarField::New - ( - IOobject::groupName("pPrime", this->name()), - this->mesh(), - dimensionedScalar(dimPressure, 0) - ); + FatalErrorInFunction + << "Cannot access the pPrime of a stationary phase" + << abort(FatalError); + + return volScalarField::null(); } diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C index 49bee4cc04..700fda5595 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C @@ -34,7 +34,7 @@ License #include "fvCorrectPhi.H" #include "fvcMeshPhi.H" #include "correctContactAngle.H" -#include "dragModel.H" +#include "fixedValueFvsPatchFields.H" #include "movingWallVelocityFvPatchVectorField.H" #include "pressureReference.H" @@ -51,30 +51,6 @@ const Foam::word Foam::phaseSystem::propertiesName("phaseProperties"); // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // -Foam::tmp Foam::phaseSystem::calcPhi -( - const phaseModelList& phaseModels -) const -{ - tmp tmpPhi - ( - surfaceScalarField::New - ( - "phi", - fvc::interpolate(phaseModels[0])*phaseModels[0].phi() - ) - ); - - for (label phasei=1; phasei Foam::phaseSystem::sumAlphaMoving() const { tmp sumAlphaMoving @@ -237,7 +213,17 @@ Foam::phaseSystem::phaseSystem phaseModel::iNew(*this, referencePhaseName_) ), - phi_("phi", calcPhi(phaseModels_)), + phi_ + ( + IOobject + ( + "phi", + mesh.time().name(), + mesh + ), + mesh, + dimensionedScalar(dimFlux, 0) + ), dpdt_ ( @@ -300,6 +286,13 @@ Foam::phaseSystem::phaseSystem } } + forAll(movingPhaseModels_, movingPhasei) + { + phi_ += + fvc::interpolate(movingPhaseModels_[movingPhasei]) + *movingPhaseModels_[movingPhasei].phi(); + } + // Write phi phi_.writeOpt() = IOobject::AUTO_WRITE; @@ -332,7 +325,7 @@ Foam::phaseSystem::phaseSystem } } - forAll(phases(), phasei) + forAll(phaseModels_, phasei) { const volScalarField& alphai = phases()[phasei]; mesh_.schemes().setFluxRequired(alphai.name()); diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H index 31953f14c3..04e875b4f8 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H @@ -172,20 +172,8 @@ protected: interfaceSurfaceTensionModelTable interfaceSurfaceTensionModels_; - //- Flag to indicate that returned lists of fields are "complete"; i.e., - // that an absence of force is returned as a constructed list of zeros, - // rather than a null pointer - static const bool fillFields_ = false; - - // Protected member functions - //- Calculate and return the mixture flux - tmp calcPhi - ( - const phaseModelList& phaseModels - ) const; - //- Return the sum of the phase fractions of the moving phases tmp sumAlphaMoving() const; @@ -436,38 +424,6 @@ public: ) const; - // Field construction - - //- Fill up gaps in a phase-indexed list of fields with zeros - template - < - class Type, - template class PatchField, - class GeoMesh - > - void fillFields - ( - const word& name, - const dimensionSet& dims, - PtrList>& fieldList - ) const; - - //- Fill up gaps in a phase-indexed table of fields with zeros - template - < - class Type, - template class PatchField, - class GeoMesh - > - void fillFields - ( - const word& name, - const dimensionSet& dims, - HashPtrTable>& - fieldTable - ) const; - - // Properties //- Return the mixture density diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C index c266eaa0e8..2484c499d9 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemSolve.C @@ -182,9 +182,9 @@ void Foam::phaseSystem::solve(const PtrList& rAs) ) ); - forAll(phases(), phasej) + forAll(movingPhases(), movingPhasej) { - const phaseModel& phase2 = phases()[phasej]; + const phaseModel& phase2 = movingPhases()[movingPhasej]; const volScalarField& alpha2 = phase2; if (&phase2 != &phase) diff --git a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C index b08a727f58..f5ede1d10c 100644 --- a/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C +++ b/applications/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemTemplates.C @@ -97,78 +97,6 @@ Foam::dictionary Foam::phaseSystem::interfacialDict(const word& name) const // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -template class PatchField, class GeoMesh> -void Foam::phaseSystem::fillFields -( - const word& name, - const dimensionSet& dims, - PtrList>& fieldList -) const -{ - forAll(this->phaseModels_, phasei) - { - if (fieldList.set(phasei)) - { - continue; - } - - const phaseModel& phase = this->phaseModels_[phasei]; - - fieldList.set - ( - phasei, - new GeometricField - ( - IOobject - ( - IOobject::groupName(name, phase.name()), - this->mesh_.time().name(), - this->mesh_ - ), - this->mesh_, - dimensioned("zero", dims, pTraits::zero) - ) - ); - } -} - - -template class PatchField, class GeoMesh> -void Foam::phaseSystem::fillFields -( - const word& name, - const dimensionSet& dims, - HashPtrTable>& fieldTable -) const -{ - forAll(this->phaseModels_, phasei) - { - const phaseModel& phase = this->phaseModels_[phasei]; - - if (fieldTable.set(phase.name())) - { - continue; - } - - fieldTable.set - ( - phase.name(), - new GeometricField - ( - IOobject - ( - IOobject::groupName(name, phase.name()), - this->mesh_.time().name(), - this->mesh_ - ), - this->mesh_, - dimensioned("zero", dims, pTraits::zero) - ) - ); - } -} - - template Foam::word Foam::phaseSystem::modelName() const {