songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

reactingEulerFoam: Use PtrListDictionary for list/table of phases

Browse Source 
This makes looping over the phases much simpler which maintaining
support for phase-name lookup.
This commit is contained in:
Henry Weller 2015-09-16 21:29:09 +01:00
parent 7b5d6114ad
commit c31789c34c
18 changed files with 93 additions and 420 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatAndMassTransferPhaseSystem/HeatAndMassTransferPhaseSystem.C
View File

@ -270,7 +270,7 @@ Foam::HeatAndMassTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
forAllConstIter
(
phaseSystem::phaseModelTable,
phaseSystem::phaseModelList,
this->phaseModels_,
phaseModelIter
)

2
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/HeatTransferPhaseSystem/HeatTransferPhaseSystem.C
View File

@ -122,7 +122,7 @@ Foam::HeatTransferPhaseSystem<BasePhaseSystem>::heatTransfer() const
forAllConstIter
(
phaseSystem::phaseModelTable,
phaseSystem::phaseModelList,
this->phaseModels_,
phaseModelIter
)

2
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/InterfaceCompositionPhaseChangePhaseSystem/InterfaceCompositionPhaseChangePhaseSystem.C
View File

@ -71,7 +71,7 @@ massTransfer() const
forAllConstIter
(
phaseSystem::phaseModelTable,
phaseSystem::phaseModelList,
this->phaseModels_,
phaseModelIter
)

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

@ -415,7 +415,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::momentumTransfer() const
forAllConstIter
(
phaseSystem::phaseModelTable,
phaseSystem::phaseModelList,
this->phaseModels_,
phaseModelIter
)

2
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/PhaseSystems/ThermalPhaseChangePhaseSystem/ThermalPhaseChangePhaseSystem.C
View File

@ -66,7 +66,7 @@ Foam::ThermalPhaseChangePhaseSystem<BasePhaseSystem>::massTransfer() const
forAllConstIter
(
phaseSystem::phaseModelTable,
phaseSystem::phaseModelList,
this->phaseModels_,
phaseModelIter
)

22
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystem.C
View File

@ -43,10 +43,10 @@ const Foam::word Foam::phaseSystem::propertiesName("phaseProperties");
Foam::tmp<Foam::surfaceScalarField> Foam::phaseSystem::calcPhi
(
const phaseModelTable& phaseModels
const phaseModelList& phaseModels
) const
{
phaseModelTable::const_iterator phaseModelIter = phaseModels.begin();
phaseModelList::const_iterator phaseModelIter = phaseModels.begin();
tmp<surfaceScalarField> tmpPhi
(
@ -196,7 +196,7 @@ Foam::phaseSystem::~phaseSystem()
Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
{
phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
phaseModelList::const_iterator phaseModelIter = phaseModels_.begin();
tmp<volScalarField> tmpRho
(
@ -214,7 +214,7 @@ Foam::tmp<Foam::volScalarField> Foam::phaseSystem::rho() const
Foam::tmp<Foam::volVectorField> Foam::phaseSystem::U() const
{
phaseModelTable::const_iterator phaseModelIter = phaseModels_.begin();
phaseModelList::const_iterator phaseModelIter = phaseModels_.begin();
tmp<volVectorField> tmpU
(
@ -266,7 +266,7 @@ void Foam::phaseSystem::solve()
void Foam::phaseSystem::correct()
{
forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllIter(phaseModelList, phaseModels_, phaseModelIter)
{
phaseModelIter().correct();
}
@ -277,7 +277,7 @@ void Foam::phaseSystem::correctKinematics()
{
bool updateDpdt = false;
forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllIter(phaseModelList, phaseModels_, phaseModelIter)
{
phaseModelIter().correctKinematics();
@ -287,14 +287,14 @@ void Foam::phaseSystem::correctKinematics()
// Update the pressure time-derivative if required
if (updateDpdt)
{
dpdt_ = fvc::ddt(phaseModels_.begin()()->thermo().p());
dpdt_ = fvc::ddt(phaseModels_.begin()().thermo().p());
}
}
void Foam::phaseSystem::correctThermo()
{
forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllIter(phaseModelList, phaseModels_, phaseModelIter)
{
phaseModelIter().correctThermo();
}
@ -303,7 +303,7 @@ void Foam::phaseSystem::correctThermo()
void Foam::phaseSystem::correctTurbulence()
{
forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllIter(phaseModelList, phaseModels_, phaseModelIter)
{
phaseModelIter().correctTurbulence();
}
@ -312,7 +312,7 @@ void Foam::phaseSystem::correctTurbulence()
void Foam::phaseSystem::correctEnergyTransport()
{
forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllIter(phaseModelList, phaseModels_, phaseModelIter)
{
phaseModelIter().correctEnergyTransport();
}
@ -325,7 +325,7 @@ bool Foam::phaseSystem::read()
{
bool readOK = true;
forAllIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllIter(phaseModelList, phaseModels_, phaseModelIter)
{
readOK &= phaseModelIter().read();
}

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

@ -42,7 +42,7 @@ SourceFiles
#include "phasePair.H"
#include "orderedPhasePair.H"
#include "HashPtrTable.H"
#include "PtrDictionary.H"
#include "PtrListDictionary.H"
#include "IOMRFZoneList.H"
#include "fvIOoptionList.H"
@ -118,7 +118,7 @@ public:
>
massTransferTable;
typedef PtrDictionary<phaseModel> phaseModelTable;
typedef PtrListDictionary<phaseModel> phaseModelList;
protected:
@ -162,7 +162,7 @@ protected:
const fvMesh& mesh_;
//- Phase models
phaseModelTable phaseModels_;
phaseModelList phaseModels_;
//- Phase pairs
phasePairTable phasePairs_;
@ -197,7 +197,7 @@ protected:
//- Calculate and return the mixture flux
tmp<surfaceScalarField> calcPhi
(
const phaseModelTable& phaseModels
const phaseModelList& phaseModels
) const;
//- Generate pairs
@ -284,10 +284,10 @@ public:
inline const fvMesh& mesh() const;
//- Constant access the phase models
inline const phaseModelTable& phases() const;
inline const phaseModelList& phases() const;
//- Access the phase models
inline phaseModelTable& phases();
inline phaseModelList& phases();
//- Constant access the phase pairs
inline const phasePairTable& phasePairs() const;

4
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemI.H
View File

@ -31,14 +31,14 @@ inline const Foam::fvMesh& Foam::phaseSystem::mesh() const
}
inline const Foam::phaseSystem::phaseModelTable&
inline const Foam::phaseSystem::phaseModelList&
Foam::phaseSystem::phases() const
{
return phaseModels_;
}
inline Foam::phaseSystem::phaseModelTable&
inline Foam::phaseSystem::phaseModelList&
Foam::phaseSystem::phases()
{
return phaseModels_;

2
applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.H
View File

@ -151,7 +151,7 @@ void Foam::phaseSystem::generatePairsAndSubModels
HashTable<autoPtr<modelType>, phasePairKey, phasePairKey::hash>
modelTypeTable;
forAllConstIter(phaseModelTable, phaseModels_, phaseModelIter)
forAllConstIter(phaseModelList, phaseModels_, phaseModelIter)
{
modelTypeTable tempModels;
generatePairsAndSubModels

4
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/CourantNo.H
View File

@ -39,12 +39,12 @@ if (mesh.nInternalFaces())
fvc::surfaceSum(mag(phi))().internalField()
);
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
sumPhi = max
(
sumPhi,
fvc::surfaceSum(mag(iter().phi()))().internalField()
fvc::surfaceSum(mag(fluid.phases()[phasei].phi()))().internalField()
);
}

9
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/EEqns.H
View File

@ -6,9 +6,10 @@
phaseSystem::heatTransferTable& heatTransfer = heatTransferPtr();
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
const volScalarField& alpha = phase;
const volScalarField& rho = phase.rho();
const volVectorField& U = phase.U();
@ -35,9 +36,9 @@
fluid.correctThermo();
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
Info<< phase.name() << " min/max T "
<< min(phase.thermo().T()).value()

5
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/YEqns.H
View File

@ -5,9 +5,10 @@
phaseSystem::massTransferTable&
massTransfer(massTransferPtr());
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
PtrList<volScalarField>& Y = phase.Y();
const volScalarField& alpha = phase;
const volScalarField& rho = phase.rho();

67
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C
View File

@ -59,9 +59,9 @@ void Foam::multiphaseSystem::calcAlphas()
scalar level = 0.0;
alphas_ == 0.0;
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), i)
{
alphas_ += level*iter();
alphas_ += level*phases()[i];
level += 1.0;
}
@ -72,11 +72,9 @@ void Foam::multiphaseSystem::calcAlphas()
void Foam::multiphaseSystem::solveAlphas()
{
PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = phases()[phasei];
volScalarField& alpha1 = phase;
phase.alphaPhi() =
@ -99,9 +97,9 @@ void Foam::multiphaseSystem::solveAlphas()
surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
forAllIter(PtrDictionary<phaseModel>, phases(), iter2)
forAll(phases(), phasej)
{
phaseModel& phase2 = iter2();
phaseModel& phase2 = phases()[phasej];
volScalarField& alpha2 = phase2;
if (&phase2 == &phase) continue;
@ -170,8 +168,6 @@ void Foam::multiphaseSystem::solveAlphas()
0,
true
);
phasei++;
}
MULES::limitSum(alphaPhiCorrs);
@ -191,10 +187,9 @@ void Foam::multiphaseSystem::solveAlphas()
volScalarField divU(fvc::div(fvc::absolute(phi_, phases().first().U())));
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = phases()[phasei];
volScalarField& alpha = phase;
surfaceScalarField& alphaPhic = alphaPhiCorrs[phasei];
@ -245,9 +240,9 @@ void Foam::multiphaseSystem::solveAlphas()
}
}
forAllConstIter(PtrDictionary<phaseModel>, phases(), iter2)
forAll(phases(), phasej)
{
const phaseModel& phase2 = iter2();
const phaseModel& phase2 = phases()[phasej];
const volScalarField& alpha2 = phase2;
if (&phase2 == &phase) continue;
@ -294,8 +289,6 @@ void Foam::multiphaseSystem::solveAlphas()
<< endl;
sumAlpha += phase;
phasei++;
}
Info<< "Phase-sum volume fraction, min, max = "
@ -506,9 +499,9 @@ Foam::multiphaseSystem::multiphaseSystem
1e-8/pow(average(mesh.V()), 1.0/3.0)
)
{
forAllIter(phaseSystem::phaseModelTable, phases(), iter)
forAll(phases(), phasei)
{
volScalarField& alphai = iter();
volScalarField& alphai = phases()[phasei];
mesh.setFluxRequired(alphai.name());
}
}
@ -547,9 +540,9 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
)
);
forAllConstIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasej)
{
const phaseModel& phase2 = iter();
const phaseModel& phase2 = phases()[phasej];
if (&phase2 != &phase1)
{
@ -591,9 +584,13 @@ Foam::multiphaseSystem::nearInterface() const
)
);
forAllConstIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
tnearInt() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
tnearInt() = max
(
tnearInt(),
pos(phases()[phasei] - 0.01)*pos(0.99 - phases()[phasei])
);
}
return tnearInt;
@ -625,10 +622,9 @@ void Foam::multiphaseSystem::solve()
PtrList<volScalarField> alpha0s(phases().size());
PtrList<surfaceScalarField> phiSums(phases().size());
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = phases()[phasei];
volScalarField& alpha = phase;
alpha0s.set
@ -652,8 +648,6 @@ void Foam::multiphaseSystem::solve()
dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
)
);
phasei++;
}
for
@ -668,18 +662,15 @@ void Foam::multiphaseSystem::solve()
{
solveAlphas();
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
phiSums[phasei] += iter().phi();
phasei++;
phiSums[phasei] += phases()[phasei].phi();
}
}
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = phases()[phasei];
volScalarField& alpha = phase;
phase.phi() = phiSums[phasei]/nAlphaSubCycles;
@ -691,8 +682,6 @@ void Foam::multiphaseSystem::solve()
// Reset the old-time field value
alpha.oldTime() = alpha0s[phasei];
alpha.oldTime().timeIndex() = runTime.timeIndex();
phasei++;
}
}
else
@ -700,9 +689,9 @@ void Foam::multiphaseSystem::solve()
solveAlphas();
}
forAllIter(PtrDictionary<phaseModel>, phases(), iter)
forAll(phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = phases()[phasei];
phase.alphaRhoPhi() = fvc::interpolate(phase.rho())*phase.alphaPhi();
}

8
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/UEqns.H
View File

@ -9,10 +9,10 @@ PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
phaseSystem::momentumTransferTable&
momentumTransfer(momentumTransferPtr());
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
const volScalarField& alpha = phase;
const volScalarField& rho = phase.rho();
volVectorField& U = phase.U();
@ -31,7 +31,5 @@ PtrList<fvVectorMatrix> UEqns(fluid.phases().size());
UEqns[phasei].relax();
fvOptions.constrain(UEqns[phasei]);
phasei++;
}
}

77
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqn.H
View File

@ -2,10 +2,9 @@ PtrList<surfaceScalarField> alphafs(fluid.phases().size());
PtrList<volScalarField> rAUs(fluid.phases().size());
PtrList<surfaceScalarField> alpharAUfs(fluid.phases().size());
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
const volScalarField& alpha = phase;
alphafs.set(phasei, fvc::interpolate(alpha).ptr());
@ -33,8 +32,6 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
fvc::interpolate(max(alpha, phase.residualAlpha())*rAUs[phasei])
).ptr()
);
phasei++;
}
// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
@ -42,10 +39,9 @@ PtrList<surfaceScalarField> phiFs(fluid.phases().size());
{
autoPtr<PtrList<volVectorField> > Fs = fluid.Fs();
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
phiFs.set
(
@ -56,8 +52,6 @@ PtrList<surfaceScalarField> phiFs(fluid.phases().size());
(fvc::interpolate(rAUs[phasei]*Fs()[phasei]) & mesh.Sf())
)
);
phasei++;
}
}
@ -69,10 +63,9 @@ while (pimple.correct())
PtrList<volVectorField> HbyAs(fluid.phases().size());
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
const volScalarField& alpha = phase;
// Correct fixed-flux BCs to be consistent with the velocity BCs
@ -95,8 +88,6 @@ while (pimple.correct())
+ max(phase.residualAlpha() - alpha, scalar(0))
*phase.rho()*phase.U().oldTime()/runTime.deltaT()
);
phasei++;
}
// Mean density for buoyancy force and p_rgh -> p
@ -109,10 +100,9 @@ while (pimple.correct())
);
PtrList<surfaceScalarField> phigs(fluid.phases().size());
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
phigs.set
(
@ -126,8 +116,6 @@ while (pimple.correct())
)
).ptr()
);
phasei++;
}
PtrList<surfaceScalarField> phiHbyAs(fluid.phases().size());
@ -146,10 +134,9 @@ while (pimple.correct())
dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
const volScalarField& alpha = phase;
// ddtPhiCorr filter -- only apply in pure(ish) phases
@ -223,8 +210,6 @@ while (pimple.correct())
}
phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
phasei++;
}
MRF.makeRelative(phiHbyA);
@ -242,11 +227,9 @@ while (pimple.correct())
dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
rAUf += alphafs[phasei]*alpharAUfs[phasei];
phasei++;
}
rAUf = mag(rAUf);
@ -255,12 +238,10 @@ while (pimple.correct())
{
surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
phib = 0;
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField();
phasei++;
}
setSnGrad<fixedFluxPressureFvPatchScalarField>
@ -273,10 +254,9 @@ while (pimple.correct())
}
PtrList<fvScalarMatrix> pEqnComps(fluid.phases().size());
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
if (phase.compressible())
{
@ -338,8 +318,6 @@ while (pimple.correct())
);
}
}
phasei++;
}
// Cache p prior to solve for density update
@ -358,17 +336,14 @@ while (pimple.correct())
{
fvScalarMatrix pEqn(pEqnIncomp);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
if (phase.compressible())
{
pEqn += pEqnComps[phasei];
}
phasei++;
}
solve
@ -385,10 +360,9 @@ while (pimple.correct())
surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
@ -397,8 +371,6 @@ while (pimple.correct())
{
phase.divU(-pEqnComps[phasei] & p_rgh);
}
phasei++;
}
// Optionally relax pressure for velocity correction
@ -406,10 +378,9 @@ while (pimple.correct())
mSfGradp = pEqnIncomp.flux()/rAUf;
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
phase.U() =
HbyAs[phasei]
@ -421,8 +392,6 @@ while (pimple.correct())
);
phase.U().correctBoundaryConditions();
fvOptions.correct(phase.U());
phasei++;
}
}
}
@ -434,9 +403,9 @@ while (pimple.correct())
p_rgh = p - rho*gh;
// Update densities from change in p_rgh
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
phaseModel& phase = iter();
phaseModel& phase = fluid.phases()[phasei];
phase.rho()() += phase.thermo().psi()*(p_rgh - p_rgh_0);
}

