functionObjects: Added phaseScalarTransport function

This is like the scalarTrasport function except that the transported
scalar is confined to a single phase of a multiphase simulation. In
addition to the usual specification for the scalarTransport function
(i.e., a field, schemes and solution parameters), the user needs to
specify the phase-flux or a pressure field which can be used to generate
it.

Example usage for interFoam:

    phaseScalarTransport1
    {
        type            phaseScalarTransport;
        libs            ("libsolverFunctionObjects.so");

        field           s.water;
        p               p_rgh;
    }

Example usage for reactingTwoPhaseEulerFoam:

    phaseScalarTransport1
    {
        type            phaseScalarTransport;
        libs            ("libsolverFunctionObjects.so");

        field           s.water;
        alphaPhi        alphaRhoPhi.water;
        rho             thermo:rho.water;
    }

The function will write out both the per-unit-phase field that is solved
for (s.water in the above examples) and also the mixture-total field
(alphaS.water), which is often more convenient for post-processing.
This commit is contained in:
Will Bainbridge
2019-02-14 15:46:00 +00:00
parent e9f3811218
commit 7b1840c7d3
6 changed files with 773 additions and 2 deletions

View File

@ -1,3 +1,4 @@
scalarTransport/scalarTransport.C
phaseScalarTransport/phaseScalarTransport.C
LIB = $(FOAM_LIBBIN)/libsolverFunctionObjects

View File

@ -0,0 +1,500 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "addToRunTimeSelectionTable.H"
#include "fixedValueFvPatchField.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
#include "fvmDdt.H"
#include "fvmDiv.H"
#include "fvmLaplacian.H"
#include "nonOrthogonalSolutionControl.H"
#include "phaseScalarTransport.H"
#include "surfaceFields.H"
#include "turbulenceModel.H"
#include "wallFvPatch.H"
#include "zeroGradientFvPatchField.H"
#define PhiDimensionErrorInFunction(phi) \
FatalErrorInFunction \
<< "Incompatible dimensions for " << phi.name() << ": " \
<< phi.dimensions() << nl \
<< "Dimensions should be " << dimMass/dimTime << " or " \
<< dimVolume/dimTime << exit(FatalError)
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
defineTypeNameAndDebug(phaseScalarTransport, 0);
addToRunTimeSelectionTable
(
functionObject,
phaseScalarTransport,
dictionary
);
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::volScalarField& Foam::functionObjects::phaseScalarTransport::Phi()
{
if (!PhiPtr_.valid())
{
const surfaceScalarField& phi =
mesh_.lookupObject<surfaceScalarField>(phiName_);
const volScalarField& p =
mesh_.lookupObject<volScalarField>(pName_);
wordList PhiPatchFieldTypes(mesh_.boundaryMesh().size());
forAll(p.boundaryField(), patchi)
{
PhiPatchFieldTypes[patchi] =
p.boundaryField()[patchi].fixesValue()
? fixedValueFvPatchField<scalar>::typeName
: zeroGradientFvPatchField<scalar>::typeName;
}
PhiPtr_.set
(
new volScalarField
(
IOobject
(
"Phi" + s_.name(),
mesh_.time().timeName(),
mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(phi.dimensions()/dimLength, Zero),
PhiPatchFieldTypes
)
);
mesh_.setFluxRequired(PhiPtr_->name());
}
return PhiPtr_();
}
Foam::tmp<Foam::surfaceScalarField>
Foam::functionObjects::phaseScalarTransport::alphaPhi()
{
// If alphaPhi exists then return it
if (mesh_.foundObject<surfaceScalarField>(alphaPhiName_))
{
return mesh_.lookupObject<surfaceScalarField>(alphaPhiName_);
}
// Otherwise generate it ...
Info<< type() << ": " << surfaceScalarField::typeName << " "
<< alphaPhiName_ << " was not found, so generating it" << endl;
const volScalarField& alpha =
mesh_.lookupObject<volScalarField>(alphaName_);
const surfaceScalarField& phi =
mesh_.lookupObject<surfaceScalarField>(phiName_);
// Make a crude guess of the phase flux using default interpolation
tmp<surfaceScalarField> tAlphaPhi
(
new surfaceScalarField
(
alphaPhiName_,
phi*fvc::interpolate(alpha)
)
);
surfaceScalarField& alphaPhi = tAlphaPhi.ref();
// Get the potential field
volScalarField& Phi(this->Phi());
// Construct the scheme names
const word laplacianScheme = "laplacian(" + pName_ + ")";
// Debug writing. Write the material derivative of alpha, before and after
// the solution of the potential and the correction of alphaPhi. Before
// correction the field should be non-zero, and after it should be
// comparable to the solution tolerance.
auto writeDDt = [&](const label i)
{
const volScalarField DDtAlpha
(
"DDt("
+ IOobject::groupName
(
IOobject::member(alpha.name()) + Foam::name(i),
IOobject::group(alpha.name())
)
+ ")",
fvc::ddt(alpha) + fvc::div(alphaPhi)
);
Info<< type() << ": Writing " << DDtAlpha.name() << endl;
DDtAlpha.write();
};
if (debug && mesh_.time().writeTime())
{
writeDDt(0);
}
// Construct a non-orthogonal solution control
nonOrthogonalSolutionControl control
(
const_cast<fvMesh&>(mesh_),
mesh_.lookupObjectRef<nonOrthogonalSolutionControl>
(
solutionControl::typeName
)
.algorithmName()
);
// Solve for the potential and correct alphaPhi with the resulting flux
if (phi.dimensions() == dimVolume/dimTime)
{
while (control.correctNonOrthogonal())
{
fvScalarMatrix PhiEqn
(
fvm::laplacian(Phi, laplacianScheme)
+ fvc::ddt(alpha)
+ fvc::div(alphaPhi)
);
PhiEqn.solve(pName_);
if (control.finalNonOrthogonalIter())
{
alphaPhi += PhiEqn.flux();
}
}
}
else if (phi.dimensions() == dimMass/dimTime)
{
const volScalarField& rho =
mesh_.lookupObject<volScalarField>(rhoName_);
while (control.correctNonOrthogonal())
{
fvScalarMatrix PhiEqn
(
fvm::laplacian(Phi, laplacianScheme)
+ fvc::ddt(rho, alpha)
+ fvc::div(alphaPhi)
);
PhiEqn.solve(pName_);
if (control.finalNonOrthogonalIter())
{
alphaPhi += PhiEqn.flux();
}
}
}
else
{
PhiDimensionErrorInFunction(phi);
}
// Debug writing
if (debug && mesh_.time().writeTime())
{
writeDDt(1);
}
return tAlphaPhi;
}
Foam::tmp<Foam::volScalarField>
Foam::functionObjects::phaseScalarTransport::D
(
const surfaceScalarField& alphaPhi
) const
{
if (constantD_)
{
return volScalarField::New
(
"D" + s_.name(),
mesh_,
dimensionedScalar(alphaPhi.dimensions()/dimLength, D_)
);
}
const word& nameNoPhase = turbulenceModel::propertiesName;
const word namePhase = IOobject::groupName(nameNoPhase, phaseName_);
const word& name =
mesh_.foundObject<turbulenceModel>(namePhase)
? namePhase
: mesh_.foundObject<turbulenceModel>(nameNoPhase)
? nameNoPhase
: word::null;
if (name == word::null)
{
return volScalarField::New
(
"D" + s_.name(),
mesh_,
dimensionedScalar(alphaPhi.dimensions()/dimLength, 0)
);
}
const turbulenceModel& turbulence =
mesh_.lookupObject<turbulenceModel>(name);
if (alphaPhi.dimensions() == dimVolume/dimTime)
{
return alphaD_*turbulence.nu() + alphaDt_*turbulence.nut();
}
else if (alphaPhi.dimensions() == dimMass/dimTime)
{
return alphaD_*turbulence.mu() + alphaDt_*turbulence.mut();
}
else
{
PhiDimensionErrorInFunction(alphaPhi);
return tmp<volScalarField>(nullptr);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::functionObjects::phaseScalarTransport::phaseScalarTransport
(
const word& name,
const Time& runTime,
const dictionary& dict
)
:
fvMeshFunctionObject(name, runTime, dict),
fieldName_(dict.lookup("field")),
phaseName_(IOobject::group(fieldName_)),
nCorr_(0),
residualAlpha_(rootSmall),
fvOptions_(mesh_),
s_
(
IOobject
(
fieldName_,
mesh_.time().timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),
alphaS_
(
IOobject
(
"alpha"
+ word(toupper(fieldName_[0]))
+ fieldName_(1, fieldName_.size() - 1),
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh_,
dimensionedScalar(s_.dimensions(), Zero)
),
PhiPtr_(nullptr)
{
if (phaseName_ == word::null)
{
FatalErrorInFunction
<< "Field \"" << fieldName_ << "\" does not have a phase extension "
<< "in its name. If it is associated with \"phaseA\" then it "
<< "should be named \"" << fieldName_ << ".phaseA\"."
<< exit(FatalError);
}
read(dict);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::functionObjects::phaseScalarTransport::~phaseScalarTransport()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::functionObjects::phaseScalarTransport::read(const dictionary& dict)
{
fvMeshFunctionObject::read(dict);
alphaName_ =
dict.lookupOrDefault<word>
(
"alpha",
IOobject::groupName("alpha", phaseName_)
);
alphaPhiName_ =
dict.lookupOrDefault<word>
(
"alphaPhi",
IOobject::groupName("alphaPhi", phaseName_)
);
phiName_ = dict.lookupOrDefault<word>("phi", "phi");
rhoName_ =
dict.lookupOrDefault<word>
(
"rho",
IOobject::groupName("rho", phaseName_)
);
pName_ = dict.lookupOrDefault<word>("p", "p");
schemesField_ = dict.lookupOrDefault<word>("schemesField", fieldName_);
constantD_ = dict.readIfPresent("D", D_);
alphaD_ = dict.lookupOrDefault<scalar>("alphaD", 1);
alphaDt_ = dict.lookupOrDefault<scalar>("alphaDt", 1);
dict.readIfPresent("nCorr", nCorr_);
dict.readIfPresent("residualAlpha", residualAlpha_);
if (dict.found("fvOptions"))
{
fvOptions_.reset(dict.subDict("fvOptions"));
}
return true;
}
bool Foam::functionObjects::phaseScalarTransport::execute()
{
Info<< type() << ": Executing" << endl;
const volScalarField& alpha =
mesh_.lookupObject<volScalarField>(alphaName_);
// Get the phase flux
tmp<surfaceScalarField> tAlphaPhi(this->alphaPhi());
const surfaceScalarField& alphaPhi = tAlphaPhi();
// Get the diffusivity
const volScalarField D(this->D(alphaPhi));
// Construct the scheme names
const word divScheme =
"div(" + alphaPhi.name() + "," + schemesField_ + ")";
const word laplacianScheme =
"laplacian(" + D.name() + "," + schemesField_ + ")";
// Get the relaxation coefficient
const scalar relaxCoeff =
mesh_.relaxEquation(schemesField_)
? mesh_.equationRelaxationFactor(schemesField_)
: 0;
// Solve
if (alphaPhi.dimensions() == dimVolume/dimTime)
{
for (label i = 0; i <= nCorr_; ++ i)
{
fvScalarMatrix fieldEqn
(
fvm::ddt(alpha, s_)
+ fvm::div(alphaPhi, s_, divScheme)
- fvm::laplacian
(
fvc::interpolate(alpha)*fvc::interpolate(D),
s_,
laplacianScheme
)
==
fvOptions_(alpha, s_)
- fvm::ddt(residualAlpha_, s_)
+ fvc::ddt(residualAlpha_, s_)
);
fieldEqn.relax(relaxCoeff);
fvOptions_.constrain(fieldEqn);
fieldEqn.solve(schemesField_);
}
}
else if (alphaPhi.dimensions() == dimMass/dimTime)
{
const volScalarField& rho =
mesh_.lookupObject<volScalarField>(rhoName_);
for (label i = 0; i <= nCorr_; ++ i)
{
fvScalarMatrix fieldEqn
(
fvm::ddt(alpha, rho, s_)
+ fvm::div(alphaPhi, s_, divScheme)
- fvm::laplacian
(
fvc::interpolate(alpha)*fvc::interpolate(D),
s_,
laplacianScheme
)
==
fvOptions_(rho, s_)
- fvm::ddt(residualAlpha_*rho, s_)
+ fvc::ddt(residualAlpha_*rho, s_)
);
fieldEqn.relax(relaxCoeff);
fvOptions_.constrain(fieldEqn);
fieldEqn.solve(schemesField_);
}
}
else
{
PhiDimensionErrorInFunction(alphaPhi);
}
// Update
alphaS_ = alpha*s_;
Info<< endl;
return true;
}
bool Foam::functionObjects::phaseScalarTransport::write()
{
return true;
}
// ************************************************************************* //

View File

@ -0,0 +1,233 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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/>.
Class
Foam::functionObjects::phaseScalarTransport
Description
Evolves a passive scalar transport equation within one phase of a
multiphase simulation. The scalar is considered to be a phase-intensive
property; i.e., its value represents an amount per-unit of the phase. In
addition to the scalar, the function also writes out the product of the
volume fraction and the scalar, as this provides a phase-extensive field
which is often more convenient to post-process.
Most entries are the same as for the \c scalarTransport function. Refer to
its documentation for details. Entries specific to this function are
detailed below. Note that the phase-name will be determined by stripping
the extension from the supplied field name.
If \c alphaPhi is specified and found in the database then it will be used
to transport the field. This is the likely mode of operation if this
function is used with Euler-Euler solvers, in which \c alphaPhi -like
fluxes are available. If \c alphaPhi is not found, then a pressure-like
equation will be solved in order to construct it so that it exactly matches
the time-derivative of the phase volume or mass. This is likely to be
necessary in volume-of-fluid solvers where \c alphaPhi is not part of the
solution procedure. The pressure field name will be required in this case.
Usage
\table
Property | Description | Req'd? | Default
alpha | Name of the volume-fraction field | no \\
| alpha.<phase-name>
alphaPhi | Name of the phase-flux field | no \\
| alphaPhi.<phase-name>
p | Name of the pressure field | no | p
residualAlpha | Small volume fraction used to stabilise the solution \\
| no | rootSmall
\endtable
Example specification for interFoam:
\verbatim
phaseScalarTransport1
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
field s.water;
p p_rgh;
}
\endverbatim
Example specification for reactingTwoPhaseEulerFoam:
\verbatim
phaseScalarTransport1
{
type phaseScalarTransport;
libs ("libsolverFunctionObjects.so");
field s.water;
alphaPhi alphaRhoPhi.water;
rho thermo:rho.water;
}
\endverbatim
See also
Foam::functionObjects::fvMeshFunctionObject
Foam::functionObjects::scalarTransport
SourceFiles
phaseScalarTransport.C
\*---------------------------------------------------------------------------*/
#ifndef functionObjects_phaseScalarTransport_H
#define functionObjects_phaseScalarTransport_H
#include "fvMeshFunctionObject.H"
#include "volFields.H"
#include "fvOptionList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
/*---------------------------------------------------------------------------*\
Class phaseScalarTransport Declaration
\*---------------------------------------------------------------------------*/
class phaseScalarTransport
:
public fvMeshFunctionObject
{
// Private data
//- Name of field to process
const word fieldName_;
//- Name of the phase in which to solve
const word phaseName_;
//- Name of phase volume-fraction field (optional)
word alphaName_;
//- Name of the phase flux field (optional)
word alphaPhiName_;
//- Name of the mixture flux field (optional)
word phiName_;
//- Name of density field (optional)
word rhoName_;
//- Name of the pressure field (optional)
word pName_;
//- Diffusion coefficient (optional)
scalar D_;
//- Flag to indicate whether a constant, uniform D_ is specified
bool constantD_;
//- Laminar diffusion coefficient (optional)
scalar alphaD_;
//- Turbulent diffusion coefficient (optional)
scalar alphaDt_;
//- Number of corrector iterations (optional)
label nCorr_;
//- Residual volume-fraction
scalar residualAlpha_;
//- Name of field whose schemes are used (optional)
word schemesField_;
//- Run-time selectable finite volume options, e.g. sources, constraints
fv::optionList fvOptions_;
//- The field
volScalarField s_;
//- The field multiplied by the phase fraction
volScalarField alphaS_;
//- Potential field used to generate the phase flux
autoPtr<volScalarField> PhiPtr_;
// Private Member Functions
//- Return the potential field used to generate the phase flux.
// Constructed on demand.
volScalarField& Phi();
//- Return the phase flux. Tries to look it up, and generates it if the
// lookup fails. The generation creates a crude guess for alphaPhi
// then solves a Laplacian to correct the flux to match the time
// derivative of alpha.
tmp<surfaceScalarField> alphaPhi();
//- Return the diffusivity field
tmp<volScalarField> D(const surfaceScalarField& alphaPhi) const;
public:
//- Runtime type information
TypeName("phaseScalarTransport");
// Constructors
//- Construct from Time and dictionary
phaseScalarTransport
(
const word& name,
const Time& runTime,
const dictionary& dict
);
//- Destructor
virtual ~phaseScalarTransport();
// Member Functions
//- Read the settings from the given dictionary
virtual bool read(const dictionary&);
//- Solve for the evolution of the field
virtual bool execute();
//- Do nothing. The field is registered and written automatically.
virtual bool write();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace functionObjects
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //