diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
index f43a3a03544e86e992c80f99193c6de83bcb8ab9..f50bd61f9ea540d88f3b4822eb5dbc5395f98128 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C
@@ -511,17 +511,12 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
template<class BasePhaseSystem>
-Foam::autoPtr<Foam::PtrList<Foam::volVectorField> >
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
+Foam::volVectorField& Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::setF
+(
+ PtrList<volVectorField>& Fs, const label phasei
+) const
{
- autoPtr<PtrList<volVectorField> > tFs
- (
- new PtrList<volVectorField>(this->phases().size())
- );
-
- PtrList<volVectorField>& Fs = tFs();
-
- forAll(Fs, phasei)
+ if (!Fs.set(phasei))
{
Fs.set
(
@@ -543,6 +538,20 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
);
}
+ return Fs[phasei];
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::PtrList<Foam::volVectorField> >
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
+{
+ autoPtr<PtrList<volVectorField> > tFs
+ (
+ new PtrList<volVectorField>(this->phases().size())
+ );
+ PtrList<volVectorField>& Fs = tFs();
+
// Add the lift force
forAllConstIter
(
@@ -555,8 +564,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
const phasePair& pair(this->phasePairs_[liftModelIter.key()]);
- Fs[pair.phase1().index()] += F;
- Fs[pair.phase2().index()] -= F;
+ setF(Fs, pair.phase1().index()) += F;
+ setF(Fs, pair.phase2().index()) -= F;
}
// Add the wall lubrication force
@@ -572,8 +581,8 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
const phasePair&
pair(this->phasePairs_[wallLubricationModelIter.key()]);
- Fs[pair.phase1().index()] += F;
- Fs[pair.phase2().index()] -= F;
+ setF(Fs, pair.phase1().index()) += F;
+ setF(Fs, pair.phase2().index()) -= F;
}
return tFs;
@@ -581,20 +590,13 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::Fs() const
template<class BasePhaseSystem>
-Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField> >
-Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
+Foam::surfaceScalarField&
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::setPhiD
(
- const PtrList<volScalarField>& rAUs
+ PtrList<surfaceScalarField>& phiDs, const label phasei
) const
{
- autoPtr<PtrList<surfaceScalarField> > tphiDs
- (
- new PtrList<surfaceScalarField>(this->phases().size())
- );
-
- PtrList<surfaceScalarField>& phiDs = tphiDs();
-
- forAll(phiDs, phasei)
+ if (!phiDs.set(phasei))
{
phiDs.set
(
@@ -603,7 +605,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
(
IOobject
(
- liftModel::typeName + ":F",
+ turbulentDispersionModel::typeName + ":phiD",
this->mesh_.time().timeName(),
this->mesh_,
IOobject::NO_READ,
@@ -621,6 +623,23 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
);
}
+ return phiDs[phasei];
+}
+
+
+template<class BasePhaseSystem>
+Foam::autoPtr<Foam::PtrList<Foam::surfaceScalarField> >
+Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
+(
+ const PtrList<volScalarField>& rAUs
+) const
+{
+ autoPtr<PtrList<surfaceScalarField> > tphiDs
+ (
+ new PtrList<surfaceScalarField>(this->phases().size())
+ );
+ PtrList<surfaceScalarField>& phiDs = tphiDs();
+
// Add the turbulent dispersion force
forAllConstIter
(
@@ -638,9 +657,9 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiDs
fvc::snGrad(pair.phase1())*this->mesh_.magSf()
);
- phiDs[pair.phase1().index()] +=
+ setPhiD(phiDs, pair.phase1().index()) +=
fvc::interpolate(rAUs[pair.phase1().index()]*D)*snGradAlpha1;
- phiDs[pair.phase2().index()] -=
+ setPhiD(phiDs, pair.phase2().index()) -=
fvc::interpolate(rAUs[pair.phase2().index()]*D)*snGradAlpha1;
}
diff --git a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
index 980f854a9e406f452dcfebbee7b2661228a331dc..2f787ed004a96329fc1cd8b11fd81d517286ac85 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H
@@ -133,6 +133,20 @@ private:
//- Turbulent dispersion models
turbulentDispersionModelTable turbulentDispersionModels_;
+ //- Construct element phasei of Fs if not set and return
+ // Used by Fs()
+ volVectorField& setF
+ (
+ PtrList<volVectorField>& Fs, const label phasei
+ ) const;
+
+ //- Construct element phasei of phiDs if not set and return
+ // Used by phiDs()
+ surfaceScalarField& setPhiD
+ (
+ PtrList<surfaceScalarField>& phiDs, const label phasei
+ ) const;
+
public:
diff --git a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
index 136d595418b27936ec0feefa9a3ba50b2840ae91..e23746a3adfa9040909055a546912d4b3b59a948 100644
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
@@ -38,22 +38,51 @@ forAll(phases, phasei)
PtrList<surfaceScalarField> phiFs(phases.size());
{
autoPtr<PtrList<volVectorField> > Fs = fluid.Fs();
- autoPtr<PtrList<surfaceScalarField> > phiDs = fluid.phiDs(rAUs);
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
- phiFs.set
- (
- phasei,
- new surfaceScalarField
+ if (Fs().set(phasei))
+ {
+ phiFs.set
(
- IOobject::groupName("phiF", phase.name()),
- (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
- + phiDs()[phasei]
- )
- );
+ phasei,
+ new surfaceScalarField
+ (
+ IOobject::groupName("phiF", phase.name()),
+ (fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
+ )
+ );
+ }
+ }
+}
+{
+ autoPtr<PtrList<surfaceScalarField> > phiDs = fluid.phiDs(rAUs);
+
+ 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]
+ )
+ );
+ }
+ }
}
}
@@ -101,12 +130,12 @@ while (pimple.correct())
ghf*fvc::snGrad(rho)*mesh.magSf()
);
- PtrList<surfaceScalarField> phigs(phases.size());
+ PtrList<surfaceScalarField> phigFs(phases.size());
forAll(phases, phasei)
{
phaseModel& phase = phases[phasei];
- phigs.set
+ phigFs.set
(
phasei,
(
@@ -118,6 +147,11 @@ while (pimple.correct())
)
).ptr()
);
+
+ if (phiFs.set(phasei))
+ {
+ phigFs[phasei] += phiFs[phasei];
+ }
}
PtrList<surfaceScalarField> phiHbyAs(phases.size());
@@ -177,8 +211,7 @@ while (pimple.correct())
MRF.absolute(phase.phi().oldTime())
- (fvc::interpolate(phase.U().oldTime()) & mesh.Sf())
)/runTime.deltaT()
- - phiFs[phasei]
- - phigs[phasei]
+ - phigFs[phasei]
)
);
@@ -389,8 +422,7 @@ while (pimple.correct())
+ fvc::reconstruct
(
alpharAUfs[phasei]*mSfGradp
- - phiFs[phasei]
- - phigs[phasei]
+ - phigFs[phasei]
);
phase.U().correctBoundaryConditions();
fvOptions.correct(phase.U());
diff --git a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
index f7590f98f3f7afd7694aa2552d49342c5f26fc9a..98192e477f8097c3691a1d9a90b6a882a1212002 100644
--- a/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/reactingMultiphaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
@@ -151,16 +151,13 @@ heatTransfer
);
lift
-(
-);
+();
wallLubrication
-(
-);
+();
turbulentDispersion
-(
-);
+();
// Minimum allowable pressure
pMin 10000;