ENH: Added new enthalpyPorosity fvOption, e.g. for windshield defrost

This commit is contained in:
andy
2014-03-05 10:30:36 +00:00
parent 5f6f5eb453
commit efb367a4ca
4 changed files with 749 additions and 0 deletions

View File

@ -12,6 +12,8 @@ $(generalSources)/semiImplicitSource/semiImplicitSource.C
derivedSources=sources/derived
$(derivedSources)/actuationDiskSource/actuationDiskSource.C
$(derivedSources)/enthalpyPorositySource/enthalpyPorositySource.C
$(derivedSources)/enthalpyPorositySource/enthalpyPorositySourceIO.C
$(derivedSources)/explicitPorositySource/explicitPorositySource.C
$(derivedSources)/MRFSource/MRFSource.C
$(derivedSources)/pressureGradientExplicitSource/pressureGradientExplicitSource.C

View File

@ -0,0 +1,428 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 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 "enthalpyPorositySource.H"
#include "fvMatrices.H"
#include "specie.H"
#include "basicThermo.H"
#include "uniformDimensionedFields.H"
#include "fixedValueFvPatchFields.H"
#include "zeroGradientFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
namespace Foam
{
template<>
const char* NamedEnum
<
fv::enthalpyPorositySource::thermoMode,
2
>::names[] =
{
"thermo",
"lookup"
};
namespace fv
{
defineTypeNameAndDebug(enthalpyPorositySource, 0);
addToRunTimeSelectionTable
(
option,
enthalpyPorositySource,
dictionary
);
}
}
const Foam::NamedEnum<Foam::fv::enthalpyPorositySource::thermoMode, 2>
Foam::fv::enthalpyPorositySource::thermoModeTypeNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::wordList Foam::fv::enthalpyPorositySource::alpha1BoundaryTypes() const
{
const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
wordList bTypes(T.boundaryField().size());
forAll(bTypes, patchI)
{
const fvPatchField<scalar>& pf = T.boundaryField()[patchI];
if (isA<fixedValueFvPatchScalarField>(pf))
{
bTypes[patchI] = fixedValueFvPatchScalarField::typeName;
}
else
{
bTypes[patchI] = zeroGradientFvPatchScalarField::typeName;
}
}
return bTypes;
}
bool Foam::fv::enthalpyPorositySource::solveField(const word& fieldName) const
{
bool result = true;
switch (mode_)
{
case mdThermo:
{
const basicThermo& thermo =
mesh_.lookupObject<basicThermo>("thermophysicalProperties");
if (fieldName != thermo.he().name())
{
result = false;
}
break;
}
case mdLookup:
{
if (fieldName != TName_)
{
result = false;
}
break;
}
default:
{
FatalErrorIn
(
"bool Foam::fv::enthalpyPorositySource::solveField"
"("
"const word&"
") const"
)
<< "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
<< abort(FatalError);
}
}
return result;
}
Foam::tmp<Foam::volScalarField> Foam::fv::enthalpyPorositySource::rho() const
{
switch (mode_)
{
case mdThermo:
{
const basicThermo& thermo =
mesh_.lookupObject<basicThermo>("thermophysicalProperties");
return thermo.rho();
break;
}
case mdLookup:
{
if (rhoName_ == "rhoRef")
{
scalar rhoRef(readScalar(coeffs_.lookup("rhoRef")));
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
name_ + ":rho",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("rho", dimDensity, rhoRef),
zeroGradientFvPatchScalarField::typeName
)
);
}
else
{
return mesh_.lookupObject<volScalarField>(rhoName_);
}
break;
}
default:
{
FatalErrorIn
(
"Foam::tmp<Foam::volScalarField> "
"Foam::fv::enthalpyPorositySource::rho() const"
)
<< "Unhandled thermo mode: " << thermoModeTypeNames_[mode_]
<< abort(FatalError);
}
}
return tmp<volScalarField>(NULL);
}
Foam::vector Foam::fv::enthalpyPorositySource::g() const
{
if (mesh_.foundObject<uniformDimensionedVectorField>("g"))
{
const uniformDimensionedVectorField& value =
mesh_.lookupObject<uniformDimensionedVectorField>("g");
return value.value();
}
else
{
return coeffs_.lookup("g");
}
}
void Foam::fv::enthalpyPorositySource::update()
{
if (curTimeIndex_ == mesh_.time().timeIndex())
{
return;
}
const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
forAll(cells_, i)
{
label cellI = cells_[i];
scalar alpha1New = 0.0;
scalar Tc = T[cellI];
if (Tc > Tliquidus_)
{
alpha1New = 1.0;
}
else if (Tc < Tsolidus_)
{
alpha1New = 0.0;
}
else
{
// lever rule
alpha1New = (Tc - Tsolidus_)/(Tliquidus_ - Tsolidus_);
}
alpha1New = (1.0 - relax_)*alpha1_[cellI] + relax_*alpha1New;
dAlpha1_[i] = alpha1New - alpha1_[cellI];
alpha1_[cellI] = alpha1New;
curTimeIndex_ = mesh_.time().timeIndex();
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::fv::enthalpyPorositySource::enthalpyPorositySource
(
const word& sourceName,
const word& modelType,
const dictionary& dict,
const fvMesh& mesh
)
:
option(sourceName, modelType, dict, mesh),
Tliquidus_(readScalar(coeffs_.lookup("Tliquidus"))),
Tsolidus_(readScalar(coeffs_.lookup("Tsolidus"))),
L_(readScalar(coeffs_.lookup("L"))),
relax_(coeffs_.lookupOrDefault("relax", 0.9)),
mode_(thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"))),
rhoName_(coeffs_.lookupOrDefault<word>("rhoName", "rho")),
TName_(coeffs_.lookupOrDefault<word>("TName", "T")),
UName_(coeffs_.lookupOrDefault<word>("UName", "U")),
Cu_(coeffs_.lookupOrDefault<scalar>("Cu", 100000)),
q_(coeffs_.lookupOrDefault("q", 0.001)),
beta_(readScalar(coeffs_.lookup("beta"))),
alpha1_
(
IOobject
(
name_ + ":alpha1",
mesh.time().timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("alpha1", dimless, 0.0),
alpha1BoundaryTypes()
),
dAlpha1_(cells_.size(), 0.0),
curTimeIndex_(-1)
{
fieldNames_.setSize(1, "source");
applied_.setSize(1, false);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::fv::enthalpyPorositySource::alwaysApply() const
{
return true;
}
void Foam::fv::enthalpyPorositySource::addSup
(
fvMatrix<scalar>& eqn,
const label fieldI
)
{
if (!solveField(eqn.psi().name()))
{
return;
}
if (debug)
{
Info<< type() << ": applying source to " << eqn.psi().name() << endl;
}
update();
volScalarField dH
(
IOobject
(
name_ + ":dH",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("dH", dimEnergy/dimMass, 0.0)
);
scalarField& dHi = dH.internalField();
forAll(cells_, i)
{
label cellI = cells_[i];
dHi[cellI] = L_*dAlpha1_[i];
}
const volScalarField rho(this->rho());
const surfaceScalarField& phi =
mesh_.lookupObject<surfaceScalarField>("phi");
dimensionedScalar rhoScale("rhoScale", dimless, 1.0);
if
(
(phi.dimensions() == dimVolume/dimTime)
&& (rho.dimensions() == dimDensity)
)
{
rhoScale.dimensions() /= dimDensity;
}
// contributions added to rhs
if (eqn.psi().dimensions() == dimTemperature)
{
dimensionedScalar Cp
(
"Cp",
dimEnergy/dimMass/dimTemperature,
readScalar(coeffs_.lookup("CpRef"))
);
eqn -=
fvc::ddt((rho*rhoScale*dH/Cp)())
+ fvc::div(phi*fvc::interpolate(dH/Cp));
}
else
{
eqn -=
fvc::ddt((rho*rhoScale*dH)())
+ fvc::div(phi*fvc::interpolate(dH));
}
}
void Foam::fv::enthalpyPorositySource::addSup
(
fvMatrix<vector>& eqn,
const label fieldI
)
{
const volVectorField& U = eqn.psi();
if (U.name() != UName_)
{
return;
}
if (debug)
{
Info<< type() << ": applying source to " << UName_ << endl;
}
update();
scalarField& Sp = eqn.diag();
vectorField& Su = eqn.source();
const scalarField& V = mesh_.V();
const volScalarField rho(this->rho());
const volScalarField& T = mesh_.lookupObject<volScalarField>(TName_);
vector g = this->g();
forAll(cells_, i)
{
label cellI = cells_[i];
scalar Vc = V[cellI];
scalar Tc = T[cellI];
scalar rhoc = rho[cellI];
scalar alpha1c = alpha1_[cellI];
Sp[cellI] -= Vc*rhoc*Cu_*sqr(1.0 - alpha1c)/(pow3(alpha1c) + q_);
Su[cellI] += Vc*rhoc*g*beta_*max(0, (Tc - Tsolidus_));
}
}
// ************************************************************************* //

View File

@ -0,0 +1,250 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 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::fv::enthalpyPorositySource
Description
This source is designed to model the effect of solidification and melting
processes, e.g. windhield defrosting. The phase change occurs between the
limiting temperatures of Tliquidus and Tsolidus.
The presence of the solid phase in the flow field is incorporated into the
model as a momentum porosity contribution; the energy associated with the
phase change is added as an enthalpy contribution.
Based on the references:
1. V.R. Voller and C. Prakash, A fixed grid numerical modelling
methodology for convection-diffusion mushy phase-change problems,
Int. J. Heat Mass Transfer 30(8):17091719, 1987.
2. C.R. Swaminathan. and V.R. Voller, A general enthalpy model for
modeling solidification processes, Metallurgical Transactions
23B:651664, 1992.
\heading Source usage
Example usage:
\verbatim
enthalpyPorositySourceCoeffs
{
type enthalpyPorositySource;
active on;
selectionMode cellZone;
cellZone iceZone;
enthalpyPorositySourceCoeffs
{
Tliquidus 288;
Tsolidus 268;
L 334000;
thermoMode thermo;
beta 50e-6;
// only for incompressible solvers:
// rhoName rhoRef;
// rhoRef 1;
// only for solvers that do not define gravity:
g (0 -9.81 0);
}
}
\endverbatim
Where:
\table
Property | Description | Required | Default value
Tliquidus | Temperature when liquid [K] | yes |
Tsolidus | Temperature when solid [K] | yes |
L | Latent heat of fusion [J/kg] | yes |
thermoMode | Thermo mode [thermo|lookup] | yes |
beta | Thermal expansion coefficient [1/K] | yes |
rhoName | Name of density field | no | rho
rhoRef | Reference density | no |
g | Accelerartion due to gravity | no |
\endtable
SourceFiles
enthalpyPorositySource.C
enthalpyPorositySourceIO.C
\*---------------------------------------------------------------------------*/
#ifndef enthalpyPorositySource_H
#define enthalpyPorositySource_H
#include "fvMesh.H"
#include "volFields.H"
#include "fvOption.H"
#include "NamedEnum.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
/*---------------------------------------------------------------------------*\
Class enthalpyPorositySource Declaration
\*---------------------------------------------------------------------------*/
class enthalpyPorositySource
:
public option
{
public:
enum thermoMode
{
mdThermo,
mdLookup
};
static const NamedEnum<thermoMode, 2> thermoModeTypeNames_;
private:
// Private data
//- Solidification initiated when T > Tliquidus_ [K]
scalar Tliquidus_;
//- Solidification temperature [K]
scalar Tsolidus_;
//- Latent heat of fusion [J/kg]
scalar L_;
//- Phase fraction under-relaxation coefficient
scalar relax_;
//- Thermodynamics mode
thermoMode mode_;
//- Name of density field - default = "rho" (optional)
word rhoName_;
//- Name of temperature field - default = "T" (optional)
word TName_;
//- Name of velocity field - default = "U" (optional)
word UName_;
//- Mushy region momentum sink coefficient [1/s]; default = 10^5
scalar Cu_;
//- Coefficient used in porosity calc - default = 0.001
scalar q_;
//- Thermal expansion coefficient [1/K]
scalar beta_;
//- Phase fraction indicator field
volScalarField alpha1_;
//- Change in phase fraction indicator field
scalarField dAlpha1_;
//- Current time index (used for updating)
label curTimeIndex_;
// Private Member Functions
//- Return the list of alpha1 boundary types
wordList alpha1BoundaryTypes() const;
//- Flag to indicate whether to solve for given field
bool solveField(const word& fieldName) const;
//- Return the density field
tmp<volScalarField> rho() const;
//- Return the gravity vector
vector g() const;
//- Update the model
void update();
//- Disallow default bitwise copy construct
enthalpyPorositySource(const enthalpyPorositySource&);
//- Disallow default bitwise assignment
void operator=(const enthalpyPorositySource&);
public:
//- Runtime type information
TypeName("enthalpyPorositySource");
// Constructors
//- Construct from explicit source name and mesh
enthalpyPorositySource
(
const word& sourceName,
const word& modelType,
const dictionary& dict,
const fvMesh& mesh
);
// Member Functions
//- Flag to bypass the apply flag list checking
virtual bool alwaysApply() const;
// Evaluate
//- Add explicit contribution to enthalpy equation
virtual void addSup(fvMatrix<scalar>& eqn, const label fieldI);
//- Add implicit contribution to momentum equation
virtual void addSup(fvMatrix<vector>& eqn, const label fieldI);
// I-O
//- Write the source properties
virtual void writeData(Ostream&) const;
//- Read source dictionary
virtual bool read(const dictionary& dict);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,69 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2014 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 "enthalpyPorositySource.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::fv::enthalpyPorositySource::writeData(Ostream& os) const
{
os << indent << name_ << endl;
dict_.write(os);
}
bool Foam::fv::enthalpyPorositySource::read(const dictionary& dict)
{
if (option::read(dict))
{
coeffs_.lookup("Tliquidus") >> Tliquidus_;
coeffs_.lookup("Tsolidus") >> Tsolidus_;
coeffs_.lookup("L") >> L_;
coeffs_.readIfPresent("relax", relax_);
mode_ = thermoModeTypeNames_.read(coeffs_.lookup("thermoMode"));
coeffs_.readIfPresent("rhoName", rhoName_);
coeffs_.readIfPresent("TName", TName_);
coeffs_.readIfPresent("UName", UName_);
coeffs_.readIfPresent("Cu", Cu_);
coeffs_.readIfPresent("q", q_);
coeffs_.readIfPresent("beta", beta_);
return true;
}
else
{
return false;
}
return false;
}
// ************************************************************************* //