Compare commits

...

2 Commits

14 changed files with 530 additions and 101 deletions

View File

@ -0,0 +1,26 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::div(phi, he)
+ (
he.name() == "e"
? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
: fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
)
- fvm::laplacian(turbulence->alphaEff(), he)
==
fvOptions(rho, he)
);
EEqn.relax();
fvOptions.constrain(EEqn);
EEqn.solve();
fvOptions.correct(he);
thermo.correct();
}

View File

@ -0,0 +1,3 @@
SRFrhoSimpleFoam.C
EXE = $(FOAM_APPBIN)/SRFrhoSimpleFoam

View File

@ -0,0 +1,21 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
EXE_LIBS = \
-lfiniteVolume \
-lfvOptions \
-lmeshTools \
-lsampling \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-latmosphericModels

View File

@ -0,0 +1,99 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 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/>.
Application
SRFSimpleFoam
Group
grpIncompressibleSolvers
Description
Steady-state solver for incompressible, turbulent flow of non-Newtonian
fluids in a single rotating frame.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "simpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "SRFModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"Steady-state solver for compressible, turbulent flow"
" of non-Newtonian fluids in a single rotating frame."
);
#include "postProcess.H"
#include "addCheckCaseOptions.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "initContinuityErrs.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity SIMPLE corrector
{
#include "UrelEqn.H"
#include "EEqn.H"
#include "pEqn.H"
}
U = Urel + SRF->U();
turbulence->correct();
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -0,0 +1,22 @@
// Relative momentum predictor
tmp<fvVectorMatrix> tUrelEqn
(
fvm::div(phi, Urel)
+ turbulence->divDevRhoReff(Urel)
+ rho*SRF->Su()
==
fvOptions(rho, Urel)
);
fvVectorMatrix& UrelEqn = tUrelEqn.ref();
UrelEqn.relax();
fvOptions.constrain(UrelEqn);
if (simple.momentumPredictor())
{
solve(UrelEqn == -fvc::grad(p));
fvOptions.correct(Urel);
}

View File

@ -0,0 +1 @@
const volScalarField& psi = thermo.psi();

View File

@ -0,0 +1,89 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidThermo> pThermo
(
fluidThermo::New(mesh)
);
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField& p = thermo.p();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
thermo.rho()
);
Info<< "Reading field Urel\n" << endl;
volVectorField Urel
(
IOobject
(
"Urel",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rho*Urel) & mesh.Sf()
);
pressureControl pressureControl(p, rho, simple.dict());
mesh.setFluxRequired(p.name());
Info<< "Creating SRF model\n" << endl;
autoPtr<SRF::SRFModel> SRF(SRF::SRFModel::New(Urel));
// Construct the absolute velocity
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Urel + SRF->U()
);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
#include "createFvOptions.H"

View File

@ -0,0 +1,109 @@
volScalarField rAUrel(1.0/UrelEqn.A());
surfaceScalarField rhorAUrelf("rhorAUf", fvc::interpolate(rho*rAUrel));
volVectorField HbyA(constrainHbyA(rAUrel*UrelEqn.H(), Urel, p));
tUrelEqn.clear();
bool closedVolume = false;
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, Urel, phiHbyA, rhorAUrelf);
if (simple.transonic())
{
surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUrelf, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
closedVolume = adjustPhi(phiHbyA, Urel, p);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAUrelf, p)
==
fvOptions(psi, p, rho.name())
);
pEqn.setReference
(
pressureControl.refCell(),
pressureControl.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
Urel = HbyA - rAUrel*fvc::grad(p);
Urel.correctBoundaryConditions();
fvOptions.correct(Urel);
bool pLimited = pressureControl.limit(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
if (pLimited || closedVolume)
{
p.correctBoundaryConditions();
}
rho = thermo.rho();
if (!simple.transonic())
{
rho.relax();
}

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -63,14 +64,7 @@ Foam::SRF::SRFModel::SRFModel
mesh_(Urel_.mesh()),
origin_("origin", dimLength, get<vector>("origin")),
axis_(normalised(get<vector>("axis"))),
SRFModelCoeffs_(optionalSubDict(type + "Coeffs")),
omega_(dimensionedVector("omega", dimless/dimTime, Zero))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::SRF::SRFModel::~SRFModel()
SRFModelCoeffs_(optionalSubDict(type + "Coeffs"))
{}
@ -109,29 +103,20 @@ const Foam::vector& Foam::SRF::SRFModel::axis() const
}
const Foam::dimensionedVector& Foam::SRF::SRFModel::omega() const
{
return omega_;
}
Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh>>
Foam::SRF::SRFModel::Fcoriolis() const
{
return tmp<volVectorField::Internal>
return tmp<DimensionedField<vector, Foam::volMesh>>::New
(
new volVectorField::Internal
IOobject
(
IOobject
(
"Fcoriolis",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
2.0*omega_ ^ Urel_
)
"Fcoriolis",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
2.0*omega() ^ Urel_
);
}
@ -139,20 +124,19 @@ Foam::SRF::SRFModel::Fcoriolis() const
Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh>>
Foam::SRF::SRFModel::Fcentrifugal() const
{
return tmp<volVectorField::Internal>
const dimensionedVector omega = this->omega();
return tmp<DimensionedField<vector, Foam::volMesh>>::New
(
new volVectorField::Internal
IOobject
(
IOobject
(
"Fcentrifugal",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
omega_ ^ (omega_ ^ (mesh_.C() - origin_))
)
"Fcentrifugal",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
omega ^ (omega ^ (mesh_.C() - origin_))
);
}
@ -170,7 +154,7 @@ Foam::vectorField Foam::SRF::SRFModel::velocity
) const
{
tmp<vectorField> tfld =
omega_.value()
omega().value()
^ (
(positions - origin_.value())
- axis_*(axis_ & (positions - origin_.value()))
@ -182,21 +166,18 @@ Foam::vectorField Foam::SRF::SRFModel::velocity
Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::U() const
{
return tmp<volVectorField>
return tmp<volVectorField>::New
(
new volVectorField
IOobject
(
IOobject
(
"Usrf",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
omega_
^ ((mesh_.C() - origin_) - axis_*(axis_ & (mesh_.C() - origin_)))
)
"Usrf",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
omega()
^ ((mesh_.C() - origin_) - axis_*(axis_ & (mesh_.C() - origin_)))
);
}
@ -205,21 +186,18 @@ Foam::tmp<Foam::volVectorField> Foam::SRF::SRFModel::Uabs() const
{
tmp<volVectorField> Usrf = U();
tmp<volVectorField> tUabs
auto tUabs = tmp<volVectorField>::New
(
new volVectorField
IOobject
(
IOobject
(
"Uabs",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
Usrf
)
"Uabs",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
Usrf
);
volVectorField& Uabs = tUabs.ref();

View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -50,6 +51,7 @@ SourceFiles
#include "fvMesh.H"
#include "volFields.H"
#include "vectorField.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,9 +88,6 @@ protected:
//- SRF model coefficients dictionary
dictionary SRFModelCoeffs_;
//- Angular velocity of the frame (rad/s)
dimensionedVector omega_;
private:
@ -141,7 +140,7 @@ public:
//- Destructor
virtual ~SRFModel();
virtual ~SRFModel() = default;
// Member Functions
@ -154,15 +153,15 @@ public:
// Access
//- Return the angular velocity field [rad/s]
virtual dimensionedVector omega() const = 0;
//- Return the origin of rotation
const dimensionedVector& origin() const;
//- Return the axis of rotation
const vector& axis() const;
//- Return the angular velocity field [rad/s]
const dimensionedVector& omega() const;
//- Return the coriolis force
tmp<volVectorField::Internal> Fcoriolis() const;

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -47,16 +47,14 @@ namespace Foam
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::SRF::rpm::rpm(const volVectorField& U)
:
SRFModel(typeName, U),
rpm_(SRFModelCoeffs_.get<scalar>("rpm"))
{
// The angular velocity
omega_.value() = axis_*rpmToRads(rpm_);
}
rpmPtr_(Function1<scalar>::New("rpm", SRFModelCoeffs_))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -65,9 +63,7 @@ bool Foam::SRF::rpm::read()
{
if (SRFModel::read())
{
SRFModelCoeffs_.readEntry("rpm", rpm_);
omega_.value() = axis_*rpmToRads(rpm_);
rpmPtr_.reset(Function1<scalar>::New("rpm", SRFModelCoeffs_));
return true;
}
@ -76,4 +72,14 @@ bool Foam::SRF::rpm::read()
}
Foam::dimensionedVector Foam::SRF::rpm::omega() const
{
const scalar t = this->db().time().timeOutputValue();
const scalar rpm = rpmPtr_->value(t);
return dimensionedVector("omega", dimless/dimTime, axis_*rpmToRads(rpm));
}
// ************************************************************************* //

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011 OpenFOAM Foundation
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -40,6 +40,7 @@ SourceFiles
#define SRF_rpm_H
#include "SRFModel.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -60,7 +61,7 @@ class rpm
// Private data
//- Revolutions per minute
scalar rpm_;
autoPtr<Function1<scalar>> rpmPtr_;
// Private Member Functions
@ -85,13 +86,16 @@ public:
//- Destructor
~rpm() = default;
virtual ~rpm() = default;
// Member functions
//- Read coefficients
bool read();
//- Return the angular velocity field [rad/s]
virtual dimensionedVector omega() const;
};

View File

@ -41,7 +41,8 @@ Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
:
fixedValueFvPatchVectorField(p, iF),
relative_(false),
inletValue_(p.size(), Zero)
inletValue_(p.size(), Zero),
refValuePtr_(fvPatchVectorField::New("refValue", p, iF))
{}
@ -55,8 +56,15 @@ Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
:
fixedValueFvPatchVectorField(ptf, p, iF, mapper),
relative_(ptf.relative_),
inletValue_(ptf.inletValue_, mapper)
{}
inletValue_(ptf.inletValue_, mapper),
refValuePtr_()
{
if (ptf.refValuePtr_.valid())
{
refValuePtr_ =
fvPatchVectorField::New(ptf.refValuePtr_(), p, iF, mapper);
}
}
Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
@ -68,8 +76,22 @@ Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
:
fixedValueFvPatchVectorField(p, iF, dict),
relative_(dict.get<Switch>("relative")),
inletValue_("inletValue", dict, p.size())
{}
inletValue_(p.size(), Zero),
refValuePtr_()
{
if (dict.found("inletValue"))
{
inletValue_ = vectorField("inletValue", dict, p.size());
}
else
{
if (dict.found("refValue"))
{
refValuePtr_ =
fvPatchVectorField::New(p, iF, dict.subDict("refValue"));
}
}
}
Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
@ -79,8 +101,14 @@ Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
:
fixedValueFvPatchVectorField(srfvpvf),
relative_(srfvpvf.relative_),
inletValue_(srfvpvf.inletValue_)
{}
inletValue_(srfvpvf.inletValue_),
refValuePtr_()
{
if (srfvpvf.refValuePtr_.valid())
{
refValuePtr_ = srfvpvf.refValuePtr_().clone();
}
}
Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
@ -91,8 +119,14 @@ Foam::SRFVelocityFvPatchVectorField::SRFVelocityFvPatchVectorField
:
fixedValueFvPatchVectorField(srfvpvf, iF),
relative_(srfvpvf.relative_),
inletValue_(srfvpvf.inletValue_)
{}
inletValue_(srfvpvf.inletValue_),
refValuePtr_()
{
if (srfvpvf.refValuePtr_.valid())
{
refValuePtr_ = srfvpvf.refValuePtr_().clone();
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -104,6 +138,11 @@ void Foam::SRFVelocityFvPatchVectorField::autoMap
{
vectorField::autoMap(m);
inletValue_.autoMap(m);
if (refValuePtr_.valid())
{
refValuePtr_->autoMap(m);
}
}
@ -115,10 +154,14 @@ void Foam::SRFVelocityFvPatchVectorField::rmap
{
fixedValueFvPatchVectorField::rmap(ptf, addr);
const SRFVelocityFvPatchVectorField& tiptf =
refCast<const SRFVelocityFvPatchVectorField>(ptf);
const auto& srfptf = refCast<const SRFVelocityFvPatchVectorField>(ptf);
inletValue_.rmap(tiptf.inletValue_, addr);
inletValue_.rmap(srfptf.inletValue_, addr);
if (srfptf.refValuePtr_.valid())
{
refValuePtr_->rmap(srfptf.refValuePtr_(), addr);
}
}
@ -133,19 +176,34 @@ void Foam::SRFVelocityFvPatchVectorField::updateCoeffs()
if (!relative_)
{
// Get reference to the SRF model
const SRF::SRFModel& srf =
db().lookupObject<SRF::SRFModel>("SRFProperties");
const auto& srf = db().lookupObject<SRF::SRFModel>("SRFProperties");
// Determine patch velocity due to SRF
const vectorField SRFVelocity(srf.velocity(patch().Cf()));
operator==(-SRFVelocity + inletValue_);
if (refValuePtr_.valid())
{
refValuePtr_->evaluate();
operator==(-SRFVelocity + refValuePtr_());
}
else
{
operator==(-SRFVelocity + inletValue_);
}
}
// If already relative to the SRF simply supply the inlet value as a fixed
// value
else
{
operator==(inletValue_);
if (refValuePtr_.valid())
{
refValuePtr_->evaluate();
operator==(refValuePtr_());
}
else
{
operator==(inletValue_);
}
}
fixedValueFvPatchVectorField::updateCoeffs();
@ -156,7 +214,18 @@ void Foam::SRFVelocityFvPatchVectorField::write(Ostream& os) const
{
fvPatchVectorField::write(os);
os.writeEntry("relative", relative_);
inletValue_.writeEntry("inletValue", os);
if (refValuePtr_.valid())
{
os.beginBlock("refValue");
refValuePtr_->write(os);
os.endBlock();
}
else
{
inletValue_.writeEntry("inletValue", os);
}
writeEntry("value", os);
}

View File

@ -114,6 +114,9 @@ class SRFVelocityFvPatchVectorField
//- Inlet value [m/s]
vectorField inletValue_;
//- Optional patch to retrieve values from
tmp<fvPatchVectorField> refValuePtr_;
public: