STY: General clean up

This commit is contained in:
sergio
2020-02-19 10:40:56 -08:00
parent 83c06f1ace
commit 514751dcf7
31 changed files with 192 additions and 836 deletions

View File

@ -1,7 +1,7 @@
{ {
radiation->correct(); radiation->correct();
rhoCp = rho*fluid.Cp(); rhoCp = rho*fluid.Cp();
const surfaceScalarField rhoCpPhi(fvc::interpolate(fluid.Cp())*rhoPhi); const surfaceScalarField rhoCpPhi(fvc::interpolate(fluid.Cp())*rhoPhi);
const volScalarField kappaEff const volScalarField kappaEff
@ -31,7 +31,7 @@
fvOptions.correct(T); fvOptions.correct(T);
fluid.correct(); fluid.correct();
Info<< "min/max(T) = " Info<< "min/max(T) = "
<< min(T).value() << ", " << max(T).value() << endl; << min(T).value() << ", " << max(T).value() << endl;
} }

View File

@ -59,14 +59,14 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Kexp
( (
min(max(this->pair().from(), scalar(0)), scalar(1)) min(max(this->pair().from(), scalar(0)), scalar(1))
); );
const volScalarField coeff const volScalarField coeff
( (
C_*from*this->pair().from().rho()*pos(from - alphaMin_) C_*from*this->pair().from().rho()*pos(from - alphaMin_)
*(refValue - Tactivate_) *(refValue - Tactivate_)
/Tactivate_ /Tactivate_
); );
if (sign(C_.value()) > 0) if (sign(C_.value()) > 0)
{ {
return return
@ -99,7 +99,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSp
( (
min(max(this->pair().from(), scalar(0)), scalar(1)) min(max(this->pair().from(), scalar(0)), scalar(1))
); );
const volScalarField coeff const volScalarField coeff
( (
C_*from*this->pair().from().rho()*pos(from - alphaMin_) C_*from*this->pair().from().rho()*pos(from - alphaMin_)
@ -142,7 +142,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::KSu
( (
min(max(this->pair().from(), scalar(0)), scalar(1)) min(max(this->pair().from(), scalar(0)), scalar(1))
); );
const volScalarField coeff const volScalarField coeff
( (
C_*from*this->pair().from().rho()*pos(from - alphaMin_) C_*from*this->pair().from().rho()*pos(from - alphaMin_)
@ -179,7 +179,7 @@ Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::Tactivate() const
template<class Thermo, class OtherThermo> template<class Thermo, class OtherThermo>
bool bool
Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::includeDivU() Foam::meltingEvaporationModels::Lee<Thermo, OtherThermo>::includeDivU()
{ {
return true; return true;

View File

@ -148,14 +148,14 @@ public:
( (
const volScalarField& field const volScalarField& field
); );
//- Implicit mass transfer coefficient //- Implicit mass transfer coefficient
virtual tmp<volScalarField> KSp virtual tmp<volScalarField> KSp
( (
label modelVariable, label modelVariable,
const volScalarField& field const volScalarField& field
); );
//- Explicit mass transfer coefficient //- Explicit mass transfer coefficient
virtual tmp<volScalarField> KSu virtual tmp<volScalarField> KSu
( (
@ -165,9 +165,9 @@ public:
//- Return T transition between phases //- Return T transition between phases
virtual const dimensionedScalar& Tactivate() const; virtual const dimensionedScalar& Tactivate() const;
//- Adds and substract alpha*div(U) as a source term //- Adds and substract alpha*div(U) as a source term
// for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
virtual bool includeDivU(); virtual bool includeDivU();
}; };

View File

@ -176,25 +176,24 @@ public:
( (
const volScalarField& field const volScalarField& field
) = 0; ) = 0;
//- Implicit mass transfer //- Implicit mass transfer
virtual tmp<volScalarField> KSp virtual tmp<volScalarField> KSp
( (
label modelVariable, label modelVariable,
const volScalarField& field const volScalarField& field
) = 0; ) = 0;
//- Explicit mass transfer //- Explicit mass transfer
virtual tmp<volScalarField> KSu virtual tmp<volScalarField> KSu
( (
label modelVariable, label modelVariable,
const volScalarField& field const volScalarField& field
) = 0; ) = 0;
//- Reference value //- Reference value
virtual const dimensionedScalar& Tactivate() const = 0; virtual const dimensionedScalar& Tactivate() const = 0;
//- Adds and substract alpha*div(U) as a source term //- Adds and substract alpha*div(U) as a source term
// for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
virtual bool includeDivU(); virtual bool includeDivU();

View File

@ -187,14 +187,14 @@ public:
( (
const volScalarField& field const volScalarField& field
); );
//- Implicit mass transfer coefficient //- Implicit mass transfer coefficient
virtual tmp<volScalarField> KSp virtual tmp<volScalarField> KSp
( (
label modelVariable, label modelVariable,
const volScalarField& field const volScalarField& field
); );
//- Explicit mass transfer coefficient //- Explicit mass transfer coefficient
virtual tmp<volScalarField> KSu virtual tmp<volScalarField> KSu
( (

View File

@ -35,7 +35,7 @@
- fvm::laplacian(rAUf, p_rgh) - fvm::laplacian(rAUf, p_rgh)
+ fluid.volTransfer(p_rgh) + fluid.volTransfer(p_rgh)
); );
p_rghEqn.setReference(pRefCell, pRefValue); p_rghEqn.setReference(pRefCell, pRefValue);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

View File

@ -219,7 +219,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
) )
); );
volScalarField& dmdtNetki = tdmdtNetki.ref(); volScalarField& dmdtNetki = tdmdtNetki.ref();
tmp<volScalarField> tSp tmp<volScalarField> tSp
( (
new volScalarField new volScalarField
@ -235,7 +235,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
) )
); );
volScalarField& Sp = tSp.ref(); volScalarField& Sp = tSp.ref();
tmp<volScalarField> tSu tmp<volScalarField> tSu
( (
new volScalarField new volScalarField
@ -259,23 +259,23 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
massTransferModels_[keyik]; massTransferModels_[keyik];
dmdtNetki -= *dmdt_[keyik]; dmdtNetki -= *dmdt_[keyik];
tmp<volScalarField> KSp = tmp<volScalarField> KSp =
interfacePtr->KSp(interfaceCompositionModel::T, T); interfacePtr->KSp(interfaceCompositionModel::T, T);
if (KSp.valid()) if (KSp.valid())
{ {
Sp -= KSp.ref(); Sp -= KSp.ref();
} }
tmp<volScalarField> KSu = tmp<volScalarField> KSu =
interfacePtr->KSu(interfaceCompositionModel::T, T); interfacePtr->KSu(interfaceCompositionModel::T, T);
if (KSu.valid()) if (KSu.valid())
{ {
Su -= KSu.ref(); Su -= KSu.ref();
} }
// If linearization is not provided used full explicit // If linearization is not provided used full explicit
if (!KSp.valid() && !KSu.valid()) if (!KSp.valid() && !KSu.valid())
{ {
@ -291,23 +291,23 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::heatTransfer
dmdtNetki += *dmdt_[keyki]; dmdtNetki += *dmdt_[keyki];
tmp<volScalarField> KSp = tmp<volScalarField> KSp =
interfacePtr->KSp(interfaceCompositionModel::T, T); interfacePtr->KSp(interfaceCompositionModel::T, T);
if (KSp.valid()) if (KSp.valid())
{ {
Sp += KSp.ref(); Sp += KSp.ref();
} }
tmp<volScalarField> KSu = tmp<volScalarField> KSu =
interfacePtr->KSu(interfaceCompositionModel::T, T); interfacePtr->KSu(interfaceCompositionModel::T, T);
if (KSu.valid()) if (KSu.valid())
{ {
Su += KSu.ref(); Su += KSu.ref();
} }
// If linearization is not provided used full explicit // If linearization is not provided used full explicit
if (!KSp.valid() && !KSu.valid()) if (!KSp.valid() && !KSu.valid())
{ {
@ -337,9 +337,9 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
( (
new fvScalarMatrix(p, dimVolume/dimTime) new fvScalarMatrix(p, dimVolume/dimTime)
); );
fvScalarMatrix& eqn = tEqnPtr.ref(); fvScalarMatrix& eqn = tEqnPtr.ref();
tmp<volScalarField> tSp tmp<volScalarField> tSp
( (
new volScalarField new volScalarField
@ -355,7 +355,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
) )
); );
volScalarField& Sp = tSp.ref(); volScalarField& Sp = tSp.ref();
tmp<volScalarField> tSu tmp<volScalarField> tSu
( (
new volScalarField new volScalarField
@ -371,7 +371,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
) )
); );
volScalarField& Su = tSu.ref(); volScalarField& Su = tSu.ref();
forAllConstIters(this->totalPhasePairs(), iter) forAllConstIters(this->totalPhasePairs(), iter)
{ {
const phasePair& pair = iter()(); const phasePair& pair = iter()();
@ -385,72 +385,72 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
phase2.name(), phase2.name(),
true true
); );
if (massTransferModels_.found(key12)) if (massTransferModels_.found(key12))
{ {
autoPtr<interfaceCompositionModel>& interfacePtr = autoPtr<interfaceCompositionModel>& interfacePtr =
massTransferModels_[key12]; massTransferModels_[key12];
tmp<volScalarField> KSp = tmp<volScalarField> KSp =
interfacePtr->KSp(interfaceCompositionModel::P, p); interfacePtr->KSp(interfaceCompositionModel::P, p);
if (KSp.valid()) if (KSp.valid())
{ {
Sp += KSp.ref(); Sp += KSp.ref();
} }
tmp<volScalarField> KSu = tmp<volScalarField> KSu =
interfacePtr->KSu(interfaceCompositionModel::P, p); interfacePtr->KSu(interfaceCompositionModel::P, p);
if (KSu.valid()) if (KSu.valid())
{ {
Su += KSu.ref(); Su += KSu.ref();
} }
// If linearization is not provided used full explicit // If linearization is not provided used full explicit
if (!KSp.valid() && !KSu.valid()) if (!KSp.valid() && !KSu.valid())
{ {
Su -= Su -=
*dmdt_[key12] *dmdt_[key12]
*( *(
- this->coeffs(phase1.name()) - this->coeffs(phase1.name())
+ this->coeffs(phase2.name()) + this->coeffs(phase2.name())
); );
} }
} }
const phasePairKey key21 const phasePairKey key21
( (
phase2.name(), phase2.name(),
phase1.name(), phase1.name(),
true true
); );
if (massTransferModels_.found(key21)) if (massTransferModels_.found(key21))
{ {
autoPtr<interfaceCompositionModel>& interfacePtr = autoPtr<interfaceCompositionModel>& interfacePtr =
massTransferModels_[key21]; massTransferModels_[key21];
tmp<volScalarField> KSp = tmp<volScalarField> KSp =
interfacePtr->KSp(interfaceCompositionModel::P, p); interfacePtr->KSp(interfaceCompositionModel::P, p);
if (KSp.valid()) if (KSp.valid())
{ {
Sp += KSp.ref(); Sp += KSp.ref();
} }
tmp<volScalarField> KSu = tmp<volScalarField> KSu =
interfacePtr->KSu(interfaceCompositionModel::P, p); interfacePtr->KSu(interfaceCompositionModel::P, p);
if (KSu.valid()) if (KSu.valid())
{ {
Su += KSu.ref(); Su += KSu.ref();
} }
// If linearization is not provided used full explicit // If linearization is not provided used full explicit
if (!KSp.valid() && !KSu.valid()) if (!KSp.valid() && !KSu.valid())
{ {
Su += Su +=
*dmdt_[key21] *dmdt_[key21]
*( *(
- this->coeffs(phase1.name()) - this->coeffs(phase1.name())
@ -460,7 +460,7 @@ Foam::MassTransferPhaseSystem<BasePhaseSystem>::volTransfer
} }
} }
eqn += fvm::Sp(Sp, p) + Su; eqn += fvm::Sp(Sp, p) + Su;
return tEqnPtr; return tEqnPtr;
} }
@ -489,18 +489,18 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::correctMassSources
// Phase k to phase i // Phase k to phase i
const phasePairKey keyki(phasek.name(), phasei.name(), true); const phasePairKey keyki(phasek.name(), phasei.name(), true);
if (massTransferModels_.found(keyik)) if (massTransferModels_.found(keyik))
{ {
autoPtr<interfaceCompositionModel>& interfacePtr = autoPtr<interfaceCompositionModel>& interfacePtr =
massTransferModels_[keyik]; massTransferModels_[keyik];
tmp<volScalarField> Kexp = interfacePtr->Kexp(T); tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
*dmdt_[keyik] = Kexp.ref(); *dmdt_[keyik] = Kexp.ref();
} }
if (massTransferModels_.found(keyki)) if (massTransferModels_.found(keyki))
{ {
autoPtr<interfaceCompositionModel>& interfacePtr = autoPtr<interfaceCompositionModel>& interfacePtr =
@ -508,7 +508,7 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::correctMassSources
// Explicit temperature mass transfer rate // Explicit temperature mass transfer rate
const tmp<volScalarField> Kexp = interfacePtr->Kexp(T); const tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
*dmdt_[keyki] = Kexp.ref(); *dmdt_[keyki] = Kexp.ref();
} }
} }
@ -527,7 +527,7 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
// This term adds and substract alpha*div(U) as a source term // This term adds and substract alpha*div(U) as a source term
// for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2) // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
bool includeDivU(true); bool includeDivU(true);
forAllConstIters(this->totalPhasePairs(), iter) forAllConstIters(this->totalPhasePairs(), iter)
{ {
const phasePair& pair = iter()(); const phasePair& pair = iter()();
@ -551,26 +551,26 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
phase2.name(), phase2.name(),
true true
); );
tmp<volScalarField> tdmdt12(this->dmdt(key12)); tmp<volScalarField> tdmdt12(this->dmdt(key12));
volScalarField& dmdt12 = tdmdt12.ref(); volScalarField& dmdt12 = tdmdt12.ref();
if (massTransferModels_.found(key12)) if (massTransferModels_.found(key12))
{ {
autoPtr<interfaceCompositionModel>& interfacePtr = autoPtr<interfaceCompositionModel>& interfacePtr =
massTransferModels_[key12]; massTransferModels_[key12];
tmp<volScalarField> KSu = tmp<volScalarField> KSu =
interfacePtr->KSu(interfaceCompositionModel::alpha, phase1); interfacePtr->KSu(interfaceCompositionModel::alpha, phase1);
if (KSu.valid()) if (KSu.valid())
{ {
dmdt12 = KSu.ref(); dmdt12 = KSu.ref();
} }
includeDivU = interfacePtr->includeDivU(); includeDivU = interfacePtr->includeDivU();
} }
// Phase 2 to phase 1 // Phase 2 to phase 1
const phasePairKey key21 const phasePairKey key21
@ -582,20 +582,20 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
tmp<volScalarField> tdmdt21(this->dmdt(key21)); tmp<volScalarField> tdmdt21(this->dmdt(key21));
volScalarField& dmdt21 = tdmdt21.ref(); volScalarField& dmdt21 = tdmdt21.ref();
if (massTransferModels_.found(key21)) if (massTransferModels_.found(key21))
{ {
autoPtr<interfaceCompositionModel>& interfacePtr = autoPtr<interfaceCompositionModel>& interfacePtr =
massTransferModels_[key21]; massTransferModels_[key21];
tmp<volScalarField> KSu = tmp<volScalarField> KSu =
interfacePtr->KSu(interfaceCompositionModel::alpha, phase2); interfacePtr->KSu(interfaceCompositionModel::alpha, phase2);
if (KSu.valid()) if (KSu.valid())
{ {
dmdt21 = KSu.ref(); dmdt21 = KSu.ref();
} }
includeDivU = interfacePtr->includeDivU(); includeDivU = interfacePtr->includeDivU();
} }
@ -612,16 +612,16 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
const volScalarField coeffs12(coeffs1 - coeffs2); const volScalarField coeffs12(coeffs1 - coeffs2);
const surfaceScalarField& phi = this->phi(); const surfaceScalarField& phi = this->phi();
if (includeDivU) if (includeDivU)
{ {
SuPhase1 += SuPhase1 +=
fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1)); fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1));
SuPhase2 += SuPhase2 +=
fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1)); fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1));
} }
// NOTE: dmdtNet is distributed in terms = // NOTE: dmdtNet is distributed in terms =
// Source for phase 1 = // Source for phase 1 =
// dmdtNet/rho1 // dmdtNet/rho1
@ -715,12 +715,12 @@ void Foam::MassTransferPhaseSystem<BasePhaseSystem>::alphaTransfer
} }
// Update ddtAlphaMax // Update ddtAlphaMax
this->ddtAlphaMax_ = this->ddtAlphaMax_ =
max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)())); max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
} }
} }
template<class BasePhaseSystem> template<class BasePhaseSystem>
void Foam::MassTransferPhaseSystem<BasePhaseSystem>::massSpeciesTransfer void Foam::MassTransferPhaseSystem<BasePhaseSystem>::massSpeciesTransfer
( (

View File

@ -67,8 +67,8 @@ public:
phasePairKey::hash phasePairKey::hash
> >
massTransferModelTable; massTransferModelTable;
typedef HashTable<volScalarField::Internal> SuSpTable; typedef HashTable<volScalarField::Internal> SuSpTable;
protected: protected:
@ -128,13 +128,13 @@ public:
//- Return the heat transfer matrix //- Return the heat transfer matrix
virtual tmp<fvScalarMatrix> heatTransfer(const volScalarField& T); virtual tmp<fvScalarMatrix> heatTransfer(const volScalarField& T);
//- Return the volumetric rate transfer matrix //- Return the volumetric rate transfer matrix
virtual tmp<fvScalarMatrix> volTransfer(const volScalarField& p); virtual tmp<fvScalarMatrix> volTransfer(const volScalarField& p);
//- Correct/calculates mass sources dmdt for phases //- Correct/calculates mass sources dmdt for phases
virtual void correctMassSources(const volScalarField& T); virtual void correctMassSources(const volScalarField& T);
//- Calculate mass transfer for alpha's //- Calculate mass transfer for alpha's
virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp); virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp);

View File

@ -120,7 +120,7 @@ Foam::multiphaseSystem::multiphaseSystem
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::multiphaseSystem::calculateSuSp() void Foam::multiphaseSystem::calculateSuSp()
{ {
this->alphaTransfer(Su_, Sp_); this->alphaTransfer(Su_, Sp_);
} }
@ -129,9 +129,9 @@ void Foam::multiphaseSystem::solve()
{ {
const dictionary& alphaControls = mesh_.solverDict("alpha"); const dictionary& alphaControls = mesh_.solverDict("alpha");
label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles")); label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
volScalarField& alpha = phases_.first(); volScalarField& alpha = phases_.first();
if (nAlphaSubCycles > 1) if (nAlphaSubCycles > 1)
{ {
surfaceScalarField rhoPhiSum surfaceScalarField rhoPhiSum
@ -164,16 +164,16 @@ void Foam::multiphaseSystem::solve()
{ {
solveAlphas(); solveAlphas();
} }
} }
void Foam::multiphaseSystem::solveAlphas() void Foam::multiphaseSystem::solveAlphas()
{ {
mesh_.solverDict("alpha").readEntry("cAlphas", cAlphas_);
const dictionary& alphaControls = mesh_.solverDict("alpha"); const dictionary& alphaControls = mesh_.solverDict("alpha");
alphaControls.readEntry("cAlphas", cAlphas_);
label nAlphaCorr(alphaControls.get<label>("nAlphaCorr")); label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size()); PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
const surfaceScalarField& phi = this->phi(); const surfaceScalarField& phi = this->phi();
@ -361,7 +361,7 @@ void Foam::multiphaseSystem::solveAlphas()
alpha1Eqn.solve(); alpha1Eqn.solve();
phiAlpha += alpha1Eqn.flux(); phiAlpha += alpha1Eqn.flux();
MULES::explicitSolve MULES::explicitSolve
( (
geometricOneField(), geometricOneField(),
@ -375,7 +375,7 @@ void Foam::multiphaseSystem::solveAlphas()
); );
phase.alphaPhi() = phiAlpha; phase.alphaPhi() = phiAlpha;
++phasei; ++phasei;
} }

View File

@ -89,7 +89,7 @@ protected:
//- Calculate Sp and Su //- Calculate Sp and Su
void calculateSuSp(); void calculateSuSp();
//- Solve alphas //- Solve alphas
void solveAlphas(); void solveAlphas();

View File

@ -1,642 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2017-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "multiphaseSystem.H"
#include "fixedValueFvsPatchFields.H"
#include "Time.H"
#include "subCycle.H"
#include "fvcMeshPhi.H"
#include "surfaceInterpolate.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
#include "fvcDiv.H"
#include "fvcDdt.H"
#include "fvcFlux.H"
#include "fvmDdt.H"
#include "fvcAverage.H"
#include "fvMatrix.H"
#include "fvmSup.H"
#include "CMULES.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(multiphaseSystem, 0);
defineRunTimeSelectionTable(multiphaseSystem, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::multiphaseSystem::multiphaseSystem
(
const fvMesh& mesh
)
:
phaseSystem(mesh),
cAlphas_(),
ddtAlphaMax_(0.0),
limitedPhiAlphas_(phaseModels_.size()),
Su_(phaseModels_.size()),
Sp_(phaseModels_.size())
{
label phasei = 0;
phases_.setSize(phaseModels_.size());
forAllIters(phaseModels_, iter)
{
phaseModel& pm = iter()();
phases_.set(phasei++, &pm);
}
mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
// Initiate Su and Sp
forAllConstIters(phaseModels_, iter)
{
const phaseModel& pm = iter()();
Su_.insert
(
pm.name(),
volScalarField::Internal
(
IOobject
(
"Su" + pm.name(),
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(dimless/dimTime, Zero)
)
);
Sp_.insert
(
pm.name(),
volScalarField::Internal
(
IOobject
(
"Sp" + pm.name(),
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(dimless/dimTime, Zero)
)
);
}
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::multiphaseSystem::calculateSuSp()
{
forAllConstIters(totalPhasePairs_, iter)
{
const phasePair& pair = iter()();
const phaseModel& phase1 = pair.phase1();
const phaseModel& phase2 = pair.phase2();
const volScalarField& alpha1 = pair.phase1();
const volScalarField& alpha2 = pair.phase2();
tmp<volScalarField> tCoeffs1 = this->coeffs(phase1.name());
const volScalarField& coeffs1 = tCoeffs1();
tmp<volScalarField> tCoeffs2 = this->coeffs(phase2.name());
const volScalarField& coeffs2 = tCoeffs2();
// Phase 1 to phase 2
const phasePairKey key12
(
phase1.name(),
phase2.name(),
true
);
tmp<volScalarField> tdmdt12(this->dmdt(key12));
const volScalarField& dmdt12 = tdmdt12();
// Phase 2 to phase 1
const phasePairKey key21
(
phase2.name(),
phase1.name(),
true
);
tmp<volScalarField> tdmdt21(this->dmdt(key21));
const volScalarField& dmdt21 = tdmdt21();
volScalarField::Internal& SpPhase1 = Sp_[phase1.name()];
volScalarField::Internal& SuPhase1 = Su_[phase1.name()];
volScalarField::Internal& SpPhase2 = Sp_[phase2.name()];
volScalarField::Internal& SuPhase2 = Su_[phase2.name()];
const volScalarField dmdtNet(dmdt21 - dmdt12);
const volScalarField coeffs12(coeffs1 - coeffs2);
// NOTE: dmdtNet is distributed in terms =
// Source for phase 1 =
// dmdtNet/rho1
// - alpha1*dmdtNet(1/rho1 - 1/rho2)
forAll(dmdtNet, celli)
{
scalar dmdt21 = dmdtNet[celli];
scalar coeffs12Cell = coeffs12[celli];
scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0);
// exp.
SuPhase1[celli] += coeffs1[celli]*dmdt21;
if (dmdt21 > 0)
{
if (coeffs12Cell > 0)
{
// imp
SpPhase1[celli] -= dmdt21*coeffs12Cell;
}
else if (coeffs12Cell < 0)
{
// exp
SuPhase1[celli] -=
dmdt21*coeffs12Cell*alpha1Limited;
}
}
else if (dmdt21 < 0)
{
if (coeffs12Cell > 0)
{
// exp
SuPhase1[celli] -=
dmdt21*coeffs12Cell*alpha1Limited;
}
else if (coeffs12Cell < 0)
{
// imp
SpPhase1[celli] -= dmdt21*coeffs12Cell;
}
}
}
forAll(dmdtNet, celli)
{
scalar dmdt12 = -dmdtNet[celli];
scalar coeffs21Cell = -coeffs12[celli];
scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0);
// exp
SuPhase2[celli] += coeffs2[celli]*dmdt12;
if (dmdt12 > 0)
{
if (coeffs21Cell > 0)
{
// imp
SpPhase2[celli] -= dmdt12*coeffs21Cell;
}
else if (coeffs21Cell < 0)
{
// exp
SuPhase2[celli] -=
dmdt12*coeffs21Cell*alpha2Limited;
}
}
else if (dmdt12 < 0)
{
if (coeffs21Cell > 0)
{
// exp
SuPhase2[celli] -=
coeffs21Cell*dmdt12*alpha2Limited;
}
else if (coeffs21Cell < 0)
{
// imp
SpPhase2[celli] -= dmdt12*coeffs21Cell;
}
}
}
// Update ddtAlphaMax
ddtAlphaMax_ =
max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
}
}
void Foam::multiphaseSystem::solve()
{
const fvMesh& mesh = this->mesh();
const dictionary& alphaControls = mesh.solverDict("alpha");
label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
mesh.solverDict("alpha").readEntry("cAlphas", cAlphas_);
PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
const surfaceScalarField& phi = this->phi();
surfaceScalarField phic(mag((phi)/mesh_.magSf()));
// Do not compress interface at non-coupled boundary faces
// (inlets, outlets etc.)
surfaceScalarField::Boundary& phicBf = phic.boundaryFieldRef();
forAll(phic.boundaryField(), patchi)
{
fvsPatchScalarField& phicp = phicBf[patchi];
if (!phicp.coupled())
{
phicp == 0;
}
}
for (int acorr=0; acorr<nAlphaCorr; acorr++)
{
label phasei = 0;
for (phaseModel& phase1 : phases_)
{
const volScalarField& alpha1 = phase1;
phiAlphaCorrs.set
(
phasei,
new surfaceScalarField
(
"phi" + alpha1.name() + "Corr",
fvc::flux
(
phi,
alpha1,
"div(phi," + alpha1.name() + ')'
)
)
);
surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
for (phaseModel& phase2 : phases_)
{
const volScalarField& alpha2 = phase2;
if (&phase2 == &phase1) continue;
const phasePairKey key12(phase1.name(), phase2.name());
if (!cAlphas_.found(key12))
{
FatalErrorInFunction
<< "Phase compression factor (cAlpha) not found for : "
<< key12
<< exit(FatalError);
}
scalar cAlpha = cAlphas_.find(key12)();
phic = min(cAlpha*phic, max(phic));
surfaceScalarField phir(phic*nHatf(alpha1, alpha2));
word phirScheme
(
"div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
);
phiAlphaCorr += fvc::flux
(
-fvc::flux(-phir, alpha2, phirScheme),
alpha1,
phirScheme
);
}
// Ensure that the flux at inflow BCs is preserved
forAll(phiAlphaCorr.boundaryField(), patchi)
{
fvsPatchScalarField& phiAlphaCorrp =
phiAlphaCorr.boundaryFieldRef()[patchi];
if (!phiAlphaCorrp.coupled())
{
const scalarField& phi1p = phi.boundaryField()[patchi];
const scalarField& alpha1p =
alpha1.boundaryField()[patchi];
forAll(phiAlphaCorrp, facei)
{
if (phi1p[facei] < 0)
{
phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei];
}
}
}
}
++phasei;
}
// Set Su and Sp to zero
for (const phaseModel& phase : phases_)
{
Su_[phase.name()] = dimensionedScalar("Su", dimless/dimTime, Zero);
Sp_[phase.name()] = dimensionedScalar("Sp", dimless/dimTime, Zero);
// Add alpha*div(U)
const volScalarField& alpha = phase;
Sp_[phase.name()] +=
fvc::div(phi);//*min(max(alpha, scalar(0)), scalar(1));
}
// Fill Su and Sp
calculateSuSp();
// Limit phiAlphaCorr on each phase
phasei = 0;
for (phaseModel& phase : phases_)
{
volScalarField& alpha1 = phase;
surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
volScalarField::Internal& Su = Su_[phase.name()];
volScalarField::Internal& Sp = Sp_[phase.name()];
MULES::limit
(
1.0/mesh_.time().deltaT().value(),
geometricOneField(),
alpha1,
phi,
phiAlphaCorr,
Sp,
Su,
oneField(),
zeroField(),
true
);
++phasei;
}
MULES::limitSum(phiAlphaCorrs);
volScalarField sumAlpha
(
IOobject
(
"sumAlpha",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(dimless, Zero)
);
phasei = 0;
for (phaseModel& phase : phases_)
{
volScalarField& alpha1 = phase;
const volScalarField::Internal& Su = Su_[phase.name()];
const volScalarField::Internal& Sp = Sp_[phase.name()];
surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
// Add a bounded upwind U-mean flux
phiAlpha += upwind<scalar>(mesh_, phi).flux(alpha1);
// fvScalarMatrix alpha1Eqn
// (
// fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
// + fv::gaussConvectionScheme<scalar>
// (
// mesh,
// phi,
// upwind<scalar>(mesh, phi)
// ).fvmDiv(phi, alpha1)
// ==
// Su + fvm::Sp(Sp, alpha1)
// );
//
// alpha1Eqn.solve();
//phiAlpha += alpha1Eqn.flux();
if (nAlphaSubCycles > 1)
{
for
(
subCycle<volScalarField> alphaSubCycle
(
alpha1,
nAlphaSubCycles
);
!(++alphaSubCycle).end();
)
{
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
phiAlpha,
(alphaSubCycle.index()*Sp)(),
(Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
oneField(),
zeroField()
);
if (alphaSubCycle.index() == 1)
{
phase.alphaPhi() = phiAlpha;
}
else
{
phase.alphaPhi() += phiAlpha;
}
}
phase.alphaPhi() /= nAlphaSubCycles;
}
else
{
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
phiAlpha,
Sp,
Su,
oneField(),
zeroField()
);
phase.alphaPhi() = phiAlpha;
}
++phasei;
}
if (acorr == nAlphaCorr - 1)
{
volScalarField sumAlpha
(
IOobject
(
"sumAlpha",
mesh_.time().timeName(),
mesh_
),
mesh_,
dimensionedScalar(dimless, Zero)
);
// Reset rhoPhi
rhoPhi_ = dimensionedScalar("rhoPhi", dimMass/dimTime, Zero);
for (phaseModel& phase : phases_)
{
volScalarField& alpha1 = phase;
sumAlpha += alpha1;
// Update rhoPhi
rhoPhi_ += fvc::interpolate(phase.rho()) * phase.alphaPhi();
}
Info<< "Phase-sum volume fraction, min, max = "
<< sumAlpha.weightedAverage(mesh_.V()).value()
<< ' ' << min(sumAlpha).value()
<< ' ' << max(sumAlpha).value()
<< endl;
volScalarField sumCorr(1.0 - sumAlpha);
for (phaseModel& phase : phases_)
{
volScalarField& alpha = phase;
//alpha += alpha*sumCorr;
Info<< alpha.name() << " volume fraction = "
<< alpha.weightedAverage(mesh.V()).value()
<< " Min(alpha) = " << min(alpha).value()
<< " Max(alpha) = " << max(alpha).value()
<< endl;
}
}
}
}
const Foam::UPtrList<Foam::phaseModel>& Foam::multiphaseSystem::phases() const
{
return phases_;
}
Foam::UPtrList<Foam::phaseModel>& Foam::multiphaseSystem::phases()
{
return phases_;
}
const Foam::phaseModel& Foam::multiphaseSystem::phase(const label i) const
{
return phases_[i];
}
Foam::phaseModel& Foam::multiphaseSystem::phase(const label i)
{
return phases_[i];
}
Foam::dimensionedScalar Foam::multiphaseSystem::ddtAlphaMax() const
{
return ddtAlphaMax_;
}
Foam::scalar Foam::multiphaseSystem::maxDiffNo() const
{
auto iter = phaseModels_.cbegin();
scalar maxVal = max(iter()->diffNo()).value();
for (++iter; iter != phaseModels_.cend(); ++iter)
{
maxVal = max(maxVal, max(iter()->diffNo()).value());
}
return maxVal * mesh_.time().deltaT().value();
}
const Foam::multiphaseSystem::compressionFluxTable&
Foam::multiphaseSystem::limitedPhiAlphas() const
{
return limitedPhiAlphas_;
}
Foam::multiphaseSystem::SuSpTable& Foam::multiphaseSystem::Su()
{
return Su_;
}
Foam::multiphaseSystem::SuSpTable& Foam::multiphaseSystem::Sp()
{
return Sp_;
}
bool Foam::multiphaseSystem::read()
{
return true;
}
// ************************************************************************* //

View File

@ -281,7 +281,7 @@ Foam::phaseSystem::phaseSystem
// Total phase pair // Total phase pair
generatePairsTable(); generatePairsTable();
// Update mu_ // Update mu_
calcMu(); calcMu();
} }
@ -903,7 +903,7 @@ void Foam::phaseSystem::correct()
{ {
iter()->correct(); iter()->correct();
} }
calcMu(); calcMu();
} }

View File

@ -78,7 +78,7 @@ public:
typedef HashTable<autoPtr<phaseModel>> phaseModelTable; typedef HashTable<autoPtr<phaseModel>> phaseModelTable;
typedef HashTable<volScalarField::Internal> SuSpTable; typedef HashTable<volScalarField::Internal> SuSpTable;
@ -116,7 +116,7 @@ protected:
//- Reference to the mesh //- Reference to the mesh
const fvMesh& mesh_; const fvMesh& mesh_;
//- Dynamic viscocity //- Dynamic viscocity
volScalarField mu_; volScalarField mu_;
@ -510,13 +510,13 @@ public:
( (
const volScalarField& T const volScalarField& T
) = 0; ) = 0;
//- Return the volumetric rate transfer matrix //- Return the volumetric rate transfer matrix
virtual tmp<fvScalarMatrix> volTransfer virtual tmp<fvScalarMatrix> volTransfer
( (
const volScalarField& p const volScalarField& p
) = 0; ) = 0;
//- Calculate mass transfer for alpha's //- Calculate mass transfer for alpha's
virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp) = 0; virtual void alphaTransfer(SuSpTable& Su, SuSpTable& Sp) = 0;
@ -537,7 +537,7 @@ public:
//- Correct the mixture thermos //- Correct the mixture thermos
virtual void correct(); virtual void correct();
//- Correct mass sources //- Correct mass sources
virtual void correctMassSources(const volScalarField& T) = 0; virtual void correctMassSources(const volScalarField& T) = 0;