285
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/pU/pEqnSave.H
View File

@ -1,285 +0,0 @@
// --- Pressure corrector loop
while (pimple.correct())
{
// rho1 = rho10 + psi1*p_rgh;
// rho2 = rho20 + psi2*p_rgh;
// tmp<fvScalarMatrix> pEqnComp1;
// tmp<fvScalarMatrix> pEqnComp2;
// //if (transonic)
// //{
// //}
// //else
// {
// surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
// surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
// pEqnComp1 =
// fvc::ddt(rho1) + psi1*correction(fvm::ddt(p_rgh))
// + fvc::div(phid1, p_rgh)
// - fvc::Sp(fvc::div(phid1), p_rgh);
// pEqnComp2 =
// fvc::ddt(rho2) + psi2*correction(fvm::ddt(p_rgh))
// + fvc::div(phid2, p_rgh)
// - fvc::Sp(fvc::div(phid2), p_rgh);
// }
PtrList<surfaceScalarField> alphafs(fluid.phases().size());
PtrList<volVectorField> HbyAs(fluid.phases().size());
PtrList<surfaceScalarField> phiHbyAs(fluid.phases().size());
PtrList<volScalarField> rAUs(fluid.phases().size());
PtrList<surfaceScalarField> rAlphaAUfs(fluid.phases().size());
int phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
MRF.makeAbsolute(phase.phi().oldTime());
MRF.makeAbsolute(phase.phi());
HbyAs.set(phasei, new volVectorField(phase.U()));
phiHbyAs.set(phasei, new surfaceScalarField(1.0*phase.phi()));
phasei++;
}
surfaceScalarField phiHbyA
(
IOobject
(
"phiHbyA",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
);
volScalarField rho("rho", fluid.rho());
surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
const volScalarField& alpha = phase;
alphafs.set(phasei, fvc::interpolate(alpha).ptr());
alphafs[phasei].rename("hmm" + alpha.name());
volScalarField dragCoeffi
(
IOobject
(
"dragCoeffi",
runTime.timeName(),
mesh
),
fluid.Kd(phase),
zeroGradientFvPatchScalarField::typeName
);
dragCoeffi.correctBoundaryConditions();
rAUs.set(phasei, (1.0/(UEqns[phasei].A() + dragCoeffi)).ptr());
rAlphaAUfs.set
(
phasei,
(
alphafs[phasei]*fvc::interpolate(phase.rho())
/fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
).ptr()
);
HbyAs[phasei] = rAUs[phasei]*UEqns[phasei].H();
phiHbyAs[phasei] =
(
(fvc::interpolate(HbyAs[phasei]) & mesh.Sf())
+ rAlphaAUfs[phasei]*fvc::ddtCorr(phase.U(), phase.phi())
);
MRF.makeRelative(phiHbyAs[phasei]);
MRF.makeRelative(phase.phi().oldTime());
MRF.makeRelative(phase.phi());
phiHbyAs[phasei] +=
rAlphaAUfs[phasei]
*(
fluid.surfaceTension(phase)*mesh.magSf()
+ (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
- ghSnGradRho
)/fvc::interpolate(phase.rho());
forAllConstIter
(
phaseSystem::KdTable,
fluid.Kds(),
KdIter
)
{
const volScalarField& K(*KdIter());
const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
const phaseModel* phase1 = &pair.phase1();
const phaseModel* phase2 = &pair.phase2();
forAllConstIter(phasePair, pair, iter)
{
if (phase1 == &phase)
{
phiHbyAs[phasei] +=
fvc::interpolate(K)
/fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
*phase2->phi();
HbyAs[phasei] += rAUs[phasei]*K*phase2->U();
}
Swap(phase1, phase2);
}
}
phiHbyA += alphafs[phasei]*phiHbyAs[phasei];
phasei++;
}
surfaceScalarField rAUf
(
IOobject
(
"rAUf",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
);
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
rAUf += mag(alphafs[phasei]*rAlphaAUfs[phasei])
/fvc::interpolate(phase.rho());
phasei++;
}
// Update the fixedFluxPressure BCs to ensure flux consistency
{
surfaceScalarField::GeometricBoundaryField phib(phi.boundaryField());
phib = 0;
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
phib +=
alphafs[phasei].boundaryField()
*(mesh.Sf().boundaryField() & phase.U().boundaryField());
phasei++;
}
setSnGrad<fixedFluxPressureFvPatchScalarField>
(
p_rgh.boundaryField(),
(
phiHbyA.boundaryField() - MRF.relative(phib)
)/(mesh.magSf().boundaryField()*rAUf.boundaryField())
);
}
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqnIncomp
(
fvc::div(phiHbyA)
- fvm::laplacian(rAUf, p_rgh)
);
pEqnIncomp.setReference(pRefCell, pRefValue);
solve
(
// (
// (alpha1/rho1)*pEqnComp1()
// + (alpha2/rho2)*pEqnComp2()
// ) +
pEqnIncomp,
mesh.solver(p_rgh.select(pimple.finalInnerIter()))
);
if (pimple.finalNonOrthogonalIter())
{
surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
phasei = 0;
phi = dimensionedScalar("phi", phi.dimensions(), 0);
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
phase.phi() =
phiHbyAs[phasei]
+ rAlphaAUfs[phasei]*mSfGradp/fvc::interpolate(phase.rho());
phi +=
alphafs[phasei]*phiHbyAs[phasei]
+ mag(alphafs[phasei]*rAlphaAUfs[phasei])
*mSfGradp/fvc::interpolate(phase.rho());
phasei++;
}
// dgdt =
// (
// pos(alpha2)*(pEqnComp2 & p)/rho2
// - pos(alpha1)*(pEqnComp1 & p)/rho1
// );
p_rgh.relax();
p = p_rgh + rho*gh;
mSfGradp = pEqnIncomp.flux()/rAUf;
phasei = 0;
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
{
phaseModel& phase = iter();
phase.U() =
HbyAs[phasei]
+ fvc::reconstruct
(
rAlphaAUfs[phasei]
*(
fvc::interpolate(phase.rho() - rho)
*(g & mesh.Sf())
- ghSnGradRho
+ mSfGradp
)
)/phase.rho();
phase.U().correctBoundaryConditions();
phasei++;
}
}
}
//p = max(p, pMin);
// rho1 = rho10 + psi1*p_rgh;
// rho2 = rho20 + psi2*p_rgh;
// Dp1Dt = fvc::DDt(phi1, p_rgh);
// Dp2Dt = fvc::DDt(phi2, p_rgh);
}

4
applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/setRDeltaT.H
View File

@ -12,9 +12,9 @@
surfaceScalarField maxPhi("maxPhi", phi);
forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
forAll(fluid.phases(), phasei)
{
maxPhi = max(maxPhi, mag(iter().phi()));
maxPhi = max(maxPhi, mag(fluid.phases()[phasei].phi()));
}
// Set the reciprocal time-step from the local Courant number

4
applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
View File

@ -57,8 +57,8 @@ Foam::twoPhaseSystem::twoPhaseSystem
)
:
phaseSystem(mesh),
phase1_(phaseModels_.first()),
phase2_(phaseModels_.last())
phase1_(phaseModels_[0]),
phase2_(phaseModels_[1])
{
phase2_.volScalarField::operator=(scalar(1) - phase1_);