From 9f6c161e494668655df73346de8cbe30477bd870 Mon Sep 17 00:00:00 2001 From: Henry Weller <http://cfd.direct> Date: Thu, 25 Jun 2015 16:10:07 +0100 Subject: [PATCH] reactingTwoPhaseEulerFoam: Generalize the handling of the dilatation rate to support any number of phases --- .../reactingTwoPhaseEulerFoam/pU/pEqn.H | 15 ++-- .../reactingTwoPhaseEulerFoam/pUf/pEqn.H | 14 ++-- .../AnisothermalPhaseModel.C | 34 +++++++++ .../AnisothermalPhaseModel.H | 15 ++++ .../phaseModel/phaseModel/phaseModel.C | 20 +++++ .../phaseModel/phaseModel/phaseModel.H | 13 ++++ .../phaseSystems/phaseSystem/phaseSystem.C | 14 ---- .../phaseSystems/phaseSystem/phaseSystem.H | 9 --- .../phaseSystems/phaseSystem/phaseSystemI.H | 12 --- .../twoPhaseSystem/twoPhaseSystem.C | 74 ++++++++++++++----- 10 files changed, 155 insertions(+), 65 deletions(-) diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H index 30ce7b221f3..9a4174e8c13 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pU/pEqn.H @@ -76,7 +76,7 @@ tmp<surfaceScalarField> phiF2; + (fvc::interpolate(rAU1*F) & mesh.Sf()) ); - // Phase-1 dispersion, lift and wall-lubrication flux + // Phase-2 dispersion, lift and wall-lubrication flux phiF2 = ( - Df2*snGradAlpha1 @@ -339,11 +339,14 @@ while (pimple.correct()) } // Compressibility correction for phase-fraction equations - fluid.dgdt() = - ( - alpha1*(pEqnComp2 & p_rgh) - - alpha2*(pEqnComp1 & p_rgh) - ); + if (phase1.compressible()) + { + phase1.D(pEqnComp1 & p_rgh); + } + if (phase2.compressible()) + { + phase2.D(pEqnComp2 & p_rgh); + } // Optionally relax pressure for velocity correction p_rgh.relax(); diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H index c48f0b07c7e..05ed9e7d087 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/pUf/pEqn.H @@ -323,11 +323,15 @@ while (pimple.correct()) U2.correctBoundaryConditions(); fvOptions.correct(U2); - fluid.dgdt() = - ( - alpha1*(pEqnComp2 & p_rgh) - - alpha2*(pEqnComp1 & p_rgh) - ); + // Compressibility correction for phase-fraction equations + if (phase1.compressible()) + { + phase1.D(pEqnComp1 & p_rgh); + } + if (phase2.compressible()) + { + phase2.D(pEqnComp2 & p_rgh); + } } } diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C index d62af01b539..99e16af8a05 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C @@ -36,6 +36,17 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel ) : BasePhaseModel(fluid, phaseName), + D_ + ( + IOobject + ( + IOobject::groupName("D", this->name()), + fluid.mesh().time().timeName(), + fluid.mesh() + ), + fluid.mesh(), + dimensionedScalar("D", dimless/dimTime, 0) + ), K_ ( IOobject @@ -120,4 +131,27 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn() } +template<class BasePhaseModel> +bool Foam::AnisothermalPhaseModel<BasePhaseModel>::compressible() const +{ + return true; +} + + +template<class BasePhaseModel> +Foam::tmp<Foam::volScalarField> +Foam::AnisothermalPhaseModel<BasePhaseModel>::D() const +{ + return D_; +} + + +template<class BasePhaseModel> +void +Foam::AnisothermalPhaseModel<BasePhaseModel>::D(const volScalarField& D) +{ + D_ = D; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H index a9c39535437..c69156b356d 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.H @@ -54,6 +54,9 @@ class AnisothermalPhaseModel { // Private data + //- Dilatation + volScalarField D_; + //- Kinetic energy volScalarField K_; @@ -79,6 +82,18 @@ public: //- Return the enthalpy equation virtual tmp<fvScalarMatrix> heEqn(); + + + // Compressibility (variable density) + + //- Return true if the phase is compressible otherwise false + virtual bool compressible() const; + + //- Phase dilatation rate ((alpha/rho)*Drho/Dt) + virtual tmp<volScalarField> D() const; + + //- Set phase dilatation rate ((alpha/rho)*Drho/Dt) + virtual void D(const volScalarField& D); }; diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C index 35186501904..d1cc94e6c19 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.C @@ -134,4 +134,24 @@ bool Foam::phaseModel::read() } +bool Foam::phaseModel::compressible() const +{ + return false; +} + + +Foam::tmp<Foam::volScalarField> Foam::phaseModel::D() const +{ + return tmp<volScalarField>(); +} + + +void Foam::phaseModel::D(const volScalarField& D) +{ + WarningIn("phaseModel::D(const volScalarField& D)") + << "Attempt to set the dilatation rate of an incompressible phase" + << endl; +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H index 0afb9f78550..e42a3a70b44 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H @@ -83,6 +83,7 @@ public: //- Runtime type information ClassName("phaseModel"); + // Declare runtime construction declareRunTimeSelectionTable @@ -163,6 +164,18 @@ public: virtual bool read(); + // Compressibility (variable density) + + //- Return true if the phase is compressible otherwise false + virtual bool compressible() const; + + //- Phase dilatation rate ((alpha/rho)*Drho/Dt) + virtual tmp<volScalarField> D() const; + + //- Set phase dilatation rate ((alpha/rho)*Drho/Dt) + virtual void D(const volScalarField& D); + + // Thermo //- Return const access to the thermophysical model diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C index 9144093db5a..c89d0eec3e7 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C @@ -183,20 +183,6 @@ Foam::phaseSystem::phaseSystem ), mesh, dimensionedScalar("dpdt", dimPressure/dimTime, 0) - ), - - dgdt_ - ( - IOobject - ( - "dgdt", - mesh.time().timeName(), - mesh, - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE - ), - mesh, - dimensionedScalar("dgdt", dimless/dimTime, 0) ) { // Blending methods diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H index d7d21c430de..46b5ec20e0f 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H @@ -174,9 +174,6 @@ protected: //- Rate of change of pressure volScalarField dpdt_; - //- Dilatation - volScalarField dgdt_; - //- Blending methods blendingMethodTable blendingMethods_; @@ -366,12 +363,6 @@ public: //- Access the rate of change of the pressure inline volScalarField& dpdt(); - //- Constant access the dilatation parameter - inline const volScalarField& dgdt() const; - - //- Access the dilatation parameter - inline volScalarField& dgdt(); - //- Access a sub model between a phase pair template <class modelType> const modelType& lookupSubModel(const phasePair& key) const; diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H index 4095e3739f9..adb3813bbd3 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H @@ -54,17 +54,5 @@ inline Foam::volScalarField& Foam::phaseSystem::dpdt() return dpdt_; } -inline const Foam::volScalarField& Foam::phaseSystem::dgdt() const -{ - return dgdt_; -} - - -inline Foam::volScalarField& Foam::phaseSystem::dgdt() -{ - return dgdt_; -} - - // ************************************************************************* // diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C index 114f7416331..00ebfcffb98 100644 --- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/twoPhaseSystem/twoPhaseSystem.C @@ -170,12 +170,6 @@ void Foam::twoPhaseSystem::solve() volScalarField& alpha1 = phase1_; volScalarField& alpha2 = phase2_; - const surfaceScalarField& phi = this->phi(); - const surfaceScalarField& phi1 = phase1_.phi(); - const surfaceScalarField& phi2 = phase2_.phi(); - - const volScalarField& dgdt = this->dgdt(); - const dictionary& alphaControls = mesh.solverDict(alpha1.name()); label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles"))); @@ -184,22 +178,59 @@ void Foam::twoPhaseSystem::solve() word alphaScheme("div(phi," + alpha1.name() + ')'); word alpharScheme("div(phir," + alpha1.name() + ')'); - alpha1.correctBoundaryConditions(); + const surfaceScalarField& phi = this->phi(); + const surfaceScalarField& phi1 = phase1_.phi(); + const surfaceScalarField& phi2 = phase2_.phi(); + + // Construct the dilatation rate source term + tmp<volScalarField::DimensionedInternalField> tdgdt; + + if (phase1_.compressible() && phase2_.compressible()) + { + tdgdt = + ( + alpha1.dimensionedInternalField() + *phase2_.D()().dimensionedInternalField() + - alpha2.dimensionedInternalField() + *phase1_.D()().dimensionedInternalField() + ); + } + else if (phase1_.compressible()) + { + tdgdt = + ( + - alpha2.dimensionedInternalField() + *phase1_.D()().dimensionedInternalField() + ); + } + else if (phase2_.compressible()) + { + tdgdt = + ( + alpha1.dimensionedInternalField() + *phase2_.D()().dimensionedInternalField() + ); + } + alpha1.correctBoundaryConditions(); + surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0)))); surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phi1 - phi2); - surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0)))); + tmp<surfaceScalarField> alpha1alpha2f; if (pPrimeByA_.valid()) { + alpha1alpha2f = + fvc::interpolate(max(alpha1, scalar(0))) + *fvc::interpolate(max(alpha2, scalar(0))); + surfaceScalarField phiP ( pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf() ); - phic += alpha1f*phiP; phir += phiP; } @@ -214,7 +245,7 @@ void Foam::twoPhaseSystem::solve() mesh ), mesh, - dimensionedScalar("Sp", dgdt.dimensions(), 0.0) + dimensionedScalar("Sp", dimless/dimTime, 0.0) ); volScalarField::DimensionedInternalField Su @@ -230,16 +261,21 @@ void Foam::twoPhaseSystem::solve() fvc::div(phi)*min(alpha1, scalar(1)) ); - forAll(dgdt, celli) + if (tdgdt.valid()) { - if (dgdt[celli] > 0.0) - { - Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4); - Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4); - } - else if (dgdt[celli] < 0.0) + scalarField& dgdt = tdgdt(); + + forAll(dgdt, celli) { - Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4); + if (dgdt[celli] > 0.0) + { + Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4); + Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4); + } + else if (dgdt[celli] < 0.0) + { + Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4); + } } } @@ -336,7 +372,7 @@ void Foam::twoPhaseSystem::solve() fvScalarMatrix alpha1Eqn ( fvm::ddt(alpha1) - fvc::ddt(alpha1) - - fvm::laplacian(alpha1f*pPrimeByA_(), alpha1, "bounded") + - fvm::laplacian(alpha1alpha2f*pPrimeByA_(), alpha1, "bounded") ); alpha1Eqn.relax(); -- GitLab