View File

@ -1,14 +1,14 @@
{ {
tmp<volScalarField> tcp(thermo->Cp()); tmp<volScalarField> tcp(thermo->Cp());
const volScalarField& cp = tcp(); const volScalarField& cp = tcp();
const dimensionedScalar Cp1 = thermo->Cp1(); const dimensionedScalar Cp1 = thermo->Cp1();
const dimensionedScalar Cp2 = thermo->Cp2(); const dimensionedScalar Cp2 = thermo->Cp2();
rhoCp = rho*cp; rhoCp = rho*cp;
kappaEff = thermo->kappa() + rho*cp*turbulence->nut()/Prt; kappaEff = thermo->kappa() + rho*cp*turbulence->nut()/Prt;
const surfaceScalarField rhoCpPhi const surfaceScalarField rhoCpPhi
( (
"rhoCpPhi", "rhoCpPhi",

View File

@ -84,7 +84,7 @@ int main(int argc, char *argv[])
#include "createAlphaFluxes.H" #include "createAlphaFluxes.H"
#include "initCorrectPhi.H" #include "initCorrectPhi.H"
#include "createUfIfPresent.H" #include "createUfIfPresent.H"
#include "CourantNo.H" #include "CourantNo.H"
#include "setInitialDeltaT.H" #include "setInitialDeltaT.H"
@ -99,13 +99,13 @@ int main(int argc, char *argv[])
while (runTime.run()) while (runTime.run())
{ {
#include "readDyMControls.H" #include "readDyMControls.H"
// Store divU from the previous mesh so that it can be mapped // Store divU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the // and used in correctPhi to ensure the corrected phi has the
// same divergence // same divergence
volScalarField divU("divU", fvc::div(fvc::absolute(phi, U))); volScalarField divU("divU", fvc::div(fvc::absolute(phi, U)));
{ {
#include "CourantNo.H" #include "CourantNo.H"
#include "alphaCourantNo.H" #include "alphaCourantNo.H"
@ -157,7 +157,7 @@ int main(int argc, char *argv[])
} }
} }
} }
mixture->correct(); mixture->correct();
#include "alphaControls.H" #include "alphaControls.H"
@ -165,7 +165,7 @@ int main(int argc, char *argv[])
#include "UEqn.H" #include "UEqn.H"
#include "TEqn.H" #include "TEqn.H"
// --- Pressure corrector loop // --- Pressure corrector loop
while (pimple.correct()) while (pimple.correct())
{ {

View File

@ -58,14 +58,14 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::constant
temperaturePhaseChangeTwoPhaseMixture(mixture, mesh), temperaturePhaseChangeTwoPhaseMixture(mixture, mesh),
coeffC_ coeffC_
( (
"coeffC", "coeffC",
dimless/dimTime/dimTemperature, dimless/dimTime/dimTemperature,
optionalSubDict(type() + "Coeffs") optionalSubDict(type() + "Coeffs")
), ),
coeffE_ coeffE_
( (
"coeffE", "coeffE",
dimless/dimTime/dimTemperature, dimless/dimTime/dimTemperature,
optionalSubDict(type() + "Coeffs") optionalSubDict(type() + "Coeffs")
) )
{} {}
@ -121,7 +121,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::mDot() const
const dimensionedScalar& TSat = thermo.TSat(); const dimensionedScalar& TSat = thermo.TSat();
const dimensionedScalar T0(dimTemperature, Zero); const dimensionedScalar T0(dimTemperature, Zero);
if (mesh_.time().outputTime()) if (mesh_.time().outputTime())
{ {
volScalarField mDot volScalarField mDot
@ -194,9 +194,9 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
); );
const dimensionedScalar& TSat = thermo.TSat(); const dimensionedScalar& TSat = thermo.TSat();
dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1(); dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
volScalarField limitedAlpha1 volScalarField limitedAlpha1
( (
min(max(mixture_.alpha1(), scalar(0)), scalar(1)) min(max(mixture_.alpha1(), scalar(0)), scalar(1))
@ -206,7 +206,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
( (
min(max(mixture_.alpha2(), scalar(0)), scalar(1)) min(max(mixture_.alpha2(), scalar(0)), scalar(1))
); );
const volScalarField Vcoeff const volScalarField Vcoeff
( (
coeffE_*mixture_.rho1()*limitedAlpha1*L coeffE_*mixture_.rho1()*limitedAlpha1*L
@ -216,7 +216,7 @@ Foam::temperaturePhaseChangeTwoPhaseMixtures::constant::TSource() const
coeffC_*mixture_.rho2()*limitedAlpha2*L coeffC_*mixture_.rho2()*limitedAlpha2*L
); );
TSource = TSource =
fvm::Sp(Vcoeff, T) - Vcoeff*TSat fvm::Sp(Vcoeff, T) - Vcoeff*TSat
- fvm::Sp(Ccoeff, T) + Ccoeff*TSat; - fvm::Sp(Ccoeff, T) + Ccoeff*TSat;

View File

@ -98,7 +98,7 @@ public:
// coefficient to multiply (Tsat - T) for the condensation rate // coefficient to multiply (Tsat - T) for the condensation rate
// and a coefficient to multiply (T - Tsat) for the vaporisation rate // and a coefficient to multiply (T - Tsat) for the vaporisation rate
virtual Pair<tmp<volScalarField>> mDotDeltaT() const; virtual Pair<tmp<volScalarField>> mDotDeltaT() const;
//- Source for T equarion //- Source for T equarion
virtual tmp<fvScalarMatrix> TSource() const; virtual tmp<fvScalarMatrix> TSource() const;

View File

@ -64,10 +64,10 @@ interfaceHeatResistance
temperaturePhaseChangeTwoPhaseMixture(mixture, mesh), temperaturePhaseChangeTwoPhaseMixture(mixture, mesh),
R_ R_
( (
"R", "R",
dimPower/dimArea/dimTemperature, optionalSubDict(type() + "Coeffs") dimPower/dimArea/dimTemperature, optionalSubDict(type() + "Coeffs")
), ),
interfaceArea_ interfaceArea_
( (
IOobject IOobject
@ -81,7 +81,7 @@ interfaceHeatResistance
mesh_, mesh_,
dimensionedScalar(dimless/dimLength, Zero) dimensionedScalar(dimless/dimLength, Zero)
), ),
mDotc_ mDotc_
( (
IOobject IOobject
@ -95,7 +95,7 @@ interfaceHeatResistance
mesh_, mesh_,
dimensionedScalar(dimDensity/dimTime, Zero) dimensionedScalar(dimDensity/dimTime, Zero)
), ),
mDote_ mDote_
( (
IOobject IOobject
@ -109,7 +109,7 @@ interfaceHeatResistance
mesh_, mesh_,
dimensionedScalar(dimDensity/dimTime, Zero) dimensionedScalar(dimDensity/dimTime, Zero)
), ),
spread_ spread_
( (
optionalSubDict(type() + "Coeffs").get<scalar>("spread") optionalSubDict(type() + "Coeffs").get<scalar>("spread")
@ -148,7 +148,7 @@ mDotAlphal() const
( (
min(max(mixture_.alpha2(), scalar(0)), scalar(1)) min(max(mixture_.alpha2(), scalar(0)), scalar(1))
); );
return Pair<tmp<volScalarField>> return Pair<tmp<volScalarField>>
( (
(mDotc_/(limitedAlpha2 + SMALL)), (mDotc_/(limitedAlpha2 + SMALL)),
@ -161,10 +161,10 @@ Foam::Pair<Foam::tmp<Foam::volScalarField>>
Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance:: Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
mDot() const mDot() const
{ {
return Pair<tmp<volScalarField>> return Pair<tmp<volScalarField>>
( (
tmp<volScalarField>(mDotc_), tmp<volScalarField>(mDotc_),
tmp<volScalarField>(mDote_) tmp<volScalarField>(mDote_)
); );
} }
@ -179,20 +179,20 @@ mDotDeltaT() const
( (
mesh_.lookupObject<basicThermo>(basicThermo::dictName) mesh_.lookupObject<basicThermo>(basicThermo::dictName)
); );
const volScalarField& T = mesh_.lookupObject<volScalarField>("T"); const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
const dimensionedScalar& TSat = thermo.TSat(); const dimensionedScalar& TSat = thermo.TSat();
Pair<tmp<volScalarField>> mDotce(mDot()); Pair<tmp<volScalarField>> mDotce(mDot());
return Pair<tmp<volScalarField>> return Pair<tmp<volScalarField>>
( (
mDotc_*pos(TSat - T.oldTime())/(TSat - T.oldTime()), mDotc_*pos(TSat - T.oldTime())/(TSat - T.oldTime()),
-mDote_*pos(T.oldTime() - TSat)/(T.oldTime() - TSat) -mDote_*pos(T.oldTime() - TSat)/(T.oldTime() - TSat)
); );
} }
@ -221,7 +221,7 @@ TSource() const
); );
const dimensionedScalar& TSat = thermo.TSat(); const dimensionedScalar& TSat = thermo.TSat();
// interface heat resistance // interface heat resistance
volScalarField IHRcoeff = interfaceArea_*R_; volScalarField IHRcoeff = interfaceArea_*R_;
@ -234,10 +234,10 @@ TSource() const
void Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance:: void Foam::temperaturePhaseChangeTwoPhaseMixtures::interfaceHeatResistance::
correct() correct()
{ {
// Update Interface // Update Interface
updateInterface(); updateInterface();
// Update mDotc_ and mDote_ // Update mDotc_ and mDote_
const volScalarField& T = mesh_.lookupObject<volScalarField>("T"); const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
@ -251,16 +251,16 @@ correct()
const dimensionedScalar T0(dimTemperature, Zero); const dimensionedScalar T0(dimTemperature, Zero);
dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1(); dimensionedScalar L = mixture_.Hf2() - mixture_.Hf1();
// interface heat resistance // interface heat resistance
mDotc_ = interfaceArea_*R_*max(TSat - T, T0)/L; mDotc_ = interfaceArea_*R_*max(TSat - T, T0)/L;
mDote_ = interfaceArea_*R_*(T - TSat)/L; mDote_ = interfaceArea_*R_*(T - TSat)/L;
forAll(mDotc_, celli) forAll(mDotc_, celli)
{ {
scalar rhobyDt = mixture_.rho1().value()/mesh_.time().deltaTValue(); scalar rhobyDt = mixture_.rho1().value()/mesh_.time().deltaTValue();
scalar maxEvap = mixture_.alpha1()[celli]*rhobyDt; // positive scalar maxEvap = mixture_.alpha1()[celli]*rhobyDt; // positive
scalar maxCond = -mixture_.alpha2()[celli]*rhobyDt; // negative scalar maxCond = -mixture_.alpha2()[celli]*rhobyDt; // negative
mDote_[celli] = min(max(mDote_[celli], maxCond), maxEvap); mDote_[celli] = min(max(mDote_[celli], maxCond), maxEvap);
mDotc_[celli] = min(max(mDotc_[celli], maxCond), maxEvap); mDotc_[celli] = min(max(mDotc_[celli], maxCond), maxEvap);
} }
@ -278,7 +278,7 @@ updateInterface()
); );
const dimensionedScalar& TSat = thermo.TSat(); const dimensionedScalar& TSat = thermo.TSat();
// interface heat resistance // interface heat resistance
// Interpolating alpha1 cell centre values to mesh points (vertices) // Interpolating alpha1 cell centre values to mesh points (vertices)
scalarField ap scalarField ap
@ -294,7 +294,7 @@ updateInterface()
interfaceArea_[celli] = 0; interfaceArea_[celli] = 0;
if (status == 0) // cell is cut if (status == 0) // cell is cut
{ {
interfaceArea_[celli] = interfaceArea_[celli] =
mag(cutCell.isoFaceArea())/mesh_.V()[celli]; mag(cutCell.isoFaceArea())/mesh_.V()[celli];
} }
} }
@ -309,14 +309,14 @@ updateInterface()
forAll(pp.faceCells(),i) forAll(pp.faceCells(),i)
{ {
const label pCelli = pp.faceCells()[i]; const label pCelli = pp.faceCells()[i];
if if
( (
(TSat.value() - T[pCelli]) > 0 (TSat.value() - T[pCelli]) > 0
&& mixture_.alpha1()[pCelli] < 0.9 && mixture_.alpha1()[pCelli] < 0.9
) )
{ {
interfaceArea_[pCelli] = interfaceArea_[pCelli] =
mag(pp.faceAreas()[i])/mesh_.V()[pCelli]; mag(pp.faceAreas()[i])/mesh_.V()[pCelli];
} }
} }
@ -331,23 +331,23 @@ vDot() const
dimensionedScalar D dimensionedScalar D
( (
"D", "D",
dimArea, dimArea,
spread_/sqr(gAverage(mesh_.nonOrthDeltaCoeffs())) spread_/sqr(gAverage(mesh_.nonOrthDeltaCoeffs()))
); );
const volScalarField& alpha1 = mixture_.alpha1(); const volScalarField& alpha1 = mixture_.alpha1();
const volScalarField& alpha2 = mixture_.alpha2(); const volScalarField& alpha2 = mixture_.alpha2();
const dimensionedScalar MDotMin("MdotMin", mDotc_.dimensions(), 1e-3); const dimensionedScalar MDotMin("MdotMin", mDotc_.dimensions(), 1e-3);
Pair<tmp<volScalarField>> mDotSpread Pair<tmp<volScalarField>> mDotSpread
( (
tmp<volScalarField>(mDotc_*0.0), tmp<volScalarField>(mDotc_*0.0),
tmp<volScalarField>(mDote_*0.0) tmp<volScalarField>(mDote_*0.0)
); );
if (max(mDotc_) > MDotMin) if (max(mDotc_) > MDotMin)
{ {
fvc::spreadSource fvc::spreadSource
@ -373,18 +373,18 @@ vDot() const
1e-3 1e-3
); );
} }
dimensionedScalar pCoeff(1.0/mixture_.rho1() - 1.0/mixture_.rho2()); dimensionedScalar pCoeff(1.0/mixture_.rho1() - 1.0/mixture_.rho2());
if (mesh_.time().outputTime()) if (mesh_.time().outputTime())
{ {
volScalarField mDotS("mDotSpread", mDotSpread[1].ref()); volScalarField mDotS("mDotSpread", mDotSpread[1].ref());
mDotS.write(); mDotS.write();
} }
return Pair<tmp<volScalarField>> return Pair<tmp<volScalarField>>
( (
pCoeff*mDotSpread[0], pCoeff*mDotSpread[0],
-pCoeff*mDotSpread[1] -pCoeff*mDotSpread[1]
); );
} }

