From 986a879bef3eb2cef91a428c72b46ce5298eea22 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 12 Dec 2017 13:43:59 +0000 Subject: [PATCH] reactingMultiphaseEulerFoam: Added support for face-based momentum equation formulation The face-based momentum equation formulation introduced to twoPhaseEulerFoam by commit 16f03f8a393bd4b6e70b548812777bb10a3e1dda has proven particularly valuable for bubbly flow simulations. The formulation is also available for reactingTwoPhaseEulerFoam and this patch adds the the same capability to reactingMultiphaseEulerFoam. It be switched on by setting the optional faceMomentum entry in the PIMPLE sub-dictionary in fvSolution: PIMPLE { nOuterCorrectors 3; nCorrectors 1; nNonOrthogonalCorrectors 0; faceMomentum yes; } Patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden - Rossendorf (HZDR) and VTT Technical Research Centre of Finland Ltd. --- .../MomentumTransferPhaseSystem.C | 259 +++++++++- .../MomentumTransferPhaseSystem.H | 48 +- .../phaseSystems/phaseSystem/phaseSystem.H | 20 +- .../reactingMultiphaseEulerFoam/Make/options | 1 + .../multiphaseSystem/multiphaseSystem.H | 28 +- .../reactingMultiphaseEulerFoam/pUf/DDtU.H | 7 + .../reactingMultiphaseEulerFoam/pUf/UEqns.H | 87 ++++ .../pUf/createDDtU.H | 15 + .../reactingMultiphaseEulerFoam/pUf/pEqn.H | 465 ++++++++++++++++++ .../reactingMultiphaseEulerFoam.C | 29 +- 10 files changed, 941 insertions(+), 18 deletions(-) create mode 100644 applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/DDtU.H create mode 100644 applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H create mode 100644 applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/createDDtU.H create mode 100644 applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index b9859eb41d..da99791b64 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -100,6 +100,16 @@ MomentumTransferPhaseSystem dragModelIter()->K() ) ); + + Kdfs_.insert + ( + pair, + new surfaceScalarField + ( + IOobject::groupName("Kdf", pair.name()), + dragModelIter()->Kf() + ) + ); } forAllConstIter @@ -120,6 +130,16 @@ MomentumTransferPhaseSystem virtualMassModelIter()->K() ) ); + + Vmfs_.insert + ( + pair, + new surfaceScalarField + ( + IOobject::groupName("Vmf", pair.name()), + virtualMassModelIter()->Kf() + ) + ); } } @@ -212,6 +232,62 @@ Foam::MomentumTransferPhaseSystem::Kd } +template +Foam::tmp +Foam::MomentumTransferPhaseSystem::Kdf +( + const Foam::phaseModel& phase +) const +{ + tmp tKdf + ( + new surfaceScalarField + ( + IOobject + ( + IOobject::groupName("Kdf", phase.name()), + this->mesh_.time().timeName(), + this->mesh_ + ), + this->mesh_, + dimensionedScalar + ( + IOobject::groupName("Kdf", phase.name()), + dimensionSet(1, -3, -1, 0, 0), + 0 + ) + ) + ); + + forAllConstIter + ( + phaseSystem::KdfTable, + Kdfs_, + KdfIter + ) + { + const surfaceScalarField& Kf(*KdfIter()); + + const phasePair& pair(this->phasePairs_[KdfIter.key()]); + + const phaseModel* phase1 = &pair.phase1(); + const phaseModel* phase2 = &pair.phase2(); + + forAllConstIter(phasePair, pair, iter) + { + if (phase1 == &phase) + { + tKdf.ref() += Kf; + } + + Swap(phase1, phase2); + } + } + + return tKdf; +} + + template Foam::tmp Foam::MomentumTransferPhaseSystem::Vm @@ -280,6 +356,62 @@ Foam::MomentumTransferPhaseSystem::Vmf } +template +Foam::tmp +Foam::MomentumTransferPhaseSystem::Vmf +( + const Foam::phaseModel& phase +) const +{ + tmp tVmf + ( + new surfaceScalarField + ( + IOobject + ( + IOobject::groupName("Vmf", phase.name()), + this->mesh_.time().timeName(), + this->mesh_ + ), + this->mesh_, + dimensionedScalar + ( + IOobject::groupName("Vmf", phase.name()), + virtualMassModel::dimK, + 0 + ) + ) + ); + + forAllConstIter + ( + phaseSystem::VmfTable, + Vmfs_, + VmfIter + ) + { + const surfaceScalarField& Vmf(*VmfIter()); + + const phasePair& pair(this->phasePairs_[VmfIter.key()]); + + const phaseModel* phase1 = &pair.phase1(); + const phaseModel* phase2 = &pair.phase2(); + + forAllConstIter(phasePair, pair, iter) + { + if (phase1 == &phase) + { + tVmf.ref() += Vmf; + } + + Swap(phase1, phase2); + } + } + + return tVmf; +} + + template Foam::tmp Foam::MomentumTransferPhaseSystem::F @@ -434,6 +566,7 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const ) { *Kds_[dragModelIter.key()] = dragModelIter()->K(); + *Kdfs_[dragModelIter.key()] = dragModelIter()->Kf(); } // Add the implicit part of the drag force @@ -470,6 +603,7 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const ) { *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K(); + *Vmfs_[virtualMassModelIter.key()] = virtualMassModelIter()->Kf(); } // Add the virtual mass force @@ -542,6 +676,39 @@ Foam::volVectorField& Foam::MomentumTransferPhaseSystem::setF } +template +Foam::surfaceScalarField& +Foam::MomentumTransferPhaseSystem::setFf +( + PtrList& Ffs, const label phasei +) const +{ + if (!Ffs.set(phasei)) + { + Ffs.set + ( + phasei, + new surfaceScalarField + ( + IOobject + ( + liftModel::typeName + ":Ff", + this->mesh_.time().timeName(), + this->mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + this->mesh_, + dimensionedScalar("zero", dimArea*liftModel::dimF, Zero) + ) + ); + } + + return Ffs[phasei]; +} + + template Foam::autoPtr> Foam::MomentumTransferPhaseSystem::Fs() const @@ -589,6 +756,53 @@ Foam::MomentumTransferPhaseSystem::Fs() const } +template +Foam::autoPtr> +Foam::MomentumTransferPhaseSystem::Ffs() const +{ + autoPtr> tFfs + ( + new PtrList(this->phases().size()) + ); + PtrList& Ffs = tFfs(); + + // Add the lift force + forAllConstIter + ( + liftModelTable, + liftModels_, + liftModelIter + ) + { + const surfaceScalarField Ff(liftModelIter()->Ff()); + + const phasePair& pair(this->phasePairs_[liftModelIter.key()]); + + setFf(Ffs, pair.phase1().index()) += Ff; + setFf(Ffs, pair.phase2().index()) -= Ff; + } + + // Add the wall lubrication force + forAllConstIter + ( + wallLubricationModelTable, + wallLubricationModels_, + wallLubricationModelIter + ) + { + const surfaceScalarField Ff(wallLubricationModelIter()->Ff()); + + const phasePair& + pair(this->phasePairs_[wallLubricationModelIter.key()]); + + setFf(Ffs, pair.phase1().index()) += Ff; + setFf(Ffs, pair.phase2().index()) -= Ff; + } + + return tFfs; +} + + template Foam::surfaceScalarField& Foam::MomentumTransferPhaseSystem::setPhiD @@ -667,6 +881,49 @@ Foam::MomentumTransferPhaseSystem::phiDs } +template +Foam::autoPtr> +Foam::MomentumTransferPhaseSystem::phiDfs +( + const PtrList& rAUfs +) const +{ + autoPtr> tphiDfs + ( + new PtrList(this->phases().size()) + ); + PtrList& phiDfs = tphiDfs(); + + // Add the face based turbulent dispersion force + forAllConstIter + ( + turbulentDispersionModelTable, + turbulentDispersionModels_, + turbulentDispersionModelIter + ) + { + const phasePair& + pair(this->phasePairs_[turbulentDispersionModelIter.key()]); + + const surfaceScalarField Df + ( + fvc::interpolate(turbulentDispersionModelIter()->D()) + ); + const surfaceScalarField snGradAlpha1 + ( + fvc::snGrad(pair.phase1())*this->mesh_.magSf() + ); + + setPhiD(phiDfs, pair.phase1().index()) += + rAUfs[pair.phase1().index()]*Df*snGradAlpha1; + setPhiD(phiDfs, pair.phase2().index()) -= + rAUfs[pair.phase2().index()]*Df*snGradAlpha1; + } + + return tphiDfs; +} + + template bool Foam::MomentumTransferPhaseSystem::read() { diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H index e6a3209f1c..e6e0d02ecf 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -113,9 +113,15 @@ private: //- Drag coefficients phaseSystem::KdTable Kds_; + //- Face drag coefficients + phaseSystem::KdfTable Kdfs_; + //- Virtual mass coefficients phaseSystem::VmTable Vms_; + //- Face virtual mass coefficients + phaseSystem::VmfTable Vmfs_; + // Sub Models //- Drag models @@ -140,6 +146,13 @@ private: PtrList& Fs, const label phasei ) const; + //- Construct element phasei of Ffs if not set and return + // Used by Ffs() + surfaceScalarField& setFf + ( + PtrList& Ffs, const label phasei + ) const; + //- Construct element phasei of phiDs if not set and return // Used by phiDs() surfaceScalarField& setPhiD @@ -168,6 +181,24 @@ public: return Kds_; } + //- Constant access to face drag coefficients + virtual const phaseSystem::KdfTable& Kdfs() const + { + return Kdfs_; + } + + //- Constant access to virtual mass force coefficients + virtual const phaseSystem::VmTable& Vms() const + { + return Vms_; + } + + //- Constant access to virtual mass force coefficients + virtual const phaseSystem::VmfTable& Vmfs() const + { + return Vmfs_; + } + //- Return the drag coefficient virtual tmp Kd(const phasePairKey& key) const; @@ -177,24 +208,39 @@ public: //- Return the drag coefficient for phase virtual tmp Kd(const phaseModel& phase) const; + //- Return the face drag coefficient for phase + virtual tmp Kdf(const phaseModel& phase) const; + //- Return the virtual mass coefficient virtual tmp Vm(const phasePairKey& key) const; //- Return the face virtual mass coefficient virtual tmp Vmf(const phasePairKey& key) const; + //- Return the face virtual mass force coefficient for phase + virtual tmp Vmf(const phaseModel& phase) const; + //- Return the combined force (lift + wall-lubrication) virtual tmp F(const phasePairKey& key) const; //- Return the combined force (lift + wall-lubrication) virtual autoPtr> Fs() const; + //- Return the combined force (lift + wall-lubrication) + virtual autoPtr> Ffs() const; + //- Return the turbulent dispersion force on faces for phase pair virtual autoPtr> phiDs ( const PtrList& rAUs ) const; + //- Return the face based turbulent dispersion force for phase pair + virtual autoPtr> phiDfs + ( + const PtrList& rAUfs + ) const; + //- Return the combined face-force (lift + wall-lubrication) virtual tmp Ff(const phasePairKey& key) const; diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H index 9bfaebd9f7..9c637a428a 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2015-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -82,6 +82,15 @@ public: > KdTable; + typedef + HashPtrTable + < + surfaceScalarField, + phasePairKey, + phasePairKey::hash + > + KdfTable; + typedef HashPtrTable < @@ -91,6 +100,15 @@ public: > VmTable; + typedef + HashPtrTable + < + surfaceScalarField, + phasePairKey, + phasePairKey::hash + > + VmfTable; + typedef HashPtrTable < diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options index 656a9941a6..06caef01de 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options @@ -3,6 +3,7 @@ EXE_INC = \ -I../phaseSystems/lnInclude \ -I../interfacialModels/lnInclude \ -I../interfacialCompositionModels/lnInclude \ + -ImultiphaseCompressibleTurbulenceModels/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/transportModels/compressible/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H index dca55bf873..245450dabd 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -163,18 +163,44 @@ public: //- Return the drag coefficient for all phase-pairs virtual const phaseSystem::KdTable& Kds() const = 0; + //- Return the face drag coefficient for all phase-pairs + virtual const phaseSystem::KdfTable& Kdfs() const = 0; + + //- Return the virtual mass force coefficient for all phase-pairs + virtual const phaseSystem::VmTable& Vms() const = 0; + + //- Return the face based virtual mass force coefficient + // for all phase-pairs + virtual const phaseSystem::VmfTable& Vmfs() const = 0; + //- Return the drag coefficient for phase virtual tmp Kd(const phaseModel& phase) const = 0; + //- Return the face based drag coefficient for phase + virtual tmp Kdf(const phaseModel& phase) const = 0; + + //- Return the face based virtual mass force coefficient for phase + virtual tmp Vmf(const phaseModel& phase) const = 0; + //- Return the combined force (lift + wall-lubrication) for phase pair virtual autoPtr> Fs() const = 0; + //- Return the combined face force (lift + wall-lubrication) + virtual autoPtr> Ffs() const = 0; + //- Return the turbulent dispersion force on faces for phase pair virtual autoPtr> phiDs ( const PtrList& rAUs ) const = 0; + //- Return the face based turbulent dispersion force on faces + // for phase pair + virtual autoPtr> phiDfs + ( + const PtrList& rAUfs + ) const = 0; + //- Return true if there is mass transfer for phase virtual bool transfersMass(const phaseModel& phase) const = 0; diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/DDtU.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/DDtU.H new file mode 100644 index 0000000000..7dfa823102 --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/DDtU.H @@ -0,0 +1,7 @@ + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + ddtPhis[phasei] = + (phase.phi() - phase.phi().oldTime())/runTime.deltaT(); + } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H new file mode 100644 index 0000000000..c513c24266 --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H @@ -0,0 +1,87 @@ +Info<< "Constructing face momentum equations" << endl; + +PtrList UEqns(phases.size()); + +{ + // Update of coefficients + fluid.momentumTransfer(); + + PtrList UgradUs(phases.size()); + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + UgradUs.set + ( + phasei, + new fvVectorMatrix + ( + fvm::div(phase.phi(), phase.U()) + - fvm::Sp(fvc::div(phase.phi()), phase.U()) + + MRF.DDt(phase.U()) + ) + ); + } + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + const volScalarField& alpha = phase; + volScalarField& rho = phase.thermo().rho(); + const surfaceScalarField& alphaRhoPhi = phase.alphaRhoPhi(); + volVectorField& U = phase.U(); + + + UEqns.set + ( + phasei, + new fvVectorMatrix + ( + fvm::div(alphaRhoPhi, U) - fvm::Sp(fvc::div(alphaRhoPhi), U) + + MRF.DDt(alpha*rho, U) + + phase.turbulence().divDevRhoReff(U) + + fvm::SuSp(fvOptions(alpha, rho)&rho,U) + == + fvOptions(alpha, rho, U) + ) + ); + + // Add gradient part of virtual mass force + forAllConstIter + ( + phaseSystem::VmTable, + fluid.Vms(), + VmIter + ) + { + const volScalarField& Vm(*VmIter()); + + const phasePair& pair(fluid.phasePairs()[VmIter.key()]); + + const phaseModel* phase1 = &pair.phase1(); + const phaseModel* phase2 = &pair.phase2(); + + forAllConstIter(phasePair, pair, iter) + { + if (phase1 == &phase) + { + UEqns[phasei] += + Vm + *( + UgradUs[phasei] + - (UgradUs[phase2->index()] & phase2->U()) + ); + } + + Swap(phase1, phase2); + } + } + + UEqns[phasei].relax(); + fvOptions.constrain(UEqns[phasei]); + U.correctBoundaryConditions(); + fvOptions.correct(U); + } +} diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/createDDtU.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/createDDtU.H new file mode 100644 index 0000000000..808bd82d0a --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/createDDtU.H @@ -0,0 +1,15 @@ + PtrList ddtPhis(phases.size()); + + forAll(phases, phasei) + { + + phaseModel& phase = phases[phasei]; + + ddtPhis.set + ( + phasei, + ( + (phase.phi() - phase.phi().oldTime())/runTime.deltaT() + ).ptr() + ); + } diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H new file mode 100644 index 0000000000..f2aa61cdcd --- /dev/null +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H @@ -0,0 +1,465 @@ +PtrList alphafs(phases.size()); +PtrList alphaRho0fs(phases.size()); +PtrList rAUfs(phases.size()); +PtrList alpharAUfs(phases.size()); + +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + + 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() + ); + + // add implicit part of drag and virtual mass force + rAUfs.set + ( + phasei, + new surfaceScalarField + ( + IOobject::groupName("rAUf", phase.name()), + 1.0 + /( + fvc::interpolate(UEqns[phasei].A()) + + (alphaRho0fs[phasei] + fluid.Vmf(phase))/runTime.deltaT() + + fluid.Kdf(phase) + ) + ) + ); + + alpharAUfs.set + ( + phasei, + ( + max(alphafs[phasei], phase.residualAlpha())*rAUfs[phasei] + ).ptr() + ); +} + +// Lift, wall-lubrication and turbulent diffusion fluxes +PtrList phiFs(phases.size()); +{ + autoPtr> Fs = fluid.Ffs(); + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + if (Fs().set(phasei)) + { + phiFs.set + ( + phasei, + new surfaceScalarField + ( + IOobject::groupName("phiF", phase.name()), + rAUfs[phasei]*Fs()[phasei] + ) + ); + } + } +} +{ + autoPtr> phiDs = fluid.phiDfs(rAUfs); + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + if (phiDs().set(phasei)) + { + if (phiFs.set(phasei)) + { + phiFs[phasei] += phiDs()[phasei]; + } + else + { + phiFs.set + ( + phasei, + new surfaceScalarField + ( + IOobject::groupName("phiF", phase.name()), + phiDs()[phasei] + ) + ); + } + } + } +} + +// --- 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*gh; + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + // Correct fixed-flux BCs to be consistent with the velocity BCs + MRF.correctBoundaryFlux(phase.U(), phase.phi()); + } + + surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + PtrList phigFs(phases.size()); + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + phigFs.set + ( + phasei, + ( + alpharAUfs[phasei] + *( + ghSnGradRho + - (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf()) + - fluid.surfaceTension(phase)*mesh.magSf() + ) + ).ptr() + ); + + if (phiFs.set(phasei)) + { + phigFs[phasei] += phiFs[phasei]; + } + } + + PtrList phiHbyAs(phases.size()); + + surfaceScalarField phiHbyA + ( + IOobject + ( + "phiHbyA", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0) + ); + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + phiHbyAs.set + ( + phasei, + new surfaceScalarField + ( + IOobject::groupName("phiHbyA", phase.name()), + rAUfs[phasei] + *( + fvc::flux(UEqns[phasei].H()) + + alphaRho0fs[phasei] + *MRF.absolute(phase.phi().oldTime())/runTime.deltaT() + ) + - phigFs[phasei] + ) + ); + + // Add explicit part of the drag + forAllConstIter + ( + phaseSystem::KdfTable, + fluid.Kdfs(), + KdfIter + ) + { + const surfaceScalarField& Kf(*KdfIter()); + + const phasePair& pair(fluid.phasePairs()[KdfIter.key()]); + + const phaseModel* phase1 = &pair.phase1(); + const phaseModel* phase2 = &pair.phase2(); + + forAllConstIter(phasePair, pair, iter) + { + if (phase1 == &phase) + { + phiHbyAs[phasei] += + rAUfs[phasei]*Kf*MRF.absolute(phase2->phi()); + } + + Swap(phase1, phase2); + } + } + + // Add explicit part of the virtual mass force + forAllConstIter + ( + phaseSystem::VmfTable, + fluid.Vmfs(), + VmfIter + ) + { + const surfaceScalarField& Vmf(*VmfIter()); + + const phasePair& pair(fluid.phasePairs()[VmfIter.key()]); + + const phaseModel* phase1 = &pair.phase1(); + const phaseModel* phase2 = &pair.phase2(); + + forAllConstIter(phasePair, pair, iter) + { + if (phase1 == &phase) + { + phiHbyAs[phasei] += + rAUfs[phasei] + *( + Vmf*MRF.absolute(phase1->phi()().oldTime()) + /runTime.deltaT() + + Vmf*ddtPhis[phase2->index()] + ); + } + + Swap(phase1, phase2); + } + } + + phiHbyA += alphafs[phasei]*phiHbyAs[phasei]; + } + + MRF.makeRelative(phiHbyA); + + // Construct pressure "diffusivity" + surfaceScalarField rAUf + ( + IOobject + ( + "rAUf", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0) + ); + + forAll(phases, phasei) + { + rAUf += alphafs[phasei]*alpharAUfs[phasei]; + } + rAUf = mag(rAUf); + + + // Update the fixedFluxPressure BCs to ensure flux consistency + { + surfaceScalarField::Boundary phib(phi.boundaryField()); + phib = 0; + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField(); + } + + setSnGrad + ( + p_rgh.boundaryFieldRef(), + ( + phiHbyA.boundaryField() - phib + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); + } + + PtrList pEqnComps(phases.size()); + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + volScalarField& rho = phase.thermo().rho(); + + if (phase.compressible()) + { + if (pimple.transonic()) + { + surfaceScalarField phid + ( + IOobject::groupName("phid", phase.name()), + fvc::interpolate(phase.thermo().psi())*phase.phi() + ); + + pEqnComps.set + ( + phasei, + ( + ( + phase.continuityError() + - fvc::Sp + ( + fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), + rho + ) + )/rho + + correction + ( + (alpha/rho)* + ( + phase.thermo().psi()*fvm::ddt(p_rgh) + + fvm::div(phid, p_rgh) + - fvm::Sp(fvc::div(phid), p_rgh) + ) + ) + ).ptr() + ); + + deleteDemandDrivenData + ( + pEqnComps[phasei].faceFluxCorrectionPtr() + ); + pEqnComps[phasei].relax(); + } + else + { + pEqnComps.set + ( + phasei, + ( + ( + phase.continuityError() + - fvc::Sp + ( + (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())), + rho + ) + )/rho + + (alpha*phase.thermo().psi()/rho) + *correction(fvm::ddt(p_rgh)) + ).ptr() + ); + } + } + else + { + pEqnComps.set + ( + phasei, + fvm::Su(-(fvOptions(alpha, rho)&rho)/rho, p_rgh).ptr() + ); + } + + if (fluid.transfersMass(phase)) + { + if (pEqnComps.set(phasei)) + { + pEqnComps[phasei] -= fluid.dmdt(phase)/rho; + } + else + { + pEqnComps.set + ( + phasei, + fvm::Su(-fluid.dmdt(phase)/rho, p_rgh) + ); + } + } + } + + // Cache p prior to solve for density update + volScalarField p_rgh_0(p_rgh); + + // Iterate over the pressure equation to correct for non-orthogonality + while (pimple.correctNonOrthogonal()) + { + // Construct the transport part of the pressure equation + fvScalarMatrix pEqnIncomp + ( + fvc::div(phiHbyA) + - fvm::laplacian(rAUf, p_rgh) + ); + + { + fvScalarMatrix pEqn(pEqnIncomp); + + forAll(phases, phasei) + { + if (pEqnComps.set(phasei)) + { + pEqn += pEqnComps[phasei]; + } + } + + solve + ( + pEqn, + mesh.solver(p_rgh.select(pimple.finalInnerIter())) + ); + } + + // Correct fluxes and velocities on last non-orthogonal iteration + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqnIncomp.flux(); + + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp; + + // Set the phase dilatation rates + if (pEqnComps.set(phasei)) + { + phase.divU(-pEqnComps[phasei] & p_rgh); + } + } + + // Optionally relax pressure for velocity correction + p_rgh.relax(); + + mSfGradp = pEqnIncomp.flux()/rAUf; + + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + + phase.U() = fvc::reconstruct(MRF.absolute(phase.phi())); + phase.U().correctBoundaryConditions(); + fvOptions.correct(phase.U()); + } + } + } + + // Update and limit the static pressure + p = max(p_rgh + rho*gh, pMin); + + // Limit p_rgh + p_rgh = p - rho*gh; + + // Update densities from change in p_rgh + forAll(phases, phasei) + { + phaseModel& phase = phases[phasei]; + phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); + } + + // Correct p_rgh for consistency with p and the updated densities + rho = fluid.rho(); + p_rgh = p - rho*gh; + p_rgh.correctBoundaryConditions(); +} diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C index 870645c5c1..7e07686403 100644 --- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C +++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -35,6 +35,7 @@ Description #include "fvCFD.H" #include "multiphaseSystem.H" +#include "phaseCompressibleTurbulenceModel.H" #include "pimpleControl.H" #include "localEulerDdtScheme.H" #include "fvcSmooth.H" @@ -59,10 +60,10 @@ int main(int argc, char *argv[]) #include "setInitialDeltaT.H" } - // Switch faceMomentum - // ( - // pimple.dict().lookupOrDefault("faceMomentum", false) - // ); + Switch faceMomentum + ( + pimple.dict().lookupOrDefault("faceMomentum", false) + ); // Switch implicitPhasePressure // ( @@ -72,7 +73,7 @@ int main(int argc, char *argv[]) // ) // ); - //#include "pUf/createDDtU.H" + #include "pUf/createDDtU.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -108,14 +109,14 @@ int main(int argc, char *argv[]) #include "YEqns.H" - // if (faceMomentum) - // { - // #include "pUf/UEqns.H" - // #include "EEqns.H" - // #include "pUf/pEqn.H" - // #include "pUf/DDtU.H" - // } - // else + if (faceMomentum) + { + #include "pUf/UEqns.H" + #include "EEqns.H" + #include "pUf/pEqn.H" + #include "pUf/DDtU.H" + } + else { #include "pU/UEqns.H" #include "EEqns.H"