Commit 1a294d02 authored by Henry's avatar Henry
Browse files

Thermo: move the he BC set functions into basicThermo

parent 5cd985d7
derivedFvPatchFields/wallHeatTransfer/wallHeatTransferFvPatchScalarField.C
compressibleInterFoam.C
EXE = $(FOAM_APPBIN)/compressibleInterFoam
EXE_INC = \
-ItwoPhaseThermo \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-IphaseEquationsOfState/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-ltwoPhaseThermo \
-lfluidThermophysicalModels \
-lspecie \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
-lphaseEquationsOfState \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
......
{
volScalarField kByCv
(
"kByCv",
(alpha1*k1/Cv1 + alpha2*k2/Cv2)
+ (alpha1*rho1 + alpha2*rho2)*turbulence->nut()
);
solve
(
fvm::ddt(rho, T)
+ fvm::div(rhoPhi, T)
- fvm::laplacian(kByCv, T)
+ p*fvc::div(phi)*(alpha1/Cv1 + alpha2/Cv2)
- fvm::laplacian(thermo.alphaEff(rho*turbulence->nut()), T)
+ (
p*fvc::div(phi)
+ fvc::ddt(rho, K) + fvc::div(rhoPhi, K)
)*(alpha1/thermo.thermo1().Cv() + alpha2/thermo.thermo2().Cv())
);
// Update compressibilities
psi1 = eos1->psi(p, T);
psi2 = eos2->psi(p, T);
thermo.correct();
}
......@@ -22,4 +22,6 @@
) * mesh.magSf()
)
);
K = 0.5*magSqr(U);
}
EXE_INC = \
-I.. \
-I../twoPhaseThermo \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
......@@ -11,6 +13,9 @@ EXE_INC = \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude
EXE_LIBS = \
-ltwoPhaseThermo \
-lfluidThermophysicalModels \
-lspecie \
-linterfaceProperties \
-ltwoPhaseInterfaceProperties \
-lincompressibleTransportModels \
......
......@@ -43,7 +43,7 @@ Description
#include "subCycle.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "phaseEquationOfState.H"
#include "twoPhaseThermo.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
......
......@@ -38,9 +38,10 @@ Description
#include "fvCFD.H"
#include "MULES.H"
#include "subCycle.H"
#include "rhoThermo.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "phaseEquationOfState.H"
#include "twoPhaseThermo.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
......
......@@ -28,34 +28,6 @@
#include "createPhi.H"
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
p_rgh
);
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi);
......@@ -69,91 +41,31 @@
scalar(1) - alpha1
);
dimensionedScalar k1
(
"k",
dimensionSet(1, 1, -3, -1, 0),
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("k")
);
dimensionedScalar k2
(
"k",
dimensionSet(1, 1, -3, -1, 0),
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("k")
);
Info<< "Reading thermophysical properties\n" << endl;
dimensionedScalar Cv1
(
"Cv",
dimensionSet(0, 2, -2, -1, 0),
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
).lookup("Cv")
);
dimensionedScalar Cv2
(
"Cv",
dimensionSet(0, 2, -2, -1, 0),
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
).lookup("Cv")
);
autoPtr<phaseEquationOfState> eos1
(
phaseEquationOfState::New
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase1Name()
)
)
);
twoPhaseThermo thermo(mesh, alpha1, alpha2);
autoPtr<phaseEquationOfState> eos2
(
phaseEquationOfState::New
(
twoPhaseProperties.subDict
(
twoPhaseProperties.phase2Name()
)
)
);
volScalarField& rho = thermo.rho();
volScalarField& p = thermo.p();
volScalarField& T = thermo.T();
const volScalarField& rho1 = thermo.thermo1().rho();
const volScalarField& psi1 = thermo.thermo1().psi();
const volScalarField& rho2 = thermo.thermo2().rho();
const volScalarField& psi2 = thermo.thermo2().psi();
volScalarField psi1
(
IOobject
(
"psi1",
runTime.timeName(),
mesh
),
eos1->psi(p, T)
);
psi1.oldTime();
// volScalarField rho
// (
// IOobject
// (
// "rho",
// runTime.timeName(),
// mesh,
// IOobject::READ_IF_PRESENT,
// IOobject::AUTO_WRITE
// ),
// alpha1*rho1 + alpha2*rho2
// );
volScalarField psi2
(
IOobject
(
"psi2",
runTime.timeName(),
mesh
),
eos2->psi(p, T)
);
psi2.oldTime();
dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
......@@ -161,23 +73,6 @@
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());
volScalarField rho1("rho1", eos1->rho(p, T));
volScalarField rho2("rho2", eos2->rho(p, T));
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
alpha1*rho1 + alpha2*rho2
);
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
......@@ -207,3 +102,19 @@
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
{
rho1 = eos1->rho(p, T);
rho2 = eos2->rho(p, T);
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
......@@ -91,9 +88,15 @@
p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
rho1 = eos1->rho(p, T);
rho2 = eos2->rho(p, T);
// thermo.correct();
Info<< "max(U) " << max(mag(U)).value() << endl;
Info<< "min(p_rgh) " << min(p_rgh).value() << endl;
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
}
twoPhaseThermo.C
LIB = $(FOAM_LIBBIN)/libtwoPhaseThermo
EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
LIB_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lfiniteVolume
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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 "twoPhaseThermo.H"
#include "gradientEnergyFvPatchScalarField.H"
#include "mixedEnergyFvPatchScalarField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(twoPhaseThermo, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::twoPhaseThermo::twoPhaseThermo
(
const fvMesh& mesh,
const volScalarField& alpha1,
const volScalarField& alpha2,
const word& phaseName
)
:
rhoThermo(mesh, phaseName),
phaseName1_("1"),
phaseName2_("2"),
alpha1_(alpha1),
alpha2_(alpha2),
thermo1_(rhoThermo::New(mesh, phaseName1_)),
thermo2_(rhoThermo::New(mesh, phaseName2_)),
he_
(
IOobject
(
phasePropertyName
(
"he"
),
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimEnergy/dimMass,
heBoundaryTypes(),
heBoundaryBaseTypes()
)
{
thermo1_->validate("phaseModel 1", "e");
thermo2_->validate("phaseModel 2", "e");
rho_ = alpha1_*thermo1_->rho() + alpha2_*thermo2_->rho();
he_ =
(
alpha1_*thermo1_->rho()*thermo1_->he()
+ alpha2_*thermo2_->rho()*thermo2_->he()
)/rho_;
volScalarField::GeometricBoundaryField& hbf = he_.boundaryField();
forAll(hbf, patchi)
{
if (isA<gradientEnergyFvPatchScalarField>(hbf[patchi]))
{
refCast<gradientEnergyFvPatchScalarField>(hbf[patchi]).gradient()
= hbf[patchi].fvPatchField::snGrad();
}
else if (isA<mixedEnergyFvPatchScalarField>(hbf[patchi]))
{
refCast<mixedEnergyFvPatchScalarField>(hbf[patchi]).refGrad()
= hbf[patchi].fvPatchField::snGrad();
}
}
correct();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::twoPhaseThermo::~twoPhaseThermo()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::twoPhaseThermo::correct()
{
thermo1_->he() = thermo1_->he(p_, T_);
thermo1_->correct();
thermo2_->he() = thermo2_->he(p_, T_);
thermo2_->correct();
psi_ = alpha1_*thermo1_->psi() + alpha2_*thermo2_->psi();
mu_ = alpha1_*thermo1_->mu() + alpha2_*thermo2_->mu();
alpha_ = alpha1_*thermo1_->alpha() + alpha2_*thermo2_->alpha();
}
bool Foam::twoPhaseThermo::incompressible() const
{
return thermo1_->incompressible() && thermo2_->incompressible();
}
bool Foam::twoPhaseThermo::isochoric() const
{
return thermo1_->isochoric() && thermo2_->isochoric();
}
Foam::tmp<Foam::volScalarField> Foam::twoPhaseThermo::he
(
const volScalarField& p,
const volScalarField& T
) const
{
return alpha1_*thermo1_->he(p, T) + alpha2_*thermo2_->he(p, T);
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::he
(
const scalarField& p,
const scalarField& T,
const labelList& cells
) const
{
return
scalarField(alpha1_, cells)*thermo1_->he(p, T, cells)
+ scalarField(alpha2_, cells)*thermo2_->he(p, T, cells);
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::he
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
{
return
alpha1_.boundaryField()[patchi]*thermo1_->he(p, T, patchi)
+ alpha2_.boundaryField()[patchi]*thermo2_->he(p, T, patchi);
}
Foam::tmp<Foam::volScalarField> Foam::twoPhaseThermo::hc() const
{
return alpha1_*thermo1_->hc() + alpha2_*thermo2_->hc();
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::THE
(
const scalarField& h,
const scalarField& p,
const scalarField& T0,
const labelList& cells
) const
{
notImplemented("Foam::twoPhaseThermo::THE");
return T0;
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::THE
(
const scalarField& h,
const scalarField& p,
const scalarField& T0,
const label patchi
) const
{
notImplemented("Foam::twoPhaseThermo::THE");
return T0;
}
Foam::tmp<Foam::volScalarField> Foam::twoPhaseThermo::Cp() const
{
return alpha1_*thermo1_->Cp() + alpha2_*thermo2_->Cp();
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::Cp
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
{
return
alpha1_.boundaryField()[patchi]*thermo1_->Cp(p, T, patchi)
+ alpha2_.boundaryField()[patchi]*thermo2_->Cp(p, T, patchi);
}
Foam::tmp<Foam::volScalarField> Foam::twoPhaseThermo::Cv() const
{
return alpha1_*thermo1_->Cv() + alpha2_*thermo2_->Cv();
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::Cv
(
const scalarField& p,
const scalarField& T,
const label patchi
) const
{
return
alpha1_.boundaryField()[patchi]*thermo1_->Cv(p, T, patchi)
+ alpha2_.boundaryField()[patchi]*thermo2_->Cv(p, T, patchi);
}
Foam::tmp<Foam::volScalarField> Foam::twoPhaseThermo::gamma() const
{
return alpha1_*thermo1_->gamma() + alpha2_*thermo2_->gamma();
}
Foam::tmp<Foam::scalarField> Foam::twoPhaseThermo::gamma
(