From 6e8f0dbe761d2b2fe7b15f610ffac4dd8d4e67ce Mon Sep 17 00:00:00 2001
From: sergio <s.ferraris@opencfd.co.uk>
Date: Wed, 4 Dec 2019 11:13:45 -0800
Subject: [PATCH] INT: org integration
1) rPolynomial Eq of State
2) externalForce and softWall in rigidBodyDynamics
INT: Several minor bug fixes plus
---
.../particleTracks/particleTracks.C | 2 +-
.../steadyParticleTracks.C | 4 +-
.../porosityModel/porosityModel.C | 5 +-
.../clouds/Templates/MPPICCloud/MPPICCloud.C | 14 +
.../dragModels/segregated/segregated.C | 30 +-
.../reactionThermo/hRefConstThermos.C | 69 +++++
src/rigidBodyDynamics/Make/files | 1 +
.../restraints/externalForce/externalForce.C | 130 +++++++++
.../restraints/externalForce/externalForce.H | 146 ++++++++++
.../linearAxialAngularSpring.C | 3 +-
.../linearAxialAngularSpring.H | 3 +-
.../restraints/linearDamper/linearDamper.C | 3 +-
.../restraints/linearDamper/linearDamper.H | 3 +-
.../restraints/linearSpring/linearSpring.C | 3 +-
.../restraints/linearSpring/linearSpring.H | 3 +-
.../prescribedRotation/prescribedRotation.C | 3 +-
.../prescribedRotation/prescribedRotation.H | 3 +-
.../restraints/restraint/rigidBodyRestraint.H | 4 +-
.../restraints/softWall/softWall.C | 3 +-
.../restraints/softWall/softWall.H | 3 +-
.../sphericalAngularDamper.C | 3 +-
.../sphericalAngularDamper.H | 3 +-
.../rigidBodyModel/forwardDynamics.C | 5 +-
.../rigidBodyModel/rigidBodyModel.H | 7 +-
.../CrankNicolson/CrankNicolson.C | 2 +-
.../rigidBodySolvers/Newmark/Newmark.C | 2 +-
.../rigidBodySolvers/symplectic/symplectic.C | 2 +-
.../Chung/Chung.C | 18 +-
.../Chung/Chung.H | 2 +-
.../Wallis/Wallis.C | 18 +-
.../Wallis/Wallis.H | 1 +
.../gradientEnergyFvPatchScalarField.C | 7 +
.../gradientEnergyFvPatchScalarField.H | 3 +
.../basic/rhoThermo/rhoThermos.C | 49 ++++
.../mixtures/SpecieMixture/SpecieMixture.C | 23 +-
.../mixtures/SpecieMixture/SpecieMixture.H | 13 +-
.../basicSpecieMixture/basicSpecieMixture.H | 13 +-
.../rhoReactionThermo/rhoReactionThermos.C | 1 +
.../equationOfState/rPolynomial/rPolynomial.C | 65 +++++
.../equationOfState/rPolynomial/rPolynomial.H | 270 ++++++++++++++++++
.../rPolynomial/rPolynomialI.H | 235 +++++++++++++++
.../specie/include/thermoPhysicsTypes.H | 29 ++
42 files changed, 1142 insertions(+), 64 deletions(-)
create mode 100644 src/rigidBodyDynamics/restraints/externalForce/externalForce.C
create mode 100644 src/rigidBodyDynamics/restraints/externalForce/externalForce.H
create mode 100644 src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.C
create mode 100644 src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.H
create mode 100644 src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomialI.H
diff --git a/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C b/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
index cebd58fb637..e247d87021b 100644
--- a/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C
@@ -69,7 +69,7 @@ int main(int argc, char *argv[])
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
- fileName vtkPath(runTime.path()/"VTK");
+ const fileName vtkPath(runTime.rootPath()/runTime.globalCaseName()/"VTK");
mkDir(vtkPath);
Info<< "Scanning times to determine track data for cloud " << cloudName
diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
index 876350c9e4d..02e0b78e182 100644
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
- fileName vtkPath(runTime.path()/"VTK");
+ const fileName vtkPath(runTime.rootPath()/runTime.globalCaseName()/"VTK");
mkDir(vtkPath);
forAll(timeDirs, timeI)
@@ -145,7 +145,7 @@ int main(int argc, char *argv[])
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Time = " << runTime.timeName() << endl;
- fileName vtkTimePath(runTime.path()/"VTK"/runTime.timeName());
+ const fileName vtkTimePath(vtkPath/runTime.timeName());
mkDir(vtkTimePath);
Info<< " Reading particle positions" << endl;
diff --git a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
index f20e9de9065..716e8aaa7f7 100644
--- a/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
+++ b/src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C
@@ -42,12 +42,13 @@ namespace Foam
void Foam::porosityModel::adjustNegativeResistance(dimensionedVector& resist)
{
- scalar maxCmpt = max(0, cmptMax(resist.value()));
+ scalar maxCmpt = cmptMax(resist.value());
if (maxCmpt < 0)
{
FatalErrorInFunction
- << "Negative resistances are invalid, resistance = " << resist
+ << "Cannot have all resistances set to negative, resistance = "
+ << resist
<< exit(FatalError);
}
else
diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
index e687d4b4f99..54549cbe8cf 100644
--- a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
+++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C
@@ -202,6 +202,13 @@ void Foam::MPPICCloud<CloudType>::motion
if (dampingModel_->active())
{
+ if (this->mesh().moving())
+ {
+ FatalErrorInFunction
+ << "MPPIC damping modelling does not support moving meshes."
+ << exit(FatalError);
+ }
+
// update averages
td.updateAverages(cloud);
@@ -226,6 +233,13 @@ void Foam::MPPICCloud<CloudType>::motion
if (packingModel_->active())
{
+ if (this->mesh().moving())
+ {
+ FatalErrorInFunction
+ << "MPPIC packing modelling does not support moving meshes."
+ << exit(FatalError);
+ }
+
// same procedure as for damping
td.updateAverages(cloud);
packingModel_->cacheFields(true);
diff --git a/src/phaseSystemModels/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C b/src/phaseSystemModels/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C
index d92458543fe..3f96b5353e7 100644
--- a/src/phaseSystemModels/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C
+++ b/src/phaseSystemModels/reactingEulerFoam/interfacialModels/dragModels/segregated/segregated.C
@@ -109,7 +109,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
L.primitiveFieldRef() = cbrt(mesh.V());
L.correctBoundaryConditions();
- volScalarField I
+ const volScalarField I
(
alpha1
/max
@@ -118,7 +118,7 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
pair_.phase1().residualAlpha() + pair_.phase2().residualAlpha()
)
);
- volScalarField magGradI
+ const volScalarField magGradI
(
max
(
@@ -127,28 +127,36 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
)
);
- volScalarField muI
+ const volScalarField muI
(
rho1*nu1*rho2*nu2
/(rho1*nu1 + rho2*nu2)
);
- volScalarField muAlphaI
+
+ const volScalarField limitedAlpha1
+ (
+ max(alpha1, pair_.phase1().residualAlpha())
+ );
+
+ const volScalarField limitedAlpha2
+ (
+ max(alpha2, pair_.phase2().residualAlpha())
+ );
+
+ const volScalarField muAlphaI
(
alpha1*rho1*nu1*alpha2*rho2*nu2
- /(
- max(alpha1, pair_.phase1().residualAlpha())*rho1*nu1
- + max(alpha2, pair_.phase2().residualAlpha())*rho2*nu2
- )
+ /(limitedAlpha1*rho1*nu1 + limitedAlpha2*rho2*nu2)
);
- volScalarField ReI
+ const volScalarField ReI
(
pair_.rho()
*pair_.magUr()
- /(magGradI*muI)
+ /(magGradI*limitedAlpha1*limitedAlpha2*muI)
);
- volScalarField lambda(m_*ReI + n_*muAlphaI/muI);
+ const volScalarField lambda(m_*ReI + n_*muAlphaI/muI);
return lambda*sqr(magGradI)*muI;
}
diff --git a/src/phaseSystemModels/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C b/src/phaseSystemModels/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
index feb63aea7f0..8cd582eb718 100644
--- a/src/phaseSystemModels/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
+++ b/src/phaseSystemModels/reactingEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
@@ -33,6 +33,7 @@ License
#include "specie.H"
#include "perfectGas.H"
+#include "rPolynomial.H"
#include "perfectFluid.H"
#include "rhoConst.H"
@@ -84,6 +85,19 @@ constTransport
>
> constRefFluidHThermoPhysics;
+typedef
+constTransport
+<
+ species::thermo
+ <
+ hRefConstThermo
+ <
+ rPolynomial<specie>
+ >,
+ sensibleEnthalpy
+ >
+> constRefrPolFluidHThermoPhysics;
+
typedef
constTransport
<
@@ -110,6 +124,19 @@ constTransport
>
> constRefFluidEThermoPhysics;
+typedef
+constTransport
+<
+ species::thermo
+ <
+ eRefConstThermo
+ <
+ rPolynomial<specie>
+ >,
+ sensibleInternalEnergy
+ >
+> constRefrPolFluidEThermoPhysics;
+
typedef
constTransport
<
@@ -151,6 +178,18 @@ makeThermos
specie
);
+makeThermos
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ constTransport,
+ sensibleEnthalpy,
+ hRefConstThermo,
+ rPolynomial,
+ specie
+);
+
makeThermos
(
rhoThermo,
@@ -190,6 +229,18 @@ makeThermos
specie
);
+makeThermos
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ constTransport,
+ sensibleInternalEnergy,
+ eRefConstThermo,
+ rPolynomial,
+ specie
+);
+
makeThermos
(
rhoThermo,
@@ -235,6 +286,15 @@ makeThermoPhysicsReactionThermos
constRefFluidEThermoPhysics
);
+makeThermoPhysicsReactionThermos
+(
+ rhoThermo,
+ rhoReactionThermo,
+ heRhoThermo,
+ multiComponentMixture,
+ constRefrPolFluidEThermoPhysics
+);
+
makeThermoPhysicsReactionThermos
(
rhoThermo,
@@ -265,6 +325,15 @@ makeThermoPhysicsReactionThermos
constRefFluidHThermoPhysics
);
+makeThermoPhysicsReactionThermos
+(
+ rhoThermo,
+ rhoReactionThermo,
+ heRhoThermo,
+ multiComponentMixture,
+ constRefrPolFluidHThermoPhysics
+);
+
makeThermoPhysicsReactionThermos
(
rhoThermo,
diff --git a/src/rigidBodyDynamics/Make/files b/src/rigidBodyDynamics/Make/files
index 703c03aee83..01067cdab47 100644
--- a/src/rigidBodyDynamics/Make/files
+++ b/src/rigidBodyDynamics/Make/files
@@ -33,6 +33,7 @@ restraints/linearDamper/linearDamper.C
restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
restraints/sphericalAngularDamper/sphericalAngularDamper.C
restraints/prescribedRotation/prescribedRotation.C
+restraints/externalForce/externalForce.C
restraints/softWall/softWall.C
rigidBodyModel/rigidBodyModel.C
diff --git a/src/rigidBodyDynamics/restraints/externalForce/externalForce.C b/src/rigidBodyDynamics/restraints/externalForce/externalForce.C
new file mode 100644
index 00000000000..86f125098d8
--- /dev/null
+++ b/src/rigidBodyDynamics/restraints/externalForce/externalForce.C
@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "externalForce.H"
+#include "rigidBodyModel.H"
+#include "rigidBodyModelState.H"
+#include "OneConstant.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RBD
+{
+namespace restraints
+{
+ defineTypeNameAndDebug(externalForce, 0);
+
+ addToRunTimeSelectionTable
+ (
+ restraint,
+ externalForce,
+ dictionary
+ );
+}
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::RBD::restraints::externalForce::externalForce
+(
+ const word& name,
+ const dictionary& dict,
+ const rigidBodyModel& model
+)
+:
+ restraint(name, dict, model),
+ externalForce_(nullptr)
+{
+ read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::RBD::restraints::externalForce::~externalForce()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+void Foam::RBD::restraints::externalForce::restrain
+(
+ scalarField& tau,
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
+) const
+{
+ const vector force = externalForce_().value(state.t());
+ const vector moment(location_ ^ force);
+
+ if (model_.debug)
+ {
+ Info<< " location " << location_
+ << " force " << force
+ << " moment " << moment
+ << endl;
+ }
+
+ // Accumulate the force for the restrained body
+ fx[bodyIndex_] += spatialVector(moment, force);
+}
+
+
+bool Foam::RBD::restraints::externalForce::read
+(
+ const dictionary& dict
+)
+{
+ restraint::read(dict);
+
+ coeffs_.lookup("location") >> location_;
+
+ externalForce_ = Function1<vector>::New("force", coeffs_);
+
+ return true;
+}
+
+
+void Foam::RBD::restraints::externalForce::write
+(
+ Ostream& os
+) const
+{
+ restraint::write(os);
+
+ os.writeEntry("location", location_);
+
+ externalForce_().writeData(os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/rigidBodyDynamics/restraints/externalForce/externalForce.H b/src/rigidBodyDynamics/restraints/externalForce/externalForce.H
new file mode 100644
index 00000000000..4b5e10365b5
--- /dev/null
+++ b/src/rigidBodyDynamics/restraints/externalForce/externalForce.H
@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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::RBD::restraints::externalForce
+
+Group
+ grpRigidBodyDynamicsRestraints
+
+Description
+ Time-dependent external force restraint using Function1.
+
+Usage
+ Example applying a constant force to the floatingObject:
+ \verbatim
+ restraints
+ {
+ force
+ {
+ type externalForce;
+ body floatingObject;
+ location (0 0 0);
+ force (100 0 0);
+ }
+ }
+ \endverbatim
+
+SourceFiles
+ externalForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RBD_restraints_externalForce_H
+#define RBD_restraints_externalForce_H
+
+#include "rigidBodyRestraint.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RBD
+{
+namespace restraints
+{
+
+/*---------------------------------------------------------------------------*\
+ Class externalForce Declaration
+\*---------------------------------------------------------------------------*/
+
+class externalForce
+:
+ public restraint
+{
+ // Private data
+
+ //- External force (N)
+ autoPtr<Function1<vector>> externalForce_;
+
+ //- Point of application of the force
+ vector location_;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("externalForce");
+
+
+ // Constructors
+
+ //- Construct from components
+ externalForce
+ (
+ const word& name,
+ const dictionary& dict,
+ const rigidBodyModel& model
+ );
+
+ //- Construct and return a clone
+ virtual autoPtr<restraint> clone() const
+ {
+ return autoPtr<restraint>
+ (
+ new externalForce(*this)
+ );
+ }
+
+
+ //- Destructor
+ virtual ~externalForce();
+
+
+ // Member Functions
+
+ //- Accumulate the retraint internal joint forces into the tau field and
+ // external forces into the fx field
+ virtual void restrain
+ (
+ scalarField& tau,
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
+ ) const;
+
+ //- Update properties from given dictionary
+ virtual bool read(const dictionary& dict);
+
+ //- Write
+ virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace restraints
+} // End namespace RBD
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C b/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
index 45b1377320b..8f6e20d832c 100644
--- a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
+++ b/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.C
@@ -77,7 +77,8 @@ Foam::RBD::restraints::linearAxialAngularSpring::~linearAxialAngularSpring()
void Foam::RBD::restraints::linearAxialAngularSpring::restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
vector refDir = rotationTensor(vector(1, 0, 0), axis_) & vector(0, 1, 0);
diff --git a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H b/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
index da313611ea3..632c2abb374 100644
--- a/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
+++ b/src/rigidBodyDynamics/restraints/linearAxialAngularSpring/linearAxialAngularSpring.H
@@ -111,7 +111,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C b/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C
index b21f80110f5..3bfbdfee045 100644
--- a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C
+++ b/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.C
@@ -76,7 +76,8 @@ Foam::RBD::restraints::linearDamper::~linearDamper()
void Foam::RBD::restraints::linearDamper::restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
vector force = -coeff_*model_.v(model_.master(bodyID_)).l();
diff --git a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H b/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H
index 9859cb51ce3..a4481576618 100644
--- a/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H
+++ b/src/rigidBodyDynamics/restraints/linearDamper/linearDamper.H
@@ -102,7 +102,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C b/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C
index e53d1e8a764..148281aa81b 100644
--- a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C
+++ b/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.C
@@ -76,7 +76,8 @@ Foam::RBD::restraints::linearSpring::~linearSpring()
void Foam::RBD::restraints::linearSpring::restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
point attachmentPt = bodyPoint(refAttachmentPt_);
diff --git a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H b/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H
index 396efc567e6..a206b4d6bed 100644
--- a/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H
+++ b/src/rigidBodyDynamics/restraints/linearSpring/linearSpring.H
@@ -115,7 +115,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C
index 1a58af51f98..297dde3d650 100644
--- a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C
+++ b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.C
@@ -78,7 +78,8 @@ Foam::RBD::restraints::prescribedRotation::~prescribedRotation()
void Foam::RBD::restraints::prescribedRotation::restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
vector refDir = rotationTensor(vector(1, 0, 0), axis_) & vector(0, 1, 0);
diff --git a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H
index e3f6f3c758c..7533e46b8ae 100644
--- a/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H
+++ b/src/rigidBodyDynamics/restraints/prescribedRotation/prescribedRotation.H
@@ -123,7 +123,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H b/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H
index 80c4f037336..87ffb662ea2 100644
--- a/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H
+++ b/src/rigidBodyDynamics/restraints/restraint/rigidBodyRestraint.H
@@ -50,6 +50,7 @@ SourceFiles
#include "point.H"
#include "scalarField.H"
#include "runTimeSelectionTables.H"
+#include "rigidBodyModelState.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -164,7 +165,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const = 0;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/restraints/softWall/softWall.C b/src/rigidBodyDynamics/restraints/softWall/softWall.C
index 26a4b360b11..f13b9766864 100644
--- a/src/rigidBodyDynamics/restraints/softWall/softWall.C
+++ b/src/rigidBodyDynamics/restraints/softWall/softWall.C
@@ -77,7 +77,8 @@ Foam::RBD::restraints::softWall::~softWall()
void Foam::RBD::restraints::softWall::restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
diff --git a/src/rigidBodyDynamics/restraints/softWall/softWall.H b/src/rigidBodyDynamics/restraints/softWall/softWall.H
index 7f58b3132e5..86e2067982c 100644
--- a/src/rigidBodyDynamics/restraints/softWall/softWall.H
+++ b/src/rigidBodyDynamics/restraints/softWall/softWall.H
@@ -121,7 +121,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C b/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C
index 9ea6d8bd38a..be8b1fd9847 100644
--- a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C
+++ b/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.C
@@ -76,7 +76,8 @@ Foam::RBD::restraints::sphericalAngularDamper::~sphericalAngularDamper()
void Foam::RBD::restraints::sphericalAngularDamper::restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
vector moment = -coeff_*model_.v(model_.master(bodyID_)).w();
diff --git a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H b/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H
index c05dac59f9b..a80b479786d 100644
--- a/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H
+++ b/src/rigidBodyDynamics/restraints/sphericalAngularDamper/sphericalAngularDamper.H
@@ -103,7 +103,8 @@ public:
virtual void restrain
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const;
//- Update properties from given dictionary
diff --git a/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C b/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C
index 95d854048ef..c0dd3988c2f 100644
--- a/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C
+++ b/src/rigidBodyDynamics/rigidBodyModel/forwardDynamics.C
@@ -34,7 +34,8 @@ License
void Foam::RBD::rigidBodyModel::applyRestraints
(
scalarField& tau,
- Field<spatialVector>& fx
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
) const
{
if (restraints_.empty())
@@ -47,7 +48,7 @@ void Foam::RBD::rigidBodyModel::applyRestraints
DebugInfo << "Restraint " << restraints_[ri].name();
// Accumulate the restraint forces
- restraints_[ri].restrain(tau, fx);
+ restraints_[ri].restrain(tau, fx, state);
}
}
diff --git a/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H b/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H
index 77667385bad..046b5eda743 100644
--- a/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H
+++ b/src/rigidBodyDynamics/rigidBodyModel/rigidBodyModel.H
@@ -340,7 +340,12 @@ public:
//- Apply the restraints and accumulate the internal joint forces
// into the tau field and external forces into the fx field
- void applyRestraints(scalarField& tau, Field<spatialVector>& fx) const;
+ void applyRestraints
+ (
+ scalarField& tau,
+ Field<spatialVector>& fx,
+ const rigidBodyModelState& state
+ ) const;
//- Calculate the joint acceleration qDdot from the joint state q,
// velocity qDot, internal force tau (in the joint frame) and
diff --git a/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C b/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C
index a7542f64190..1a0f0f8f572 100644
--- a/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C
+++ b/src/rigidBodyDynamics/rigidBodySolvers/CrankNicolson/CrankNicolson.C
@@ -74,7 +74,7 @@ void Foam::RBD::rigidBodySolvers::CrankNicolson::solve
// Accumulate the restraint forces
scalarField rtau(tau);
Field<spatialVector> rfx(fx);
- model_.applyRestraints(rtau, rfx);
+ model_.applyRestraints(rtau, rfx, state());
// Calculate the accelerations for the given state and forces
model_.forwardDynamics(state(), rtau, rfx);
diff --git a/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C b/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C
index ad0e5633612..ff7a3c26ddd 100644
--- a/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C
+++ b/src/rigidBodyDynamics/rigidBodySolvers/Newmark/Newmark.C
@@ -81,7 +81,7 @@ void Foam::RBD::rigidBodySolvers::Newmark::solve
// Accumulate the restraint forces
scalarField rtau(tau);
Field<spatialVector> rfx(fx);
- model_.applyRestraints(rtau, rfx);
+ model_.applyRestraints(rtau, rfx, state());
// Calculate the accelerations for the given state and forces
model_.forwardDynamics(state(), rtau, rfx);
diff --git a/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C b/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C
index 16ca4096572..2b2a6c0de4f 100644
--- a/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C
+++ b/src/rigidBodyDynamics/rigidBodySolvers/symplectic/symplectic.C
@@ -83,7 +83,7 @@ void Foam::RBD::rigidBodySolvers::symplectic::solve
// Accumulate the restraint forces
scalarField rtau(tau);
Field<spatialVector> rfx(fx);
- model_.applyRestraints(rtau, rfx);
+ model_.applyRestraints(rtau, rfx, state());
// Calculate the body acceleration for the given state
// and restraint forces
diff --git a/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.C b/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.C
index 54ea1301078..40e640c05eb 100644
--- a/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.C
+++ b/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.C
@@ -54,6 +54,12 @@ Foam::compressibilityModels::Chung::Chung
)
:
barotropicCompressibilityModel(compressibilityProperties, gamma, psiName),
+ pSat_
+ (
+ "pSat",
+ dimPressure,
+ compressibilityProperties_
+ ),
psiv_
(
"psiv",
@@ -66,12 +72,6 @@ Foam::compressibilityModels::Chung::Chung
dimCompressibility,
compressibilityProperties_
),
- rhovSat_
- (
- "rhovSat",
- dimDensity,
- compressibilityProperties_
- ),
rholSat_
(
"rholSat",
@@ -91,8 +91,8 @@ void Foam::compressibilityModels::Chung::correct()
(
sqrt
(
- (rhovSat_/psiv_)
- /((scalar(1) - gamma_)*rhovSat_/psiv_ + gamma_*rholSat_/psil_)
+ pSat_
+ /((scalar(1) - gamma_)*pSat_ + gamma_*rholSat_/psil_)
)
);
@@ -111,9 +111,9 @@ bool Foam::compressibilityModels::Chung::read
{
barotropicCompressibilityModel::read(compressibilityProperties);
+ compressibilityProperties_.readEntry("pSat", pSat_);
compressibilityProperties_.readEntry("psiv", psiv_);
compressibilityProperties_.readEntry("psil", psil_);
- compressibilityProperties_.readEntry("rhovSat", rhovSat_);
compressibilityProperties_.readEntry("rholSat", rholSat_);
return true;
diff --git a/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.H b/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.H
index 570cf098dfd..8244c7095aa 100644
--- a/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.H
+++ b/src/thermophysicalModels/barotropicCompressibilityModel/Chung/Chung.H
@@ -57,10 +57,10 @@ class Chung
{
// Private data
+ dimensionedScalar pSat_;
dimensionedScalar psiv_;
dimensionedScalar psil_;
- dimensionedScalar rhovSat_;
dimensionedScalar rholSat_;
diff --git a/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.C b/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.C
index 1f7a5006b67..ed2a2476a50 100644
--- a/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.C
+++ b/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.C
@@ -54,6 +54,12 @@ Foam::compressibilityModels::Wallis::Wallis
)
:
barotropicCompressibilityModel(compressibilityProperties, gamma, psiName),
+ pSat_
+ (
+ "pSat",
+ dimPressure,
+ compressibilityProperties_
+ ),
psiv_
(
"psiv",
@@ -66,12 +72,7 @@ Foam::compressibilityModels::Wallis::Wallis
dimCompressibility,
compressibilityProperties_
),
- rhovSat_
- (
- "rhovSat",
- dimDensity,
- compressibilityProperties_
- ),
+ rhovSat_(psiv_*pSat_),
rholSat_
(
"rholSat",
@@ -89,7 +90,7 @@ void Foam::compressibilityModels::Wallis::correct()
{
psi_ =
(gamma_*rhovSat_ + (scalar(1) - gamma_)*rholSat_)
- *(gamma_*psiv_/rhovSat_ + (scalar(1) - gamma_)*psil_/rholSat_);
+ *(gamma_/pSat_ + (scalar(1) - gamma_)*psil_/rholSat_);
}
@@ -100,9 +101,10 @@ bool Foam::compressibilityModels::Wallis::read
{
barotropicCompressibilityModel::read(compressibilityProperties);
+ compressibilityProperties_.readEntry("pSat", pSat_);
compressibilityProperties_.readEntry("psiv", psiv_);
compressibilityProperties_.readEntry("psil", psil_);
- compressibilityProperties_.readEntry("rhovSat", rhovSat_);
+ rhovSat_ = psiv_*pSat_;
compressibilityProperties_.readEntry("rholSat", rholSat_);
return true;
diff --git a/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.H b/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.H
index 103fbf6bfae..c8abe7eef86 100644
--- a/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.H
+++ b/src/thermophysicalModels/barotropicCompressibilityModel/Wallis/Wallis.H
@@ -57,6 +57,7 @@ class Wallis
{
// Private data
+ dimensionedScalar pSat_;
dimensionedScalar psiv_;
dimensionedScalar psil_;
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.C b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.C
index cbdea7f7605..7866817302d 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.C
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.C
@@ -119,6 +119,13 @@ void Foam::gradientEnergyFvPatchScalarField::updateCoeffs()
}
+void Foam::gradientEnergyFvPatchScalarField::write(Ostream& os) const
+{
+ fixedGradientFvPatchScalarField::write(os);
+ os.writeEntry("value", *this);
+}
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
diff --git a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.H b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.H
index 586a9cfb19f..bf84e88d0bb 100644
--- a/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.H
+++ b/src/thermophysicalModels/basic/derivedFvPatchFields/gradientEnergy/gradientEnergyFvPatchScalarField.H
@@ -158,6 +158,9 @@ public:
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
+
+ //- Write
+ virtual void write(Ostream&) const;
};
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
index 44c128db354..7789ebbe91e 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermos.C
@@ -33,6 +33,7 @@ License
#include "incompressiblePerfectGas.H"
#include "Boussinesq.H"
#include "rhoConst.H"
+#include "rPolynomial.H"
#include "perfectFluid.H"
#include "PengRobinsonGas.H"
#include "adiabaticPerfectFluid.H"
@@ -122,6 +123,18 @@ makeThermos
specie
);
+makeThermos
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ constTransport,
+ sensibleEnthalpy,
+ hConstThermo,
+ rPolynomial,
+ specie
+);
+
makeThermos
(
rhoThermo,
@@ -146,6 +159,18 @@ makeThermos
specie
);
+makeThermos
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ polynomialTransport,
+ sensibleEnthalpy,
+ hPolynomialThermo,
+ rPolynomial,
+ specie
+);
+
makeThermos
(
rhoThermo,
@@ -328,6 +353,18 @@ makeThermos
specie
);
+makeThermos
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ constTransport,
+ sensibleInternalEnergy,
+ hConstThermo,
+ rPolynomial,
+ specie
+);
+
makeThermos
(
rhoThermo,
@@ -340,6 +377,18 @@ makeThermos
specie
);
+makeThermos
+(
+ rhoThermo,
+ heRhoThermo,
+ pureMixture,
+ constTransport,
+ sensibleInternalEnergy,
+ eConstThermo,
+ rPolynomial,
+ specie
+);
+
makeThermos
(
rhoThermo,
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C
index 9c442b53eb5..c5d85769332 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C
+++ b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C
@@ -60,6 +60,13 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::W
}
+template<class MixtureType>
+Foam::scalar Foam::SpecieMixture<MixtureType>::Hc(const label speciei) const
+{
+ return this->getLocalThermo(speciei).Hc();
+}
+
+
template<class MixtureType>
Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
(
@@ -85,36 +92,38 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::Cv
template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Ha
+Foam::scalar Foam::SpecieMixture<MixtureType>::HE
(
const label speciei,
const scalar p,
const scalar T
) const
{
- return this->getLocalThermo(speciei).Ha(p, T);
+ return this->getLocalThermo(speciei).HE(p, T);
}
template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Hs
+Foam::scalar Foam::SpecieMixture<MixtureType>::Ha
(
const label speciei,
const scalar p,
const scalar T
) const
{
- return this->getLocalThermo(speciei).Hs(p, T);
+ return this->getLocalThermo(speciei).Ha(p, T);
}
template<class MixtureType>
-Foam::scalar Foam::SpecieMixture<MixtureType>::Hc
+Foam::scalar Foam::SpecieMixture<MixtureType>::Hs
(
- const label speciei
+ const label speciei,
+ const scalar p,
+ const scalar T
) const
{
- return this->getLocalThermo(speciei).Hc();
+ return this->getLocalThermo(speciei).Hs(p, T);
}
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H
index 3fb34f531cd..572dff91e21 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H
@@ -87,6 +87,9 @@ public:
//- Molecular weight of the given specie [kg/kmol]
virtual scalar W(const label speciei) const;
+ //- Chemical enthalpy [J/kg]
+ virtual scalar Hc(const label speciei) const;
+
// Per specie thermo properties
@@ -106,6 +109,14 @@ public:
const scalar T
) const;
+ //- Enthalpy/Internal energy [J/kg]
+ virtual scalar HE
+ (
+ const label speciei,
+ const scalar p,
+ const scalar T
+ ) const;
+
//- Absolute enthalpy [J/kg]
virtual scalar Ha
(
@@ -122,8 +133,6 @@ public:
const scalar T
) const;
- //- Chemical enthalpy [J/kg]
- virtual scalar Hc(const label speciei) const;
//- Entropy [J/(kg K)]
virtual scalar S
diff --git a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
index 76194cc4eff..d7db29cc424 100644
--- a/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
+++ b/src/thermophysicalModels/reactionThermo/mixtures/basicSpecieMixture/basicSpecieMixture.H
@@ -89,6 +89,9 @@ public:
//- Molecular weight of the given specie [kg/kmol]
virtual scalar W(const label speciei) const = 0;
+ //- Chemical enthalpy [J/kg]
+ virtual scalar Hc(const label speciei) const = 0;
+
// Per specie thermo properties
@@ -108,6 +111,14 @@ public:
const scalar T
) const = 0;
+ //- Enthalpy/Internal energy [J/kg]
+ virtual scalar HE
+ (
+ const label speciei,
+ const scalar p,
+ const scalar T
+ ) const = 0;
+
//- Absolute enthalpy [J/kg]
virtual scalar Ha
(
@@ -124,8 +135,6 @@ public:
const scalar T
) const = 0;
- //- Chemical enthalpy [J/kg]
- virtual scalar Hc(const label speciei) const = 0;
//- Entropy [J/(kg K)]
virtual scalar S
diff --git a/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermos.C b/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermos.C
index 734daa6cd56..1f09cb42729 100644
--- a/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermos.C
+++ b/src/thermophysicalModels/reactionThermo/rhoReactionThermo/rhoReactionThermos.C
@@ -39,6 +39,7 @@ License
#include "sensibleEnthalpy.H"
#include "thermo.H"
#include "rhoConst.H"
+#include "rPolynomial.H"
#include "perfectFluid.H"
#include "adiabaticPerfectFluid.H"
#include "Boussinesq.H"
diff --git a/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.C b/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.C
new file mode 100644
index 00000000000..a67c6d7ac19
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.C
@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "rPolynomial.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Specie>
+Foam::rPolynomial<Specie>::rPolynomial(const dictionary& dict)
+:
+ Specie(dict),
+ C_(dict.subDict("equationOfState").lookup("C"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Specie>
+void Foam::rPolynomial<Specie>::write(Ostream& os) const
+{
+ Specie::write(os);
+
+ dictionary dict("equationOfState");
+ dict.add("C", C_);
+
+ os << indent << dict.dictName() << dict;
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
+
+template<class Specie>
+Foam::Ostream& Foam::operator<<(Ostream& os, const rPolynomial<Specie>& pf)
+{
+ pf.write(os);
+ return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.H b/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.H
new file mode 100644
index 00000000000..44ad5aa517d
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomial.H
@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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::rPolynomial
+
+Description
+ Reciprocal polynomial equation of state for liquids and solids
+
+ \f[
+ 1/\rho = C_0 + C_1 T + C_2 T^2 - C_3 p - C_4 p T
+ \f]
+
+ This polynomial for the reciprocal of the density provides a much better fit
+ than the equivalent polynomial for the density and has the advantage that it
+ support coefficient mixing to support liquid and solid mixtures in an
+ efficient manner.
+
+Usage
+ \table
+ Property | Description
+ C | Density polynomial coefficients
+ \endtable
+
+ Example of the specification of the equation of state for pure water:
+ \verbatim
+ equationOfState
+ {
+ C (0.001278 -2.1055e-06 3.9689e-09 4.3772e-13 -2.0225e-16);
+ }
+ \endverbatim
+ Note: This fit is based on the small amount of data which is freely
+ available for the range 20-65degC and 1-100bar.
+
+SourceFiles
+ rPolynomialI.H
+ rPolynomial.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rPolynomial_H
+#define rPolynomial_H
+
+#include "autoPtr.H"
+#include "VectorSpace.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class Specie> class rPolynomial;
+
+template<class Specie>
+inline rPolynomial<Specie> operator+
+(
+ const rPolynomial<Specie>&,
+ const rPolynomial<Specie>&
+);
+
+template<class Specie>
+inline rPolynomial<Specie> operator*
+(
+ const scalar,
+ const rPolynomial<Specie>&
+);
+
+template<class Specie>
+inline rPolynomial<Specie> operator==
+(
+ const rPolynomial<Specie>&,
+ const rPolynomial<Specie>&
+);
+
+template<class Specie>
+Ostream& operator<<
+(
+ Ostream&,
+ const rPolynomial<Specie>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+ Class rPolynomial Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Specie>
+class rPolynomial
+:
+ public Specie
+{
+ // Private Data
+
+ class coeffList
+ :
+ public VectorSpace<coeffList, scalar, 5>
+ {
+ public:
+
+ // Constructors
+
+ //- Construct null
+ inline coeffList()
+ {}
+
+ //- Construct from Istream
+ inline coeffList(Istream& is)
+ :
+ VectorSpace<coeffList, scalar, 5>(is)
+ {}
+ };
+
+
+ //- Density coefficients
+ coeffList C_;
+
+
+public:
+
+ // Constructors
+
+ //- Construct from components
+ inline rPolynomial
+ (
+ const Specie& sp,
+ const coeffList& coeffs
+ );
+
+ //- Construct from dictionary
+ rPolynomial(const dictionary& dict);
+
+ //- Construct as named copy
+ inline rPolynomial(const word& name, const rPolynomial&);
+
+ //- Construct and return a clone
+ inline autoPtr<rPolynomial> clone() const;
+
+ // Selector from dictionary
+ inline static autoPtr<rPolynomial> New(const dictionary& dict);
+
+
+ // Member Functions
+
+ //- Return the instantiated type name
+ static word typeName()
+ {
+ return "rPolynomial<" + word(Specie::typeName_()) + '>';
+ }
+
+
+ // Fundamental properties
+
+ //- Is the equation of state is incompressible i.e. rho != f(p)
+ static const bool incompressible = false;
+
+ //- Is the equation of state is isochoric i.e. rho = const
+ static const bool isochoric = false;
+
+ //- Return density [kg/m^3]
+ inline scalar rho(scalar p, scalar T) const;
+
+ //- Return enthalpy departure [J/kg]
+ inline scalar H(const scalar p, const scalar T) const;
+
+ //- Return Cp departure [J/(kg K]
+ inline scalar Cp(scalar p, scalar T) const;
+
+ //- Return internal energy departure [J/kg]
+ inline scalar E(const scalar p, const scalar T) const;
+
+ //- Return Cv departure [J/(kg K]
+ inline scalar Cv(scalar p, scalar T) const;
+
+ //- Return entropy [J/kg/K]
+ inline scalar S(const scalar p, const scalar T) const;
+
+ //- Return compressibility [s^2/m^2]
+ inline scalar psi(scalar p, scalar T) const;
+
+ //- Return compression factor []
+ inline scalar Z(scalar p, scalar T) const;
+
+ //- Return (Cp - Cv) [J/(kg K]
+ inline scalar CpMCv(scalar p, scalar T) const;
+
+
+ // IO
+
+ //- Write to Ostream
+ void write(Ostream& os) const;
+
+
+ // Member Operators
+
+ inline void operator+=(const rPolynomial&);
+ inline void operator*=(const scalar);
+
+
+ // Friend operators
+
+ friend rPolynomial operator+ <Specie>
+ (
+ const rPolynomial&,
+ const rPolynomial&
+ );
+
+ friend rPolynomial operator* <Specie>
+ (
+ const scalar s,
+ const rPolynomial&
+ );
+
+ friend rPolynomial operator== <Specie>
+ (
+ const rPolynomial&,
+ const rPolynomial&
+ );
+
+
+ // Ostream Operator
+
+ friend Ostream& operator<< <Specie>
+ (
+ Ostream&,
+ const rPolynomial&
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "rPolynomialI.H"
+
+#ifdef NoRepository
+ #include "rPolynomial.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomialI.H b/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomialI.H
new file mode 100644
index 00000000000..0f472363078
--- /dev/null
+++ b/src/thermophysicalModels/specie/equationOfState/rPolynomial/rPolynomialI.H
@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "rPolynomial.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::rPolynomial<Specie>::rPolynomial
+(
+ const Specie& sp,
+ const coeffList& coeffs
+)
+:
+ Specie(sp),
+ C_(coeffs)
+{}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::rPolynomial<Specie>::rPolynomial
+(
+ const word& name,
+ const rPolynomial<Specie>& rp
+)
+:
+ Specie(name, rp),
+ C_(rp.C_)
+{}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::rPolynomial<Specie>>
+Foam::rPolynomial<Specie>::clone() const
+{
+ return autoPtr<rPolynomial<Specie>>::New(*this);
+}
+
+
+template<class Specie>
+inline Foam::autoPtr<Foam::rPolynomial<Specie>>
+Foam::rPolynomial<Specie>::New
+(
+ const dictionary& dict
+)
+{
+ return autoPtr<rPolynomial<Specie>>::New(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::rho(scalar p, scalar T) const
+{
+ return 1/(C_[0] + (C_[1] + C_[2]*T - C_[4]*p)*T - C_[3]*p);
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::H(scalar p, scalar T) const
+{
+ return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::Cp(scalar p, scalar T) const
+{
+ return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::E(scalar p, scalar T) const
+{
+ return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::Cv(scalar p, scalar T) const
+{
+ return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::S(scalar p, scalar T) const
+{
+ return 0;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::psi(scalar p, scalar T) const
+{
+ return sqr(rho(p, T))*(C_[3] + C_[4]*T);
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::Z(scalar p, scalar T) const
+{
+ return 1;
+}
+
+
+template<class Specie>
+inline Foam::scalar Foam::rPolynomial<Specie>::CpMCv(scalar p, scalar T) const
+{
+ return 0;
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
+
+template<class Specie>
+inline void Foam::rPolynomial<Specie>::operator+=
+(
+ const rPolynomial<Specie>& rp
+)
+{
+ const scalar Y1 = this->Y();
+ Specie::operator+=(rp);
+
+ if (mag(this->Y()) > SMALL)
+ {
+ C_ = (Y1*C_ + rp.Y()*rp.C_)/this->Y();
+ }
+}
+
+
+template<class Specie>
+inline void Foam::rPolynomial<Specie>::operator*=(const scalar s)
+{
+ Specie::operator*=(s);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
+
+template<class Specie>
+inline Foam::rPolynomial<Specie> Foam::operator+
+(
+ const rPolynomial<Specie>& rp1,
+ const rPolynomial<Specie>& rp2
+)
+{
+ Specie sp
+ (
+ static_cast<const Specie&>(rp1)
+ + static_cast<const Specie&>(rp2)
+ );
+
+ if (mag(sp.Y()) < SMALL)
+ {
+ return rPolynomial<Specie>
+ (
+ sp,
+ rp1.C_
+ );
+ }
+
+ return rPolynomial<Specie>
+ (
+ sp,
+ (rp1.Y()*rp1.C_ + rp2.Y()*rp2.C_)/sp.Y()
+ );
+}
+
+
+template<class Specie>
+inline Foam::rPolynomial<Specie> Foam::operator*
+(
+ const scalar s,
+ const rPolynomial<Specie>& rp
+)
+{
+ return rPolynomial<Specie>
+ (
+ s*static_cast<const Specie&>(rp),
+ rp.C_
+ );
+}
+
+
+template<class Specie>
+inline Foam::rPolynomial<Specie> Foam::operator==
+(
+ const rPolynomial<Specie>& rp1,
+ const rPolynomial<Specie>& rp2
+)
+{
+ Specie sp
+ (
+ static_cast<const Specie&>(rp1)
+ == static_cast<const Specie&>(rp2)
+ );
+
+ return rPolynomial<Specie>
+ (
+ sp,
+ (rp1.Y()*rp1.C_ - rp2.Y()*rp2.C_)/sp.Y()
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H b/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H
index 2932ad9efba..de1b4e419c1 100644
--- a/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H
+++ b/src/thermophysicalModels/specie/include/thermoPhysicsTypes.H
@@ -38,6 +38,7 @@ Description
#include "specie.H"
#include "perfectGas.H"
#include "incompressiblePerfectGas.H"
+#include "rPolynomial.H"
#include "perfectFluid.H"
#include "adiabaticPerfectFluid.H"
#include "rhoConst.H"
@@ -155,6 +156,20 @@ namespace Foam
>
constFluidHThermoPhysics;
+ typedef
+ constTransport
+ <
+ species::thermo
+ <
+ hConstThermo
+ <
+ rPolynomial<specie>
+ >,
+ sensibleEnthalpy
+ >
+ >
+ constrPolFluidHThermoPhysics;
+
typedef
constTransport
<
@@ -280,6 +295,20 @@ namespace Foam
>
constFluidEThermoPhysics;
+ typedef
+ constTransport
+ <
+ species::thermo
+ <
+ eConstThermo
+ <
+ perfectFluid<specie>
+ >,
+ sensibleInternalEnergy
+ >
+ >
+ constrPolFluidEThermoPhysics;
+
typedef
constTransport
<
--
GitLab