From 6e8f0dbe761d2b2fe7b15f610ffac4dd8d4e67ce Mon Sep 17 00:00:00 2001 From: sergio 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 | 32 ++- .../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 | 29 +- .../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, 1146 insertions(+), 68 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 cebd58fb63..e247d87021 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 876350c9e4..02e0b78e18 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 f20e9de906..716e8aaa7f 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 e687d4b4f9..54549cbe8c 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::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::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 d92458543f..3f96b5353e 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::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::dragModels::segregated::K() const pair_.phase1().residualAlpha() + pair_.phase2().residualAlpha() ) ); - volScalarField magGradI + const volScalarField magGradI ( max ( @@ -127,28 +127,36 @@ Foam::tmp Foam::dragModels::segregated::K() const ) ); - volScalarField muI + const volScalarField muI ( rho1*nu1*rho2*nu2 /(rho1*nu1 + rho2*nu2) ); - volScalarField muAlphaI + + const volScalarField limitedAlpha1 ( - alpha1*rho1*nu1*alpha2*rho2*nu2 - /( - max(alpha1, pair_.phase1().residualAlpha())*rho1*nu1 - + max(alpha2, pair_.phase2().residualAlpha())*rho2*nu2 - ) + max(alpha1, pair_.phase1().residualAlpha()) ); - volScalarField ReI + const volScalarField limitedAlpha2 + ( + max(alpha2, pair_.phase2().residualAlpha()) + ); + + const volScalarField muAlphaI + ( + alpha1*rho1*nu1*alpha2*rho2*nu2 + /(limitedAlpha1*rho1*nu1 + limitedAlpha2*rho2*nu2) + ); + + 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 feb63aea7f..8cd582eb71 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 + >, + sensibleEnthalpy + > +> constRefrPolFluidHThermoPhysics; + typedef constTransport < @@ -110,6 +124,19 @@ constTransport > > constRefFluidEThermoPhysics; +typedef +constTransport +< + species::thermo + < + eRefConstThermo + < + rPolynomial + >, + 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 703c03aee8..01067cdab4 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 0000000000..86f125098d --- /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 . + +\*---------------------------------------------------------------------------*/ + +#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& 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::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 0000000000..4b5e10365b --- /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 . + +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> 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 clone() const + { + return autoPtr + ( + 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& 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 45b1377320..8f6e20d832 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& fx + Field& 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 da313611ea..632c2abb37 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& fx + Field& 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 b21f80110f..3bfbdfee04 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& fx + Field& 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 9859cb51ce..a448157661 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& fx + Field& 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 e53d1e8a76..148281aa81 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& fx + Field& 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 396efc567e..a206b4d6be 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& fx + Field& 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 1a58af51f9..297dde3d65 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& fx + Field& 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 e3f6f3c758..7533e46b8a 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& fx + Field& 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 80c4f03733..87ffb662ea 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& fx + Field& 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 26a4b360b1..f13b976686 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& fx + Field& fx, + const rigidBodyModelState& state ) const { diff --git a/src/rigidBodyDynamics/restraints/softWall/softWall.H b/src/rigidBodyDynamics/restraints/softWall/softWall.H index 7f58b3132e..86e2067982 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& fx + Field& 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 9ea6d8bd38..be8b1fd984 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& fx + Field& 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 c05dac59f9..a80b479786 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& fx + Field& 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 95d854048e..c0dd3988c2 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& fx + Field& 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 77667385ba..046b5eda74 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& fx) const; + void applyRestraints + ( + scalarField& tau, + Field& 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 a7542f6419..1a0f0f8f57 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 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 ad0e563361..ff7a3c26dd 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 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 16ca409657..2b2a6c0de4 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 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 54ea130107..40e640c05e 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 570cf098df..8244c7095a 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 1f7a5006b6..ed2a2476a5 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 103fbf6bfa..c8abe7eef8 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 cbdea7f760..7866817302 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 586a9cfb19..bf84e88d0b 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 44c128db35..7789ebbe91 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 9c442b53eb..c5d8576933 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::W } +template +Foam::scalar Foam::SpecieMixture::Hc(const label speciei) const +{ + return this->getLocalThermo(speciei).Hc(); +} + + template Foam::scalar Foam::SpecieMixture::Cp ( @@ -84,6 +91,18 @@ Foam::scalar Foam::SpecieMixture::Cv } +template +Foam::scalar Foam::SpecieMixture::HE +( + const label speciei, + const scalar p, + const scalar T +) const +{ + return this->getLocalThermo(speciei).HE(p, T); +} + + template Foam::scalar Foam::SpecieMixture::Ha ( @@ -108,16 +127,6 @@ Foam::scalar Foam::SpecieMixture::Hs } -template -Foam::scalar Foam::SpecieMixture::Hc -( - const label speciei -) const -{ - return this->getLocalThermo(speciei).Hc(); -} - - template Foam::scalar Foam::SpecieMixture::S ( diff --git a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H index 3fb34f531c..572dff91e2 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 76194cc4ef..d7db29cc42 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 734daa6cd5..1f09cb4272 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 0000000000..a67c6d7ac1 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#include "rPolynomial.H" +#include "IOstreams.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::rPolynomial::rPolynomial(const dictionary& dict) +: + Specie(dict), + C_(dict.subDict("equationOfState").lookup("C")) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::rPolynomial::write(Ostream& os) const +{ + Specie::write(os); + + dictionary dict("equationOfState"); + dict.add("C", C_); + + os << indent << dict.dictName() << dict; +} + + +// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // + +template +Foam::Ostream& Foam::operator<<(Ostream& os, const rPolynomial& 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 0000000000..44ad5aa517 --- /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 . + +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 rPolynomial; + +template +inline rPolynomial operator+ +( + const rPolynomial&, + const rPolynomial& +); + +template +inline rPolynomial operator* +( + const scalar, + const rPolynomial& +); + +template +inline rPolynomial operator== +( + const rPolynomial&, + const rPolynomial& +); + +template +Ostream& operator<< +( + Ostream&, + const rPolynomial& +); + + +/*---------------------------------------------------------------------------*\ + Class rPolynomial Declaration +\*---------------------------------------------------------------------------*/ + +template +class rPolynomial +: + public Specie +{ + // Private Data + + class coeffList + : + public VectorSpace + { + public: + + // Constructors + + //- Construct null + inline coeffList() + {} + + //- Construct from Istream + inline coeffList(Istream& is) + : + VectorSpace(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 clone() const; + + // Selector from dictionary + inline static autoPtr 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+ + ( + const rPolynomial&, + const rPolynomial& + ); + + friend rPolynomial operator* + ( + const scalar s, + const rPolynomial& + ); + + friend rPolynomial operator== + ( + const rPolynomial&, + const rPolynomial& + ); + + + // Ostream Operator + + friend Ostream& operator<< + ( + 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 0000000000..0f47236307 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#include "rPolynomial.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template +inline Foam::rPolynomial::rPolynomial +( + const Specie& sp, + const coeffList& coeffs +) +: + Specie(sp), + C_(coeffs) +{} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::rPolynomial::rPolynomial +( + const word& name, + const rPolynomial& rp +) +: + Specie(name, rp), + C_(rp.C_) +{} + + +template +inline Foam::autoPtr> +Foam::rPolynomial::clone() const +{ + return autoPtr>::New(*this); +} + + +template +inline Foam::autoPtr> +Foam::rPolynomial::New +( + const dictionary& dict +) +{ + return autoPtr>::New(dict); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +inline Foam::scalar Foam::rPolynomial::rho(scalar p, scalar T) const +{ + return 1/(C_[0] + (C_[1] + C_[2]*T - C_[4]*p)*T - C_[3]*p); +} + + +template +inline Foam::scalar Foam::rPolynomial::H(scalar p, scalar T) const +{ + return 0; +} + + +template +inline Foam::scalar Foam::rPolynomial::Cp(scalar p, scalar T) const +{ + return 0; +} + + +template +inline Foam::scalar Foam::rPolynomial::E(scalar p, scalar T) const +{ + return 0; +} + + +template +inline Foam::scalar Foam::rPolynomial::Cv(scalar p, scalar T) const +{ + return 0; +} + + +template +inline Foam::scalar Foam::rPolynomial::S(scalar p, scalar T) const +{ + return 0; +} + + +template +inline Foam::scalar Foam::rPolynomial::psi(scalar p, scalar T) const +{ + return sqr(rho(p, T))*(C_[3] + C_[4]*T); +} + + +template +inline Foam::scalar Foam::rPolynomial::Z(scalar p, scalar T) const +{ + return 1; +} + + +template +inline Foam::scalar Foam::rPolynomial::CpMCv(scalar p, scalar T) const +{ + return 0; +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template +inline void Foam::rPolynomial::operator+= +( + const rPolynomial& rp +) +{ + const scalar Y1 = this->Y(); + Specie::operator+=(rp); + + if (mag(this->Y()) > SMALL) + { + C_ = (Y1*C_ + rp.Y()*rp.C_)/this->Y(); + } +} + + +template +inline void Foam::rPolynomial::operator*=(const scalar s) +{ + Specie::operator*=(s); +} + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +template +inline Foam::rPolynomial Foam::operator+ +( + const rPolynomial& rp1, + const rPolynomial& rp2 +) +{ + Specie sp + ( + static_cast(rp1) + + static_cast(rp2) + ); + + if (mag(sp.Y()) < SMALL) + { + return rPolynomial + ( + sp, + rp1.C_ + ); + } + + return rPolynomial + ( + sp, + (rp1.Y()*rp1.C_ + rp2.Y()*rp2.C_)/sp.Y() + ); +} + + +template +inline Foam::rPolynomial Foam::operator* +( + const scalar s, + const rPolynomial& rp +) +{ + return rPolynomial + ( + s*static_cast(rp), + rp.C_ + ); +} + + +template +inline Foam::rPolynomial Foam::operator== +( + const rPolynomial& rp1, + const rPolynomial& rp2 +) +{ + Specie sp + ( + static_cast(rp1) + == static_cast(rp2) + ); + + return rPolynomial + ( + 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 2932ad9efb..de1b4e419c 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 + >, + sensibleEnthalpy + > + > + constrPolFluidHThermoPhysics; + typedef constTransport < @@ -280,6 +295,20 @@ namespace Foam > constFluidEThermoPhysics; + typedef + constTransport + < + species::thermo + < + eConstThermo + < + perfectFluid + >, + sensibleInternalEnergy + > + > + constrPolFluidEThermoPhysics; + typedef constTransport <