From a62529c0e9bd21e7d627db608eb344d2c62819a9 Mon Sep 17 00:00:00 2001
From: Andrew Heather <>
Date: Tue, 16 Mar 2021 21:35:20 +0000
Subject: [PATCH] ENH: new SRFrhoSimpleFoam solver
---
.../rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H | 26 +++++
.../rhoSimpleFoam/SRFrhoSimpleFoam/Make/files | 3 +
.../SRFrhoSimpleFoam/Make/options | 21 ++++
.../SRFrhoSimpleFoam/SRFrhoSimpleFoam.C | 99 ++++++++++++++++
.../rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H | 22 ++++
.../SRFrhoSimpleFoam/createFieldRefs.H | 1 +
.../SRFrhoSimpleFoam/createFields.H | 89 ++++++++++++++
.../rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H | 109 ++++++++++++++++++
8 files changed, 370 insertions(+)
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H
create mode 100644 applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H
new file mode 100644
index 00000000000..56e78ff6308
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H
@@ -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();
+}
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files
new file mode 100644
index 00000000000..8ad8dee5650
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files
@@ -0,0 +1,3 @@
+SRFrhoSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/SRFrhoSimpleFoam
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options
new file mode 100644
index 00000000000..806154972b4
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options
@@ -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
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C
new file mode 100644
index 00000000000..fec9ebe4fc7
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H
new file mode 100644
index 00000000000..c1ec71e9469
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H
@@ -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);
+ }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H
new file mode 100644
index 00000000000..502b3b42300
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H
@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H
new file mode 100644
index 00000000000..e333fc546ef
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H
@@ -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"
diff --git a/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H
new file mode 100644
index 00000000000..cfcce1d9419
--- /dev/null
+++ b/applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H
@@ -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();
+}
--
GitLab