zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

reactingMultiphaseEulerFoam: Added support for face-based momentum equation formulation

Browse Source 
The face-based momentum equation formulation introduced to twoPhaseEulerFoam by
commit 16f03f8a39 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.
This commit is contained in:
Henry Weller
2017-12-12 13:43:59 +00:00
parent 68afe78b9b
commit 986a879bef
10 changed files with 941 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

259
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
View File

@ -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<BasePhaseSystem>::Kd
}
template<class BasePhaseSystem>
Foam::tmp<Foam::surfaceScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Kdf
(
const Foam::phaseModel& phase
) const
{
tmp<surfaceScalarField> 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<class BasePhaseSystem>
Foam::tmp<Foam::volScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vm
@ -280,6 +356,62 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vmf
}
template<class BasePhaseSystem>
Foam::tmp<Foam::surfaceScalarField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Vmf
(
const Foam::phaseModel& phase
) const
{
tmp<surfaceScalarField> 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<class BasePhaseSystem>
Foam::tmp<Foam::volVectorField>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::F
@ -434,6 +566,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::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<BasePhaseSystem>::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<BasePhaseSystem>::setF
}
template<class BasePhaseSystem>
Foam::surfaceScalarField&
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::setFf
(
PtrList<surfaceScalarField>& 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<class BasePhaseSystem>
Foam::autoPtr<Foam::PtrList<Foam::volVectorField>>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
@ -589,6 +756,53 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField>>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Ffs() const
{
autoPtr<PtrList<surfaceScalarField>> tFfs
(
new PtrList<surfaceScalarField>(this->phases().size())
);
PtrList<surfaceScalarField>& 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<class BasePhaseSystem>
Foam::surfaceScalarField&
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::setPhiD
@ -667,6 +881,49 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
}
template<class BasePhaseSystem>
Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField>>
Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDfs
(
const PtrList<surfaceScalarField>& rAUfs
) const
{
autoPtr<PtrList<surfaceScalarField>> tphiDfs
(
new PtrList<surfaceScalarField>(this->phases().size())
);
PtrList<surfaceScalarField>& 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<class BasePhaseSystem>
bool Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::read()
{

48
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
View File

@ -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<volVectorField>& Fs, const label phasei
) const;
//- Construct element phasei of Ffs if not set and return
// Used by Ffs()
surfaceScalarField& setFf
(
PtrList<surfaceScalarField>& 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<volScalarField> Kd(const phasePairKey& key) const;
@ -177,24 +208,39 @@ public:
//- Return the drag coefficient for phase
virtual tmp<volScalarField> Kd(const phaseModel& phase) const;
//- Return the face drag coefficient for phase
virtual tmp<surfaceScalarField> Kdf(const phaseModel& phase) const;
//- Return the virtual mass coefficient
virtual tmp<volScalarField> Vm(const phasePairKey& key) const;
//- Return the face virtual mass coefficient
virtual tmp<surfaceScalarField> Vmf(const phasePairKey& key) const;
//- Return the face virtual mass force coefficient for phase
virtual tmp<surfaceScalarField> Vmf(const phaseModel& phase) const;
//- Return the combined force (lift + wall-lubrication)
virtual tmp<volVectorField> F(const phasePairKey& key) const;
//- Return the combined force (lift + wall-lubrication)
virtual autoPtr<PtrList<volVectorField>> Fs() const;
//- Return the combined force (lift + wall-lubrication)
virtual autoPtr<PtrList<surfaceScalarField>> Ffs() const;
//- Return the turbulent dispersion force on faces for phase pair
virtual autoPtr<PtrList<surfaceScalarField>> phiDs
(
const PtrList<volScalarField>& rAUs
) const;
//- Return the face based turbulent dispersion force for phase pair
virtual autoPtr<PtrList<surfaceScalarField>> phiDfs
(
const PtrList<surfaceScalarField>& rAUfs
) const;
//- Return the combined face-force (lift + wall-lubrication)
virtual tmp<surfaceScalarField> Ff(const phasePairKey& key) const;

20
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
View File

@ -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
<

1
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options
View File

@ -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 \

28
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.H
View File

@ -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<volScalarField> Kd(const phaseModel& phase) const = 0;
//- Return the face based drag coefficient for phase
virtual tmp<surfaceScalarField> Kdf(const phaseModel& phase) const = 0;
//- Return the face based virtual mass force coefficient for phase
virtual tmp<surfaceScalarField> Vmf(const phaseModel& phase) const = 0;
//- Return the combined force (lift + wall-lubrication) for phase pair
virtual autoPtr<PtrList<Foam::volVectorField>> Fs() const = 0;
//- Return the combined face force (lift + wall-lubrication)
virtual autoPtr<PtrList<Foam::surfaceScalarField>> Ffs() const = 0;
//- Return the turbulent dispersion force on faces for phase pair
virtual autoPtr<PtrList<Foam::surfaceScalarField>> phiDs
(
const PtrList<volScalarField>& rAUs
) const = 0;
//- Return the face based turbulent dispersion force on faces
// for phase pair
virtual autoPtr<PtrList<Foam::surfaceScalarField>> phiDfs
(
const PtrList<surfaceScalarField>& rAUfs
) const = 0;
//- Return true if there is mass transfer for phase
virtual bool transfersMass(const phaseModel& phase) const = 0;

7
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/DDtU.H Normal file
View File

@ -0,0 +1,7 @@
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
ddtPhis[phasei] =
(phase.phi() - phase.phi().oldTime())/runTime.deltaT();
}

87
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/UEqns.H Normal file
View File

@ -0,0 +1,87 @@
Info<< "Constructing face momentum equations" << endl;
PtrList<fvVectorMatrix> UEqns(phases.size());
{
// Update of coefficients
fluid.momentumTransfer();
PtrList<fvVectorMatrix> 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);
}
}

15
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/createDDtU.H Normal file
View File

@ -0,0 +1,15 @@
PtrList<surfaceScalarField> ddtPhis(phases.size());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
ddtPhis.set
(
phasei,
(
(phase.phi() - phase.phi().oldTime())/runTime.deltaT()
).ptr()
);
}

465
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pUf/pEqn.H Normal file
View File

@ -0,0 +1,465 @@
PtrList<surfaceScalarField> alphafs(phases.size());
PtrList<surfaceScalarField> alphaRho0fs(phases.size());
PtrList<surfaceScalarField> rAUfs(phases.size());
PtrList<surfaceScalarField> 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<surfaceScalarField> phiFs(phases.size());
{
autoPtr<PtrList<surfaceScalarField>> 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<PtrList<surfaceScalarField>> 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<surfaceScalarField> 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<surfaceScalarField> 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<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryFieldRef(),
(
phiHbyA.boundaryField() - phib
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
}
PtrList<fvScalarMatrix> 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();
}

29
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/reactingMultiphaseEulerFoam.C
View File

@ -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<Switch>("faceMomentum", false)
// );
Switch faceMomentum
(
pimple.dict().lookupOrDefault<Switch>("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"