View File

@ -28,9 +28,9 @@ Class
Foam::phaseChangeTwoPhaseMixtures::interfaceHeatResistance Foam::phaseChangeTwoPhaseMixtures::interfaceHeatResistance
Description Description
Interface Heat Resistance type of condensation/saturation model using Interface Heat Resistance type of condensation/saturation model using
spread source distribution following: spread source distribution following:
References: References:
\verbatim \verbatim
Hardt, S., Wondra, F. (2008). Hardt, S., Wondra, F. (2008).
@ -73,22 +73,22 @@ class interfaceHeatResistance
//- Interface area //- Interface area
volScalarField interfaceArea_; volScalarField interfaceArea_;
//- Mass condensation source //- Mass condensation source
volScalarField mDotc_; volScalarField mDotc_;
//- Mass evaporation source //- Mass evaporation source
volScalarField mDote_; volScalarField mDote_;
//- Spread for mass source //- Spread for mass source
scalar spread_; scalar spread_;
// Private member functions // Private member functions
//- Update interface area //- Update interface area
void updateInterface(); void updateInterface();
public: public:
@ -124,17 +124,17 @@ public:
// coefficient to multiply (Tsat - T) for the condensation rate // coefficient to multiply (Tsat - T) for the condensation rate
// and a coefficient to multiply (T - Tsat) for the vaporisation rate // and a coefficient to multiply (T - Tsat) for the vaporisation rate
virtual Pair<tmp<volScalarField>> mDotDeltaT() const; virtual Pair<tmp<volScalarField>> mDotDeltaT() const;
//- Source for T equarion //- Source for T equarion
virtual tmp<fvScalarMatrix> TSource() const; virtual tmp<fvScalarMatrix> TSource() const;
//- Volumetric source for alpha (used by alphaEq) //- Volumetric source for alpha (used by alphaEq)
virtual Pair<tmp<volScalarField>> vDotAlphal() const; virtual Pair<tmp<volScalarField>> vDotAlphal() const;
//- Return the volumetric condensation and vaporisation rates as //- Return the volumetric condensation and vaporisation rates as
// coefficients (used by p_rghEq) // coefficients (used by p_rghEq)
virtual Pair<tmp<volScalarField>> vDot() const; virtual Pair<tmp<volScalarField>> vDot() const;
//- Correct the interfaceHeatResistance phaseChange model //- Correct the interfaceHeatResistance phaseChange model
virtual void correct(); virtual void correct();

View File

@ -68,7 +68,7 @@ temperaturePhaseChangeTwoPhaseMixture
Foam::Pair<Foam::tmp<Foam::volScalarField>> Foam::Pair<Foam::tmp<Foam::volScalarField>>
Foam::temperaturePhaseChangeTwoPhaseMixture::vDotAlphal() const Foam::temperaturePhaseChangeTwoPhaseMixture::vDotAlphal() const
{ {
volScalarField alphalCoeff volScalarField alphalCoeff
( (
1.0/mixture_.rho1() - mixture_.alpha1() 1.0/mixture_.rho1() - mixture_.alpha1()

View File

@ -141,7 +141,7 @@ public:
// coefficient to multiply (Tsat - T) for the condensation rate // coefficient to multiply (Tsat - T) for the condensation rate
// and a coefficient to multiply (T - Tsat) for the vaporisation rate // and a coefficient to multiply (T - Tsat) for the vaporisation rate
virtual Pair<tmp<volScalarField>> mDotDeltaT() const = 0; virtual Pair<tmp<volScalarField>> mDotDeltaT() const = 0;
//- Source for T equarion //- Source for T equarion
virtual tmp<fvScalarMatrix> TSource() const = 0; virtual tmp<fvScalarMatrix> TSource() const = 0;

View File

@ -114,7 +114,7 @@
} }
rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2; rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
// Cache alphaPhi // Cache alphaPhi
alphaPhi10 = talphaPhi(); alphaPhi10 = talphaPhi();

View File

@ -332,7 +332,7 @@ void Foam::fvc::spreadSource
) )
{ {
const fvMesh& mesh = alpha1.mesh(); const fvMesh& mesh = alpha1.mesh();
volScalarField mDotSmear volScalarField mDotSmear
( (
IOobject IOobject
@ -351,43 +351,43 @@ void Foam::fvc::spreadSource
//- Smearing of source term field //- Smearing of source term field
fvScalarMatrix mSourceEqn fvScalarMatrix mSourceEqn
( (
fvm::Sp(scalar(1), mDotSmear) fvm::Sp(scalar(1), mDotSmear)
- fvm::laplacian(D, mDotSmear) - fvm::laplacian(D, mDotSmear)
== ==
mDotIn mDotIn
); );
mSourceEqn.solve(); mSourceEqn.solve();
// Cut cells with cutoff < alpha1 < 1-cutoff and rescale remaining // Cut cells with cutoff < alpha1 < 1-cutoff and rescale remaining
// source term field // source term field
dimensionedScalar intvDotLiquid("intvDotLiquid", dimMass/dimTime, 0.0); dimensionedScalar intvDotLiquid("intvDotLiquid", dimMass/dimTime, 0.0);
dimensionedScalar intvDotVapor ("intvDotVapor", dimMass/dimTime, 0.0); dimensionedScalar intvDotVapor ("intvDotVapor", dimMass/dimTime, 0.0);
const scalarField& Vol = mesh.V(); const scalarField& Vol = mesh.V();
forAll(mesh.C(), celli) forAll(mesh.C(), celli)
{ {
if (alpha1[celli] < cutoff) if (alpha1[celli] < cutoff)
{ {
intvDotVapor.value() += intvDotVapor.value() +=
alpha2[celli]*mDotSmear[celli]*Vol[celli]; alpha2[celli]*mDotSmear[celli]*Vol[celli];
} }
else if (alpha1[celli] > 1.0 - cutoff) else if (alpha1[celli] > 1.0 - cutoff)
{ {
intvDotLiquid.value() += intvDotLiquid.value() +=
alpha1[celli]*mDotSmear[celli]*Vol[celli]; alpha1[celli]*mDotSmear[celli]*Vol[celli];
} }
} }
reduce(intvDotVapor.value(), sumOp<scalar>()); reduce(intvDotVapor.value(), sumOp<scalar>());
reduce(intvDotLiquid.value(), sumOp<scalar>()); reduce(intvDotLiquid.value(), sumOp<scalar>());
//- Calculate Nl and Nv //- Calculate Nl and Nv
dimensionedScalar Nl ("Nl", dimless, Zero); dimensionedScalar Nl ("Nl", dimless, Zero);
dimensionedScalar Nv ("Nv", dimless, Zero); dimensionedScalar Nv ("Nv", dimless, Zero);
const dimensionedScalar intmSource0(fvc::domainIntegrate(mDotIn)); const dimensionedScalar intmSource0(fvc::domainIntegrate(mDotIn));
if (intvDotVapor.value() > VSMALL) if (intvDotVapor.value() > VSMALL)

View File

@ -44,11 +44,11 @@ Description
gradient of alpha is large (where the difference between the values gradient of alpha is large (where the difference between the values
in neighbouring cells is larger than alphaDiff) away from that in neighbouring cells is larger than alphaDiff) away from that
starting point of the sweep. starting point of the sweep.
spreadSource: spread a source field (mDotIn) for two phase multiphase using spreadSource: spread a source field (mDotIn) for two phase multiphase using
a laplacian operator and diffussivity D. a laplacian operator and diffussivity D.
The spread source (mDotOut) is distributed from alpha1 < cutoff The spread source (mDotOut) is distributed from alpha1 < cutoff
to alpha1 > 1 - cutoff, and it is zero across the interface to alpha1 > 1 - cutoff, and it is zero across the interface
SourceFiles SourceFiles
@ -91,7 +91,7 @@ namespace fvc
const label nLayers, const label nLayers,
const scalar alphaDiff = 0.2 const scalar alphaDiff = 0.2
); );
void spreadSource void spreadSource
( (
volScalarField& mDotOut, volScalarField& mDotOut,

View File

@ -55,7 +55,7 @@ void Foam::functionObjects::interfaceHeight::writePositions()
const uniformDimensionedVectorField& g = const uniformDimensionedVectorField& g =
mesh_.time().lookupObject<uniformDimensionedVectorField>("g"); mesh_.time().lookupObject<uniformDimensionedVectorField>("g");
vector gHat = vector::zero; vector gHat = vector::zero;
if (mag(direction_) > 0.0) if (mag(direction_) > 0.0)
{ {
gHat = direction_/mag(direction_); gHat = direction_/mag(direction_);
@ -237,7 +237,7 @@ Foam::functionObjects::interfaceHeight::interfaceHeight
{ {
read(dict); read(dict);
resetNames({"height", "position"}); resetNames({"height", "position"});
writeFileHeader(fileID::heightFile); writeFileHeader(fileID::heightFile);
writeFileHeader(fileID::positionFile); writeFileHeader(fileID::positionFile);
} }
@ -253,7 +253,7 @@ Foam::functionObjects::interfaceHeight::~interfaceHeight()
bool Foam::functionObjects::interfaceHeight::read(const dictionary& dict) bool Foam::functionObjects::interfaceHeight::read(const dictionary& dict)
{ {
dict.readIfPresent("alpha", alphaName_); dict.readIfPresent("alpha", alphaName_);
dict.readIfPresent("liquid", liquid_); dict.readIfPresent("liquid", liquid_);
dict.lookup("locations") >> locations_; dict.lookup("locations") >> locations_;

View File

@ -96,9 +96,8 @@ class interfaceHeight
//- Interpolation scheme //- Interpolation scheme
word interpolationScheme_; word interpolationScheme_;
//- Direction of interface motion
//- Direction of interface motion
vector direction_; vector direction_;

View File

@ -24,7 +24,7 @@ boundaryField
{ {
type zeroGradient; type zeroGradient;
} }
outlet outlet
{ {
type inletOutlet; type inletOutlet;

View File

@ -24,7 +24,7 @@ boundaryField
{ {
type zeroGradient; type zeroGradient;
} }
outlet outlet
{ {
type inletOutlet; type inletOutlet;

View File

@ -1,11 +1,12 @@
#!/bin/sh #!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory cd ${0%/*} || exit 1 # Run from this directory
. $WM_PROJECT_DIR/bin/tools/RunFunctions # Tutorial run functions
# Source tutorial run functions #------------------------------------------------------------------------------
. $WM_PROJECT_DIR/bin/tools/RunFunctions
runApplication blockMesh runApplication blockMesh
restore0Dir restore0Dir
# copy 0 folder to 1.36
cp -r 0 1.36 cp -r 0 1.36
cp system/setAlphaFieldDict.liquid system/setAlphaFieldDict cp system/setAlphaFieldDict.liquid system/setAlphaFieldDict
runApplication setAlphaField runApplication setAlphaField

View File

@ -1,8 +1,7 @@
#!/bin/sh #!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory cd ${0%/*} || exit 1 # Run from this directory
. $WM_PROJECT_DIR/bin/tools/RunFunctions # Tutorial run functions
# Source tutorial run functions #------------------------------------------------------------------------------
. $WM_PROJECT_DIR/bin/tools/RunFunctions
runApplication blockMesh runApplication blockMesh
restore0Dir restore0Dir

View File

@ -31,7 +31,7 @@ divSchemes
div(rhoPhi,U) Gauss vanLeerV; div(rhoPhi,U) Gauss vanLeerV;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(rhoCpPhi,T) Gauss vanLeer; div(rhoCpPhi,T) Gauss vanLeer;
div((interpolate(cp)*rhoPhi),T) Gauss linear; div((interpolate(cp)*rhoPhi),T) Gauss linear;
} }