New modular solver framework for single- and multi-region simulations

in which different solver modules can be selected in each region to for complex
conjugate heat-transfer and other combined physics problems such as FSI
(fluid-structure interaction).

For single-region simulations the solver module is selected, instantiated and
executed in the PIMPLE loop in the new foamRun application.

For multi-region simulations the set of solver modules, one for each region, are
selected, instantiated and executed in the multi-region PIMPLE loop of new the
foamMultiRun application.

This provides a very general, flexible and extensible framework for complex
coupled problems by creating more solver modules, either by converting existing
solver applications or creating new ones.

The current set of solver modules provided are:

isothermalFluid
    Solver module for steady or transient turbulent flow of compressible
    isothermal fluids with optional mesh motion and mesh topology changes.

    Created from the rhoSimpleFoam, rhoPimpleFoam and buoyantFoam solvers but
    without the energy equation, hence isothermal.  The buoyant pressure
    formulation corresponding to the buoyantFoam solver is selected
    automatically by the presence of the p_rgh pressure field in the start-time
    directory.

fluid
    Solver module for steady or transient turbulent flow of compressible fluids
    with heat-transfer for HVAC and similar applications, with optional
    mesh motion and mesh topology changes.

    Derived from the isothermalFluid solver module with the addition of the
    energy equation from the rhoSimpleFoam, rhoPimpleFoam and buoyantFoam
    solvers, thus providing the equivalent functionality of these three solvers.

multicomponentFluid
    Solver module for steady or transient turbulent flow of compressible
    reacting fluids with optional mesh motion and mesh topology changes.

    Derived from the isothermalFluid solver module with the addition of
    multicomponent thermophysical properties energy and specie mass-fraction
    equations from the reactingFoam solver, thus providing the equivalent
    functionality in reactingFoam and buoyantReactingFoam.  Chemical reactions
    and/or combustion modelling may be optionally selected to simulate reacting
    systems including fires, explosions etc.

solid
    Solver module for turbulent flow of compressible fluids for conjugate heat
    transfer, HVAC and similar applications, with optional mesh motion and mesh
    topology changes.

    The solid solver module may be selected in solid regions of a CHT case, with
    either the fluid or multicomponentFluid solver module in the fluid regions
    and executed with foamMultiRun to provide functionality equivalent
    chtMultiRegionFoam but in a flexible and extensible framework for future
    extension to more complex coupled problems.

All the usual fvModels, fvConstraints, functionObjects etc. are available with
these solver modules to support simulations including body-forces, local sources,
Lagrangian clouds, liquid films etc. etc.

Converting compressibleInterFoam and multiphaseEulerFoam into solver modules
would provide a significant enhancement to the CHT capability and incompressible
solvers like pimpleFoam run in conjunction with solidDisplacementFoam in
foamMultiRun would be useful for a range of FSI problems.  Many other
combinations of existing solvers converted into solver modules could prove
useful for a very wide range of complex combined physics simulations.

All tutorials from the rhoSimpleFoam, rhoPimpleFoam, buoyantFoam, reactingFoam,
buoyantReactingFoam and chtMultiRegionFoam solver applications replaced by
solver modules have been updated and moved into the tutorials/modules directory:

modules
├── CHT
│   ├── coolingCylinder2D
│   ├── coolingSphere
│   ├── heatedDuct
│   ├── heatExchanger
│   ├── reverseBurner
│   └── shellAndTubeHeatExchanger
├── fluid
│   ├── aerofoilNACA0012
│   ├── aerofoilNACA0012Steady
│   ├── angledDuct
│   ├── angledDuctExplicitFixedCoeff
│   ├── angledDuctLTS
│   ├── annularThermalMixer
│   ├── BernardCells
│   ├── blockedChannel
│   ├── buoyantCavity
│   ├── cavity
│   ├── circuitBoardCooling
│   ├── decompressionTank
│   ├── externalCoupledCavity
│   ├── forwardStep
│   ├── helmholtzResonance
│   ├── hotRadiationRoom
│   ├── hotRadiationRoomFvDOM
│   ├── hotRoom
│   ├── hotRoomBoussinesq
│   ├── hotRoomBoussinesqSteady
│   ├── hotRoomComfort
│   ├── iglooWithFridges
│   ├── mixerVessel2DMRF
│   ├── nacaAirfoil
│   ├── pitzDaily
│   ├── prism
│   ├── shockTube
│   ├── squareBend
│   ├── squareBendLiq
│   └── squareBendLiqSteady
└── multicomponentFluid
    ├── aachenBomb
    ├── counterFlowFlame2D
    ├── counterFlowFlame2D_GRI
    ├── counterFlowFlame2D_GRI_TDAC
    ├── counterFlowFlame2DLTS
    ├── counterFlowFlame2DLTS_GRI_TDAC
    ├── cylinder
    ├── DLR_A_LTS
    ├── filter
    ├── hotBoxes
    ├── membrane
    ├── parcelInBox
    ├── rivuletPanel
    ├── SandiaD_LTS
    ├── simplifiedSiwek
    ├── smallPoolFire2D
    ├── smallPoolFire3D
    ├── splashPanel
    ├── verticalChannel
    ├── verticalChannelLTS
    └── verticalChannelSteady

Also redirection scripts are provided for the replaced solvers which call
foamRun -solver <solver module name> or foamMultiRun in the case of
chtMultiRegionFoam for backward-compatibility.

Documentation for foamRun and foamMultiRun:

Application
    foamRun

Description
    Loads and executes an OpenFOAM solver module either specified by the
    optional \c solver entry in the \c controlDict or as a command-line
    argument.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

Usage
    \b foamRun [OPTION]

      - \par -solver <name>
        Solver name

      - \par -libs '(\"lib1.so\" ... \"libN.so\")'
        Specify the additional libraries loaded

    Example usage:
      - To run a \c rhoPimpleFoam case by specifying the solver on the
        command line:
        \verbatim
            foamRun -solver fluid
        \endverbatim

      - To update and run a \c rhoPimpleFoam case add the following entries to
        the controlDict:
        \verbatim
            application     foamRun;

            solver          fluid;
        \endverbatim
        then execute \c foamRun

Application
    foamMultiRun

Description
    Loads and executes an OpenFOAM solver modules for each region of a
    multiregion simulation e.g. for conjugate heat transfer.

    The region solvers are specified in the \c regionSolvers dictionary entry in
    \c controlDict, containing a list of pairs of region and solver names,
    e.g. for a two region case with one fluid region named
    liquid and one solid region named tubeWall:
    \verbatim
        regionSolvers
        {
            liquid          fluid;
            tubeWall        solid;
        }
    \endverbatim

    The \c regionSolvers entry is a dictionary to support name substitutions to
    simplify the specification of a single solver type for a set of
    regions, e.g.
    \verbatim
        fluidSolver     fluid;
        solidSolver     solid;

        regionSolvers
        {
            tube1             $fluidSolver;
            tubeWall1         solid;
            tube2             $fluidSolver;
            tubeWall2         solid;
            tube3             $fluidSolver;
            tubeWall3         solid;
        }
    \endverbatim

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

Usage
    \b foamMultiRun [OPTION]

      - \par -libs '(\"lib1.so\" ... \"libN.so\")'
        Specify the additional libraries loaded

    Example usage:
      - To update and run a \c chtMultiRegion case add the following entries to
        the controlDict:
        \verbatim
            application     foamMultiRun;

            regionSolvers
            {
                fluid           fluid;
                solid           solid;
            }
        \endverbatim
        then execute \c foamMultiRun
This commit is contained in:
Henry Weller
2022-08-04 21:11:35 +01:00
parent 7e2dd6dda2
commit 968e60148a
1540 changed files with 5034 additions and 5070 deletions

View File

@ -1 +1,19 @@
#include "../../compressible/rhoPimpleFoam/correctPhi.H"
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
correctUphiBCs(rho, U, phi, true);
CorrectPhi
(
phi,
p,
rho,
psi,
dimensionedScalar("rAUf", dimTime, 1),
divrhoU(),
pimple
);
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);

View File

@ -1,31 +0,0 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
: -dpdt
)
+ thermophysicalTransport->divq(he)
==
reaction->Qdot()
+ fvModels.source(rho, he)
);
EEqn.relax();
fvConstraints.constrain(EEqn);
EEqn.solve();
fvConstraints.constrain(he);
thermo.correct();
Info<< "min/max(T) = "
<< min(T).value() << ", " << max(T).value() << endl;
}

View File

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

View File

@ -1,33 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/fluidMulticomponentThermo/lnInclude \
-I$(LIB_SRC)/physicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lchemistryModel \
-lODE \
-lcombustionModels \
-lmulticomponentThermophysicalModels \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \
-lfluidMulticomponentThermophysicalTransportModels \
-lfiniteVolume \
-lfvModels \
-lfvConstraints \
-lmeshTools \
-lsampling

View File

@ -1,25 +0,0 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels.source(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvConstraints.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
}

View File

@ -1,40 +0,0 @@
tmp<fv::convectionScheme<scalar>> mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.schemes().div("div(phi,Yi_h)")
)
);
reaction->correct();
forAll(Y, i)
{
if (composition.solve(i))
{
volScalarField& Yi = Y[i];
fvScalarMatrix YiEqn
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
+ thermophysicalTransport->divj(Yi)
==
reaction->R(Yi)
+ fvModels.source(rho, Yi)
);
YiEqn.relax();
fvConstraints.constrain(YiEqn);
YiEqn.solve("Yi");
fvConstraints.constrain(Yi);
}
}
composition.normalise();

View File

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

View File

@ -1,36 +0,0 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels.source(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvConstraints.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
}

View File

@ -1,211 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
buoyantReactingFoam
Description
Transient solver for turbulent flow of compressible reacting fluids with
enhanced buoyancy treatment and optional mesh motion and mesh topology
changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidMulticomponentThermo.H"
#include "combustionModel.H"
#include "compressibleMomentumTransportModels.H"
#include "fluidMulticomponentThermophysicalTransportModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
#include "pressureReference.H"
#include "hydrostaticInitialisation.H"
#include "CorrectPhi.H"
#include "fvModels.H"
#include "fvConstraints.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
turbulence->validate();
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
divrhoU = new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
}
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
fvModels.preUpdateMesh();
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU = new volVectorField("rhoU", rho*U);
}
// Update the mesh for topology change, mesh to mesh mapping
mesh.update();
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (!pimple.flow())
{
if (pimple.models())
{
fvModels.correct();
}
if (pimple.thermophysics())
{
#include "YEqn.H"
#include "EEqn.H"
}
}
else
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
// Move the mesh
mesh.move();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
MRF.update();
if (correctPhi)
{
#include "correctPhi.H"
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
if (pimple.firstPimpleIter() && !pimple.simpleRho())
{
#include "rhoEqn.H"
}
if (pimple.models())
{
fvModels.correct();
}
#include "UEqn.H"
if (pimple.thermophysics())
{
#include "YEqn.H"
#include "EEqn.H"
}
// --- Pressure corrector loop
while (pimple.correct())
{
#include "../../../heatTransfer/buoyantFoam/pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
thermophysicalTransport->correct();
}
}
}
rho = thermo.rho();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

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

View File

@ -1,143 +0,0 @@
#include "createRDeltaT.H"
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidMulticomponentThermo> pThermo
(
fluidMulticomponentThermo::New(mesh)
);
fluidMulticomponentThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh
),
thermo.rho()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField& p = thermo.p();
#include "compressibleCreatePhi.H"
mesh.schemes().setFluxRequired(p.name());
Info << "Creating turbulence model.\n" << nl;
autoPtr<compressible::momentumTransportModel> turbulence
(
compressible::momentumTransportModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Creating thermophysical transport model\n" << endl;
autoPtr<fluidMulticomponentThermophysicalTransportModel>
thermophysicalTransport
(
fluidMulticomponentThermophysicalTransportModel::New
(
turbulence(),
thermo
)
);
Info<< "Creating reaction model\n" << endl;
autoPtr<combustionModel> reaction(combustionModel::New(thermo, turbulence()));
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
#include "readpRef.H"
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
pressureReference pressureReference
(
p,
p_rgh,
pimple.dict(),
thermo.incompressible()
);
mesh.schemes().setFluxRequired(p_rgh.name());
hydrostaticInitialisation
(
p_rgh,
p,
rho,
U,
gh,
ghf,
pRef,
thermo,
pimple.dict()
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(thermo.he());
#include "createMRF.H"
#include "createFvModels.H"
#include "createFvConstraints.H"

View File

@ -1 +0,0 @@
#include "../../compressible/rhoPimpleFoam/correctPhi.H"

View File

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

View File

@ -1,105 +0,0 @@
#include "createRDeltaT.H"
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidMulticomponentThermo> pThermo
(
fluidMulticomponentThermo::New(mesh)
);
fluidMulticomponentThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh
),
thermo.rho()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
volScalarField& p = thermo.p();
#include "compressibleCreatePhi.H"
pressureReference pressureReference(p, pimple.dict(), false);
mesh.schemes().setFluxRequired(p.name());
Info << "Creating turbulence model.\n" << nl;
autoPtr<compressible::momentumTransportModel> turbulence
(
compressible::momentumTransportModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Creating thermophysical transport model\n" << endl;
autoPtr<fluidMulticomponentThermophysicalTransportModel>
thermophysicalTransport
(
fluidMulticomponentThermophysicalTransportModel::New
(
turbulence(),
thermo
)
);
Info<< "Creating reaction model\n" << endl;
autoPtr<combustionModel> reaction(combustionModel::New(thermo, turbulence()));
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(thermo.he());
#include "createMRF.H"
#include "createFvModels.H"
#include "createFvConstraints.H"

View File

@ -1,206 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
reactingFoam
Description
Transient solver for turbulent flow of compressible reacting fluids with
optional mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidMulticomponentThermo.H"
#include "combustionModel.H"
#include "compressibleMomentumTransportModels.H"
#include "fluidMulticomponentThermophysicalTransportModel.H"
#include "multivariateScheme.H"
#include "pimpleControl.H"
#include "pressureReference.H"
#include "CorrectPhi.H"
#include "fvModels.H"
#include "fvConstraints.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
turbulence->validate();
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
divrhoU = new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
}
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
fvModels.preUpdateMesh();
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU = new volVectorField("rhoU", rho*U);
}
// Update the mesh for topology change, mesh to mesh mapping
mesh.update();
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (!pimple.flow())
{
if (pimple.models())
{
fvModels.correct();
}
if (pimple.thermophysics())
{
#include "YEqn.H"
#include "EEqn.H"
}
}
else
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
// Move the mesh
mesh.move();
if (mesh.changing())
{
MRF.update();
if (correctPhi)
{
#include "correctPhi.H"
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
if (pimple.firstPimpleIter() && !pimple.simpleRho())
{
#include "rhoEqn.H"
}
if (pimple.models())
{
fvModels.correct();
}
#include "UEqn.H"
if (pimple.thermophysics())
{
#include "YEqn.H"
#include "EEqn.H"
}
// --- Pressure corrector loop
while (pimple.correct())
{
#include "../../compressible/rhoPimpleFoam/pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
thermophysicalTransport->correct();
}
}
}
rho = thermo.rho();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,27 +0,0 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + fvm::div(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? fvc::div(fvc::absolute(phi, rho, U), p/rho)
: -dpdt
)
+ thermophysicalTransport->divq(he)
==
fvModels.source(rho, he)
);
EEqn.relax();
fvConstraints.constrain(EEqn);
EEqn.solve();
fvConstraints.constrain(he);
thermo.correct();
}

View File

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

View File

@ -1,25 +0,0 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels.source(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvConstraints.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
}

View File

@ -1,19 +0,0 @@
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
correctUphiBCs(rho, U, phi, true);
CorrectPhi
(
phi,
p,
rho,
psi,
dimensionedScalar("rAUf", dimTime, 1),
divrhoU(),
pimple
);
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);

View File

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

View File

@ -1,90 +0,0 @@
#include "createRDeltaT.H"
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 U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
pressureReference pressureReference
(
p,
pimple.dict(),
thermo.incompressible()
);
mesh.schemes().setFluxRequired(p.name());
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::momentumTransportModel> turbulence
(
compressible::momentumTransportModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Creating thermophysical transport model\n" << endl;
autoPtr<fluidThermophysicalTransportModel> thermophysicalTransport
(
fluidThermophysicalTransportModel::New(turbulence(), thermo)
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
#include "createMRF.H"
#include "createFvModels.H"
#include "createFvConstraints.H"

View File

@ -1,200 +0,0 @@
rho = thermo.rho();
rho.relax();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
const volScalarField rAU("rAU", 1.0/UEqn.A());
const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
tmp<volScalarField> rAtU
(
pimple.consistent()
? volScalarField::New("rAtU", 1.0/(1.0/rAU - UEqn.H1()))
: tmp<volScalarField>(nullptr)
);
tmp<surfaceScalarField> rhorAtUf
(
pimple.consistent()
? surfaceScalarField::New("rhoRAtUf", fvc::interpolate(rho*rAtU()))
: tmp<surfaceScalarField>(nullptr)
);
const volScalarField& rAAtU = pimple.consistent() ? rAtU() : rAU;
const surfaceScalarField& rhorAAtUf =
pimple.consistent() ? rhorAtUf() : rhorAUf;
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
if (pimple.nCorrPiso() <= 1)
{
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
bool adjustMass = false;
if (pimple.transonic())
{
const surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
const fvScalarMatrix divPhidp(fvm::div(phid, p));
phiHbyA -= divPhidp.flux();
if (pimple.consistent())
{
phiHbyA += (rhorAAtUf - rhorAUf)*fvc::snGrad(p)*mesh.magSf();
HbyA += (rAAtU - rAU)*fvc::grad(p);
}
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA) + divPhidp
==
fvModels.source(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAAtUf, p));
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
if (mesh.schemes().steady())
{
adjustMass = adjustPhi(phiHbyA, U, p);
}
if (pimple.consistent())
{
phiHbyA += (rhorAAtUf - rhorAUf)*fvc::snGrad(p)*mesh.magSf();
HbyA += (rAAtU - rAU)*fvc::grad(p);
}
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
fvModels.source(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAAtUf, p));
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
if (mesh.schemes().steady())
{
#include "incompressible/continuityErrs.H"
}
else
{
const bool constrained = fvConstraints.constrain(p);
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (constrained)
{
rho = thermo.rho();
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
}
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAAtU*fvc::grad(p);
U.correctBoundaryConditions();
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
if (mesh.schemes().steady())
{
fvConstraints.constrain(p);
}
// For steady compressible closed-volume cases adjust the pressure level
// to obey overall mass continuity
if (adjustMass && !thermo.incompressible())
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/fvc::domainIntegrate(psi);
p.correctBoundaryConditions();
}
if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass)
{
rho = thermo.rho();
}
// Correct rhoUf if the mesh is moving
fvc::correctRhoUf(rhoUf, rho, U, phi, MRF);
if (mesh.schemes().steady() || pimple.simpleRho())
{
rho.relax();
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}

View File

@ -1,188 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
rhoPimpleFoam
Description
Transient solver for turbulent flow of compressible fluids for HVAC and
similar applications, with optional mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidThermo.H"
#include "compressibleMomentumTransportModels.H"
#include "fluidThermophysicalTransportModel.H"
#include "pimpleControl.H"
#include "pressureReference.H"
#include "CorrectPhi.H"
#include "fvModels.H"
#include "fvConstraints.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
turbulence->validate();
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
divrhoU = new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
}
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
fvModels.preUpdateMesh();
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU = new volVectorField("rhoU", rho*U);
}
// Update the mesh for topology change, mesh to mesh mapping
mesh.update();
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
// Move the mesh
mesh.move();
if (mesh.changing())
{
MRF.update();
if (correctPhi)
{
#include "correctPhi.H"
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
if
(
!mesh.schemes().steady()
&& !pimple.simpleRho()
&& pimple.firstPimpleIter()
)
{
#include "rhoEqn.H"
}
fvModels.correct();
#include "UEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
thermophysicalTransport->correct();
}
}
if (!mesh.schemes().steady())
{
rho = thermo.rho();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

View File

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

View File

@ -1,21 +0,0 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels.source(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvConstraints.constrain(UEqn);
solve(UEqn == -fvc::grad(p));
fvConstraints.constrain(U);

View File

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

View File

@ -1,141 +0,0 @@
rho = thermo.rho();
rho.relax();
const volScalarField rAU("rAU", 1.0/UEqn.A());
const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
tmp<volScalarField> rAtU
(
simple.consistent()
? volScalarField::New("rAtU", 1.0/(1.0/rAU - UEqn.H1()))
: tmp<volScalarField>(nullptr)
);
tmp<surfaceScalarField> rhorAtUf
(
simple.consistent()
? surfaceScalarField::New("rhoRAtUf", fvc::interpolate(rho*rAtU()))
: tmp<surfaceScalarField>(nullptr)
);
const volScalarField& rAAtU = simple.consistent() ? rAtU() : rAU;
const surfaceScalarField& rhorAAtUf =
simple.consistent() ? rhorAtUf() : rhorAUf;
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
tUEqn.clear();
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
bool closedVolume = false;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
if (simple.transonic())
{
const surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*phiHbyA/fvc::interpolate(rho)
);
const fvScalarMatrix divPhidp(fvm::div(phid, p));
phiHbyA -= divPhidp.flux();
if (simple.consistent())
{
phiHbyA += (rhorAAtUf - rhorAUf)*fvc::snGrad(p)*mesh.magSf();
HbyA += (rAAtU - rAU)*fvc::grad(p);
}
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA) + divPhidp
- fvm::laplacian(rhorAAtUf, p)
==
fvModels.source(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
closedVolume = adjustPhi(phiHbyA, U, p);
if (simple.consistent())
{
phiHbyA += (rhorAAtUf - rhorAUf)*fvc::snGrad(p)*mesh.magSf();
HbyA += (rAAtU - rAU)*fvc::grad(p);
}
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rhorAAtUf, p)
==
fvModels.source(psi, p, rho.name())
);
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAAtU*fvc::grad(p);
U.correctBoundaryConditions();
fvConstraints.constrain(U);
fvConstraints.constrain(p);
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume && !thermo.incompressible())
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
p.correctBoundaryConditions();
}
rho = thermo.rho();
rho.relax();

View File

@ -1,88 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
rhoSimpleFoam
Description
Steady-state solver for turbulent flow of compressible fluids.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidThermo.H"
#include "compressibleMomentumTransportModels.H"
#include "fluidThermophysicalTransportModel.H"
#include "simpleControl.H"
#include "pressureReference.H"
#include "fvModels.H"
#include "fvConstraints.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.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(runTime))
{
Info<< "Time = " << runTime.userTimeName() << nl << endl;
fvModels.correct();
// Pressure-velocity SIMPLE corrector
#include "UEqn.H"
#include "EEqn.H"
#include "pEqn.H"
turbulence->correct();
thermophysicalTransport->correct();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,28 +0,0 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + fvm::div(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? fvc::div(fvc::absolute(phi, rho, U), p/rho)
: -dpdt
)
+ thermophysicalTransport->divq(he)
==
rho*(U&g)
+ fvModels.source(rho, he)
);
EEqn.relax();
fvConstraints.constrain(EEqn);
EEqn.solve();
fvConstraints.constrain(he);
thermo.correct();
}

View File

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

View File

@ -1,24 +0,0 @@
EXE_INC = \
-I. \
-I../../compressible/rhoPimpleFoam \
-I$(LIB_SRC)/physicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvModels \
-lfvConstraints

View File

@ -1,36 +0,0 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels.source(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvConstraints.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
}

View File

@ -1,207 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
buoyantFoam
Description
Solver for steady or transient buoyant, turbulent flow of compressible
fluids for ventilation and heat-transfer, with optional mesh motion and mesh
topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidThermo.H"
#include "compressibleMomentumTransportModels.H"
#include "fluidThermophysicalTransportModel.H"
#include "pimpleControl.H"
#include "pressureReference.H"
#include "hydrostaticInitialisation.H"
#include "CorrectPhi.H"
#include "fvModels.H"
#include "fvConstraints.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDyMControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
turbulence->validate();
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
divrhoU = new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
}
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
fvModels.preUpdateMesh();
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU = new volVectorField("rhoU", rho*U);
}
// Update the mesh for topology change, mesh to mesh mapping
mesh.update();
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (!pimple.flow())
{
if (pimple.models())
{
fvModels.correct();
}
if (pimple.thermophysics())
{
#include "EEqn.H"
}
}
else
{
if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
{
// Move the mesh
mesh.move();
if (mesh.changing())
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
MRF.update();
if (correctPhi)
{
#include "correctPhi.H"
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
if (pimple.firstPimpleIter() && !pimple.simpleRho())
{
#include "rhoEqn.H"
}
if (pimple.models())
{
fvModels.correct();
}
#include "UEqn.H"
if (pimple.thermophysics())
{
#include "EEqn.H"
}
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
thermophysicalTransport->correct();
}
}
}
rho = thermo.rho();
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,19 +0,0 @@
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
correctUphiBCs(rho, U, phi, true);
CorrectPhi
(
phi,
p_rgh,
rho,
psi,
dimensionedScalar("rAUf", dimTime, 1),
divrhoU(),
pimple
);
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);

View File

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

View File

@ -1,123 +0,0 @@
#include "createRDeltaT.H"
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::NO_READ,
IOobject::NO_WRITE
),
thermo.rho()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::momentumTransportModel> turbulence
(
compressible::momentumTransportModel::New
(
rho,
U,
phi,
thermo
)
);
Info<< "Creating thermophysical transport model\n" << endl;
autoPtr<fluidThermophysicalTransportModel> thermophysicalTransport
(
fluidThermophysicalTransportModel::New(turbulence(), thermo)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
#include "readpRef.H"
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
pressureReference pressureReference
(
p,
p_rgh,
pimple.dict(),
thermo.incompressible()
);
mesh.schemes().setFluxRequired(p_rgh.name());
hydrostaticInitialisation
(
p_rgh,
p,
rho,
U,
gh,
ghf,
pRef,
thermo,
pimple.dict()
);
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(p.dimensions()/dimTime, 0)
);
Info<< "Creating field kinetic energy K\n" << endl;
volScalarField K("K", 0.5*magSqr(U));
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
#include "createMRF.H"
#include "createFvModels.H"
#include "createFvConstraints.H"
#include "checkRadiationModel.H"

View File

@ -1,174 +0,0 @@
rho = thermo.rho();
rho.relax();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
const volScalarField rAU("rAU", 1.0/UEqn.A());
const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
if (pimple.nCorrPiso() <= 1)
{
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
bool adjustMass = mesh.schemes().steady() && adjustPhi(phiHbyA, U, p_rgh);
const surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
phiHbyA += phig;
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
fvScalarMatrix p_rghEqn(p_rgh, dimMass/dimTime);
if (pimple.transonic())
{
const surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
const fvScalarMatrix divPhidp(fvm::div(phid, p));
phiHbyA -= divPhidp.flux();
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA) + divPhidp
==
fvModels.source(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
p_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh);
// Relax the pressure equation to ensure diagonal-dominance
p_rghEqn.relax();
p_rghEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
p_rghEqn.solve();
}
}
else
{
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
==
fvModels.source(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
p_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh);
p_rghEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
p_rghEqn.solve();
}
}
phi = phiHbyA + p_rghEqn.flux();
if (mesh.schemes().steady())
{
#include "incompressible/continuityErrs.H"
}
else
{
p = p_rgh + rho*gh + pRef;
const bool constrained = fvConstraints.constrain(p);
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (constrained)
{
rho = thermo.rho();
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
}
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
p = p_rgh + rho*gh + pRef;
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
U.correctBoundaryConditions();
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
if (mesh.schemes().steady())
{
if (fvConstraints.constrain(p))
{
p_rgh = p - rho*gh - pRef;
p_rgh.correctBoundaryConditions();
}
}
// For steady closed-volume compressible cases adjust the pressure level
// to obey overall mass continuity
if (adjustMass && !thermo.incompressible())
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh - pRef;
p_rgh.correctBoundaryConditions();
}
if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass)
{
rho = thermo.rho();
}
if (mesh.schemes().steady() || pimple.simpleRho())
{
rho.relax();
}
// Correct rhoUf if the mesh is moving
fvc::correctRhoUf(rhoUf, rho, U, phi, MRF);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}

View File

@ -1,5 +0,0 @@
fluid/compressibleCourantNo.C
solid/solidRegionDiffNo.C
chtMultiRegionFoam.C
EXE = $(FOAM_APPBIN)/chtMultiRegionFoam

View File

@ -1,44 +0,0 @@
EXE_INC = \
-I. \
-I./fluid \
-I./solid \
-I./porousFluid \
-I./porousSolid \
-I./include \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/physicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/solidThermo/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/fluidMulticomponentThermo/lnInclude \
-I$(LIB_SRC)/regionModels/regionModel/lnInclude
EXE_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lmulticomponentThermophysicalModels \
-lsolidThermo \
-lchemistryModel \
-lODE \
-lcombustionModels \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \
-lfluidMulticomponentThermophysicalTransportModels \
-lmeshTools \
-lfiniteVolume \
-lfvModels \
-lfvConstraints \
-lregionModels \
-lsampling

View File

@ -1,129 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
chtMultiRegionFoam
Description
Solver for steady or transient fluid flow and solid heat conduction, with
conjugate heat transfer between regions, buoyancy effects, turbulence,
reactions and radiation modelling.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "compressibleMomentumTransportModels.H"
#include "fluidMulticomponentThermophysicalTransportModel.H"
#include "fluidMulticomponentThermo.H"
#include "combustionModel.H"
#include "fixedGradientFvPatchFields.H"
#include "regionProperties.H"
#include "compressibleCourantNo.H"
#include "solidRegionDiffNo.H"
#include "solidThermo.H"
#include "fvModels.H"
#include "fvConstraints.H"
#include "coordinateSystem.H"
#include "pimpleMultiRegionControl.H"
#include "pressureReference.H"
#include "hydrostaticInitialisation.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#define NO_CONTROL
#define CREATE_MESH createMeshesPostProcess.H
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMeshes.H"
pimpleMultiRegionControl pimples(fluidRegions, solidRegions);
#include "createFields.H"
#include "initContinuityErrs.H"
#include "createFluidPressureControls.H"
#include "createTimeControls.H"
#include "readSolidTimeControls.H"
#include "compressibleMultiRegionCourantNo.H"
#include "solidRegionDiffusionNo.H"
#include "setInitialMultiRegionDeltaT.H"
while (pimples.run(runTime))
{
#include "readTimeControls.H"
#include "readSolidTimeControls.H"
#include "compressibleMultiRegionCourantNo.H"
#include "solidRegionDiffusionNo.H"
#include "setMultiRegionDeltaT.H"
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// Optional number of energy correctors
const int nEcorr = pimples.dict().lookupOrDefault<int>
(
"nEcorrectors",
1
);
// --- PIMPLE loop
while (pimples.loop())
{
List<tmp<fvVectorMatrix>> UEqns(fluidRegions.size());
for(int Ecorr=0; Ecorr<nEcorr; Ecorr++)
{
forAll(solidRegions, i)
{
Info<< "\nSolving for solid region "
<< solidRegions[i].name() << endl;
#include "setRegionSolidFields.H"
#include "solveSolid.H"
}
forAll(fluidRegions, i)
{
Info<< "\nSolving for fluid region "
<< fluidRegions[i].name() << endl;
#include "setRegionFluidFields.H"
#include "solveFluid.H"
}
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,2 +0,0 @@
#include "createFluidFields.H"
#include "createSolidFields.H"

View File

@ -1,4 +0,0 @@
regionProperties rp(runTime);
#include "createFluidMeshes.H"
#include "createSolidMeshes.H"

View File

@ -1,9 +0,0 @@
#include "createMeshes.H"
if (!fluidRegions.size() && !solidRegions.size())
{
FatalErrorIn(args.executable())
<< "No region meshes present" << exit(FatalError);
}
fvMesh& mesh = fluidRegions.size() ? fluidRegions[0] : solidRegions[0];

View File

@ -1,39 +0,0 @@
{
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he)
+ (
mvConvection.valid()
? mvConvection->fvmDiv(phi, he)
: fvm::div(phi, he)
)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? mvConvection.valid()
? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
: fvc::div(fvc::absolute(phi, rho, U), p/rho)
: -dpdt
)
+ thermophysicalTransport.divq(he)
==
rho*(U&g)
+ reaction.Qdot()
+ fvModels.source(rho, he)
);
EEqn.relax();
fvConstraints.constrain(EEqn);
EEqn.solve();
fvConstraints.constrain(he);
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;
}

View File

@ -1,38 +0,0 @@
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
UEqns[i] =
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence.divDevTau(U)
==
fvModels.source(rho, U)
);
fvVectorMatrix& UEqn = UEqns[i].ref();
UEqn.relax();
fvConstraints.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
fvConstraints.constrain(U);
K = 0.5*magSqr(U);
}
fvConstraints.constrain(U);

View File

@ -1,42 +0,0 @@
if (Y.size())
{
mvConvection = tmp<fv::convectionScheme<scalar>>
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.schemes().div("div(phi,Yi_h)")
)
);
}
reaction.correct();
forAll(Y, i)
{
if (composition.solve(i))
{
volScalarField& Yi = Y[i];
fvScalarMatrix YiEqn
(
fvm::ddt(rho, Yi) + mvConvection->fvmDiv(phi, Yi)
+ thermophysicalTransport.divj(Yi)
==
reaction.R(Yi)
+ fvModels.source(rho, Yi)
);
YiEqn.relax();
fvConstraints.constrain(YiEqn);
YiEqn.solve("Yi");
fvConstraints.constrain(Yi);
}
}
composition.normalise();

View File

@ -1,55 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "compressibleCourantNo.H"
#include "fvc.H"
Foam::scalar Foam::compressibleCourantNo
(
const fvMesh& mesh,
const Time& runTime,
const volScalarField& rho,
const surfaceScalarField& phi
)
{
scalarField sumPhi
(
fvc::surfaceSum(mag(phi))().primitiveField()
/ rho.primitiveField()
);
scalar CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
scalar meanCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
<< " max: " << CoNum << endl;
return CoNum;
}
// ************************************************************************* //

View File

@ -1,16 +0,0 @@
scalar CoNum = -great;
forAll(fluidRegions, regionI)
{
CoNum = max
(
compressibleCourantNo
(
fluidRegions[regionI],
runTime,
rhoFluid[regionI],
phiFluid[regionI]
),
CoNum
);
}

View File

@ -1,287 +0,0 @@
// Initialise fluid field pointer lists
PtrList<fluidMulticomponentThermo> thermoFluid(fluidRegions.size());
PtrList<volScalarField> rhoFluid(fluidRegions.size());
PtrList<volVectorField> UFluid(fluidRegions.size());
PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
PtrList<uniformDimensionedScalarField> hRefFluid(fluidRegions.size());
PtrList<uniformDimensionedScalarField> pRefFluid(fluidRegions.size());
PtrList<volScalarField> ghFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
PtrList<compressible::momentumTransportModel>
turbulenceFluid(fluidRegions.size());
PtrList<fluidMulticomponentThermophysicalTransportModel>
thermophysicalTransportFluid(fluidRegions.size());
PtrList<combustionModel> reactionFluid(fluidRegions.size());
PtrList<volScalarField> p_rghFluid(fluidRegions.size());
PtrList<volScalarField> KFluid(fluidRegions.size());
PtrList<volScalarField> dpdtFluid(fluidRegions.size());
PtrList<multivariateSurfaceInterpolationScheme<scalar>::fieldTable>
fieldsFluid(fluidRegions.size());
List<scalar> initialMassFluid(fluidRegions.size());
PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
PtrList<Foam::fvModels> fvModelsFluid(fluidRegions.size());
PtrList<fvConstraints> fvConstraintsFluid(fluidRegions.size());
// Populate fluid field pointer lists
forAll(fluidRegions, i)
{
Info<< "*** Reading fluid mesh thermophysical properties for region "
<< fluidRegions[i].name() << nl << endl;
Info<< " Adding to thermoFluid\n" << endl;
thermoFluid.set(i, fluidMulticomponentThermo::New(fluidRegions[i]).ptr());
Info<< " Adding to rhoFluid\n" << endl;
rhoFluid.set
(
i,
new volScalarField
(
IOobject
(
"rho",
runTime.timeName(),
fluidRegions[i],
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermoFluid[i].rho()
)
);
Info<< " Adding to UFluid\n" << endl;
UFluid.set
(
i,
new volVectorField
(
IOobject
(
"U",
runTime.timeName(),
fluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
fluidRegions[i]
)
);
Info<< " Adding to phiFluid\n" << endl;
phiFluid.set
(
i,
new surfaceScalarField
(
IOobject
(
"phi",
runTime.timeName(),
fluidRegions[i],
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rhoFluid[i]*UFluid[i])
& fluidRegions[i].Sf()
)
);
Info<< " Adding to gFluid\n" << endl;
gFluid.set
(
i,
new uniformDimensionedVectorField
(
IOobject
(
"g",
runTime.constant(),
fluidRegions[i],
IOobject::MUST_READ,
IOobject::NO_WRITE
)
)
);
Info<< " Adding to hRefFluid\n" << endl;
hRefFluid.set
(
i,
new uniformDimensionedScalarField
(
IOobject
(
"hRef",
runTime.constant(),
fluidRegions[i],
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
dimensionedScalar(dimLength, 0)
)
);
Info<< " Adding to pRefFluid\n" << endl;
pRefFluid.set
(
i,
new uniformDimensionedScalarField
(
IOobject
(
"pRef",
runTime.constant(),
fluidRegions[i],
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
dimensionedScalar(dimPressure, 0)
)
);
dimensionedScalar ghRef(- mag(gFluid[i])*hRefFluid[i]);
Info<< " Adding to ghFluid\n" << endl;
ghFluid.set
(
i,
new volScalarField
(
"gh",
(gFluid[i] & fluidRegions[i].C()) - ghRef
)
);
Info<< " Adding to ghfFluid\n" << endl;
ghfFluid.set
(
i,
new surfaceScalarField
(
"ghf",
(gFluid[i] & fluidRegions[i].Cf()) - ghRef
)
);
Info<< " Adding to turbulenceFluid\n" << endl;
turbulenceFluid.set
(
i,
compressible::momentumTransportModel::New
(
rhoFluid[i],
UFluid[i],
phiFluid[i],
thermoFluid[i]
).ptr()
);
Info<< " Adding to thermophysicalTransport\n" << endl;
thermophysicalTransportFluid.set
(
i,
fluidMulticomponentThermophysicalTransportModel::New
(
turbulenceFluid[i],
thermoFluid[i]
).ptr()
);
Info<< " Adding to reactionFluid\n" << endl;
reactionFluid.set
(
i,
combustionModel::New(thermoFluid[i], turbulenceFluid[i])
);
p_rghFluid.set
(
i,
new volScalarField
(
IOobject
(
"p_rgh",
runTime.timeName(),
fluidRegions[i],
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
fluidRegions[i]
)
);
fluidRegions[i].schemes().setFluxRequired(p_rghFluid[i].name());
initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
Info<< " Adding to KFluid\n" << endl;
KFluid.set
(
i,
new volScalarField
(
"K",
0.5*magSqr(UFluid[i])
)
);
Info<< " Adding to dpdtFluid\n" << endl;
dpdtFluid.set
(
i,
new volScalarField
(
IOobject
(
"dpdt",
runTime.timeName(),
fluidRegions[i]
),
fluidRegions[i],
dimensionedScalar
(
thermoFluid[i].p().dimensions()/dimTime,
0
)
)
);
Info<< " Adding to fieldsFluid\n" << endl;
fieldsFluid.set
(
i,
new multivariateSurfaceInterpolationScheme<scalar>::fieldTable
);
forAll(thermoFluid[i].composition().Y(), j)
{
fieldsFluid[i].add(thermoFluid[i].composition().Y()[j]);
}
fieldsFluid[i].add(thermoFluid[i].he());
Info<< " Adding MRF\n" << endl;
MRFfluid.set
(
i,
new IOMRFZoneList(fluidRegions[i])
);
Info<< " Adding fvModelsFluid\n" << endl;
fvModelsFluid.set
(
i,
new Foam::fvModels(fluidRegions[i])
);
Info<< " Adding fvConstraintsFluid\n" << endl;
fvConstraintsFluid.set
(
i,
new fvConstraints(fluidRegions[i])
);
turbulenceFluid[i].validate();
}

View File

@ -1,27 +0,0 @@
const wordList fluidNames
(
rp.found("fluid") ? rp["fluid"] : wordList(0)
);
PtrList<fvMesh> fluidRegions(fluidNames.size());
forAll(fluidNames, i)
{
Info<< "Create fluid mesh for region " << fluidNames[i]
<< " for time = " << runTime.timeName() << nl << endl;
fluidRegions.set
(
i,
new fvMesh
(
IOobject
(
fluidNames[i],
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
}

View File

@ -1,28 +0,0 @@
PtrList<pressureReference> pressureReferenceFluid(fluidRegions.size());
forAll(fluidRegions, i)
{
pressureReferenceFluid.set
(
i,
new pressureReference
(
p_rghFluid[i],
pimples.pimple(i).dict(),
false
)
);
hydrostaticInitialisation
(
p_rghFluid[i],
thermoFluid[i].p(),
rhoFluid[i],
UFluid[i],
ghFluid[i],
ghfFluid[i],
pRefFluid[i],
thermoFluid[i],
pimples.pimple(i).dict()
);
}

View File

@ -1 +0,0 @@
List<scalar> cumulativeContErrs(fluidRegions.size(), 0.0);

View File

@ -1,56 +0,0 @@
const fvMesh& mesh = fluidRegions[i];
combustionModel& reaction = reactionFluid[i];
fluidMulticomponentThermo& thermo = thermoFluid[i];
thermo.validate(args.executable(), "h", "e");
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
volScalarField& rho = rhoFluid[i];
volVectorField& U = UFluid[i];
surfaceScalarField& phi = phiFluid[i];
compressible::momentumTransportModel& turbulence = turbulenceFluid[i];
fluidMulticomponentThermophysicalTransportModel& thermophysicalTransport =
thermophysicalTransportFluid[i];
volScalarField& K = KFluid[i];
volScalarField& dpdt = dpdtFluid[i];
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& p_rgh = p_rghFluid[i];
const dimensionedVector& g = gFluid[i];
const dimensionedScalar& pRef = pRefFluid[i];
const volScalarField& gh = ghFluid[i];
const surfaceScalarField& ghf = ghfFluid[i];
multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields =
fieldsFluid[i];
IOMRFZoneList& MRF = MRFfluid[i];
Foam::fvModels& fvModels = fvModelsFluid[i];
Foam::fvConstraints& fvConstraints = fvConstraintsFluid[i];
#include "checkRadiationModel.H"
const dimensionedScalar initialMass
(
"initialMass",
dimMass,
initialMassFluid[i]
);
pimpleNoLoopControl& pimple = pimples.pimple(i);
pressureReference& pressureReference = pressureReferenceFluid[i];
scalar cumulativeContErr = cumulativeContErrs[i];
// This solver does not support moving mesh but it uses the pressure
// equation of one which does, so we need a dummy face-momentum field
autoPtr<surfaceVectorField> rhoUf(nullptr);

View File

@ -1,69 +0,0 @@
if (!pimple.flow())
{
if (pimple.models())
{
fvModels.correct();
}
if (pimple.thermophysics())
{
tmp<fv::convectionScheme<scalar>> mvConvection(nullptr);
if (Ecorr == 0)
{
#include "YEqn.H"
}
#include "EEqn.H"
}
}
else
{
if (Ecorr == 0)
{
if (!mesh.schemes().steady() && pimples.firstPimpleIter())
{
#include "rhoEqn.H"
}
if (pimple.models())
{
fvModels.correct();
}
#include "UEqn.H"
}
if (pimple.thermophysics())
{
tmp<fv::convectionScheme<scalar>> mvConvection(nullptr);
if (Ecorr == 0)
{
#include "YEqn.H"
}
#include "EEqn.H"
}
if (Ecorr == nEcorr - 1)
{
tmp<fvVectorMatrix>& tUEqn = UEqns[i];
fvVectorMatrix& UEqn = tUEqn.ref();
// --- PISO loop
while (pimple.correct())
{
#include "../../buoyantFoam/pEqn.H"
}
if (pimples.pimpleTurbCorr(i))
{
turbulence.correct();
thermophysicalTransport.correct();
}
if (!mesh.schemes().steady() && pimples.finalPimpleIter())
{
rho = thermo.rho();
}
}
}

View File

@ -1,28 +0,0 @@
// Initialise solid field pointer lists
PtrList<solidThermo> thermoSolid(solidRegions.size());
PtrList<Foam::fvModels> fvModelsSolid(solidRegions.size());
PtrList<fvConstraints> fvConstraintsSolid(solidRegions.size());
// Populate solid field pointer lists
forAll(solidRegions, i)
{
Info<< "*** Reading solid mesh thermophysical properties for region "
<< solidRegions[i].name() << nl << endl;
Info<< " Adding to thermoSolid\n" << endl;
thermoSolid.set(i, solidThermo::New(solidRegions[i]));
Info<< " Adding to fvModelsSolid\n" << endl;
fvModelsSolid.set
(
i,
new Foam::fvModels(solidRegions[i])
);
Info<< " Adding fvConstraintsSolid\n" << endl;
fvConstraintsSolid.set
(
i,
new fvConstraints(solidRegions[i])
);
}

View File

@ -1,32 +0,0 @@
const wordList solidNames
(
rp.found("solid") ? rp["solid"] : wordList(0)
);
PtrList<fvMesh> solidRegions(solidNames.size());
forAll(solidNames, i)
{
Info<< "Create solid mesh for region " << solidNames[i]
<< " for time = " << runTime.timeName() << nl << endl;
solidRegions.set
(
i,
new fvMesh
(
IOobject
(
solidNames[i],
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
// Force calculation of geometric properties to prevent it being done
// later in e.g. some boundary evaluation
//(void)solidRegions[i].weights();
//(void)solidRegions[i].deltaCoeffs();
}

View File

@ -1,15 +0,0 @@
const fvMesh& mesh = solidRegions[i];
solidThermo& thermo = thermoSolid[i];
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
volScalarField& e = thermo.he();
const Foam::fvModels& fvModels = fvModelsSolid[i];
Foam::fvConstraints& fvConstraints = fvConstraintsSolid[i];
#include "checkRadiationModel.H"
solidNoLoopControl& pimple = pimples.solid(i);

View File

@ -1,37 +0,0 @@
scalar DiNum = -great;
forAll(solidRegions, i)
{
// Note: do not use setRegionSolidFields.H to avoid double registering Cp
//#include "setRegionSolidFields.H"
const solidThermo& thermo = thermoSolid[i];
tmp<volScalarField> magKappa;
if (thermo.isotropic())
{
magKappa = thermo.kappa();
}
else
{
magKappa = mag(thermo.Kappa());
}
tmp<volScalarField> tcp = thermo.Cp();
const volScalarField& cp = tcp();
tmp<volScalarField> trho = thermo.rho();
const volScalarField& rho = trho();
DiNum = max
(
solidRegionDiffNo
(
solidRegions[i],
runTime,
rho*cp,
magKappa()
),
DiNum
);
}

View File

@ -1,25 +0,0 @@
{
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix eEqn
(
fvm::ddt(rho, e)
+ thermo.divq(e)
==
fvModels.source(rho, e)
);
eEqn.relax();
fvConstraints.constrain(eEqn);
eEqn.solve();
fvConstraints.constrain(e);
}
}
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;

View File

@ -0,0 +1,11 @@
#!/bin/sh
cd ${0%/*} || exit 1 # Run from this directory
# Parse arguments for library compilation
. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
wmake $targetType isothermalFluid
wmake $targetType fluid
wmake $targetType multicomponentFluid
#------------------------------------------------------------------------------

View File

@ -0,0 +1,4 @@
fluid.C
thermophysicalPredictor.C
LIB = $(FOAM_LIBBIN)/libfluid

View File

@ -1,4 +1,5 @@
EXE_INC = \
-I../isothermalFluid/lnInclude \
-I$(LIB_SRC)/physicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
@ -9,9 +10,9 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
LIB_LIBS = \
-lisothermalFluid \
-lfluidThermophysicalModels \
-lspecie \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2018-2019 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,61 +23,48 @@ License
\*---------------------------------------------------------------------------*/
#include "solidNoLoopControl.H"
#include "fluid.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(solidNoLoopControl, 0);
namespace solvers
{
defineTypeNameAndDebug(fluid, 0);
addToRunTimeSelectionTable(solver, fluid, fvMesh);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solidNoLoopControl::solidNoLoopControl
(
fvMesh& mesh,
const word& algorithmName,
const pimpleLoop& loop
)
Foam::solvers::fluid::fluid(fvMesh& mesh)
:
nonOrthogonalSolutionControl(mesh, algorithmName),
singleRegionConvergenceControl
isothermalFluid(mesh),
thermophysicalTransport
(
static_cast<singleRegionSolutionControl&>(*this)
),
singleRegionCorrectorConvergenceControl
(
static_cast<singleRegionSolutionControl&>(*this),
"outerCorrector"
),
loop_(loop)
{
read();
}
fluidThermophysicalTransportModel::New(turbulence(), thermo)
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::solidNoLoopControl::~solidNoLoopControl()
Foam::solvers::fluid::~fluid()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::solidNoLoopControl::read()
void Foam::solvers::fluid::thermophysicalTransportCorrector()
{
return
nonOrthogonalSolutionControl::read()
&& readResidualControls()
&& readCorrResidualControls();
}
bool Foam::solidNoLoopControl::isFinal() const
{
return loop_.finalPimpleIter() && finalNonOrthogonalIter();
if (pimple.transportCorr())
{
thermophysicalTransport->correct();
}
}

View File

@ -0,0 +1,122 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::solvers::fluid
Description
Solver module for steady or transient turbulent flow of compressible fluids
with heat-transfer for HVAC and similar applications, with optional
mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
Reference:
\verbatim
Greenshields, C. J., & Weller, H. G. (2022).
Notes on Computational Fluid Dynamics: General Principles.
CFD Direct Ltd.: Reading, UK.
\endverbatim
SourceFiles
fluid.C
\*---------------------------------------------------------------------------*/
#ifndef fluid_H
#define fluid_H
#include "isothermalFluid.H"
#include "fluidThermophysicalTransportModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
/*---------------------------------------------------------------------------*\
Class fluid Declaration
\*---------------------------------------------------------------------------*/
class fluid
:
public isothermalFluid
{
protected:
// Thermophysical transport
autoPtr<fluidThermophysicalTransportModel> thermophysicalTransport;
public:
//- Runtime type information
TypeName("fluid");
// Constructors
//- Construct from region mesh
fluid(fvMesh& mesh);
//- Disallow default bitwise copy construction
fluid(const fluid&) = delete;
//- Destructor
virtual ~fluid();
// Member Functions
//- Construct and solve the energy equation,
// convert to temperature
// and update thermophysical and transport properties
virtual void thermophysicalPredictor();
//- Correct the thermophysical transport modelling
virtual void thermophysicalTransportCorrector();
// Member Operators
//- Disallow default bitwise assignment
void operator=(const fluid&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,34 +23,42 @@ License
\*---------------------------------------------------------------------------*/
#include "solidRegionDiffNo.H"
#include "surfaceInterpolate.H"
#include "fluid.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar Foam::solidRegionDiffNo
(
const fvMesh& mesh,
const Time& runTime,
const volScalarField& Cprho,
const volScalarField& kappa
)
void Foam::solvers::fluid::thermophysicalPredictor()
{
surfaceScalarField kapparhoCpbyDelta
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
sqr(mesh.surfaceInterpolation::deltaCoeffs())
*fvc::interpolate(kappa)
/fvc::interpolate(Cprho)
fvm::ddt(rho, he) + fvm::div(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? fvc::div(fvc::absolute(phi, rho, U), p/rho)
: -dpdt
)
+ thermophysicalTransport->divq(he)
==
(
buoyancy.valid()
? fvModels().source(rho, he) + rho*(U & buoyancy->g)
: fvModels().source(rho, he)
)
);
const scalar DiNum = max(kapparhoCpbyDelta).value()*runTime.deltaTValue();
const scalar meanDiNum =
average(kapparhoCpbyDelta).value()*runTime.deltaTValue();
EEqn.relax();
Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
<< " max: " << DiNum << endl;
fvConstraints().constrain(EEqn);
return DiNum;
EEqn.solve();
fvConstraints().constrain(he);
thermo.correct();
}
// ************************************************************************* //

View File

@ -0,0 +1,11 @@
buoyancy/buoyancy.C
setRDeltaT.C
moveMesh.C
prePredictor.C
correctDensity.C
momentumPredictor.C
correctPressure.C
correctBuoyantPressure.C
isothermalFluid.C
LIB = $(FOAM_LIBBIN)/libisothermalFluid

View File

@ -3,20 +3,16 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
LIB_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \
-lfiniteVolume \
-lsampling \
-lmeshTools \
-lsampling \
-lfvModels \
-lfvConstraints

View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "buoyancy.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solvers::buoyancy::buoyancy(const fvMesh& mesh_)
:
mesh(mesh_),
runTime(mesh.time()),
g
(
IOobject
(
"g",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
hRef
(
IOobject
(
"hRef",
runTime.constant(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
dimensionedScalar(dimLength, 0)
),
pRef
(
IOobject
(
"pRef",
runTime.constant(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
dimensionedScalar(dimPressure, 0)
),
ghRef(-mag(g)*hRef),
gh("gh", (g & mesh.C()) - ghRef),
ghf("ghf", (g & mesh.Cf()) - ghRef),
p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
{
mesh.schemes().setFluxRequired(p_rgh.name());
}
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::solvers::buoyancy> Foam::solvers::buoyancy::New
(
const fvMesh& mesh
)
{
return typeIOobject<volScalarField>
(
"p_rgh",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ
).headerOk()
? autoPtr<buoyancy>(new solvers::buoyancy(mesh))
: autoPtr<buoyancy>(nullptr);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::solvers::buoyancy::~buoyancy()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::buoyancy::moveMesh()
{
gh = (g & mesh.C()) - ghRef;
ghf = (g & mesh.Cf()) - ghRef;
}
// ************************************************************************* //

View File

@ -0,0 +1,146 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::solvers::buoyancy
Description
Buoyancy related data for the Foam::solvers::isothermalFluid solver module
when solving buoyant cases with \c p_rgh and is selected based on the
presence of the \c p_rgh field file.
Provides:
g : gravitational acceleration
hRef : optional reference height
pRef : optional reference pressure
ghRef : -mag(g)*hRef
gh : (g & h) - ghRef
ghf : (g & hf) - ghRef
p_rgh : p - rho*gh - pRef
The \c Foam::solvers::buoyancy::New function returns an \c autoPtr<buoyancy>
pointer containing either a \c buoyancy class pointer if the \c p_rgh field
file is available othewise a \c nullptr which can be checked using
the \c autoPtr::valid() member function.
SourceFiles
buoyancy.C
\*---------------------------------------------------------------------------*/
#ifndef buoyancy_H
#define buoyancy_H
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "uniformDimensionedFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
/*---------------------------------------------------------------------------*\
Class buoyancy Declaration
\*---------------------------------------------------------------------------*/
class buoyancy
{
// Private member data
const fvMesh& mesh;
const Time& runTime;
public:
// Thermophysical properties
//- Gravitational acceleration
uniformDimensionedVectorField g;
//- Reference height
uniformDimensionedScalarField hRef;
//- Reference pressure
uniformDimensionedScalarField pRef;
//- -mag(g)*hRef
dimensionedScalar ghRef;
//- (g & h) - ghRef
volScalarField gh;
//- (g & hf) - ghRef
surfaceScalarField ghf;
//- Buoyant pressure p - rho*gh - pRef
volScalarField p_rgh;
// Constructors
//- Construct from the region mesh
buoyancy(const fvMesh& mesh);
//- Disallow default bitwise copy construction
buoyancy(const buoyancy&) = delete;
// Selectors
//- Select, construct and return the buoyancy
static autoPtr<buoyancy> New(const fvMesh& mesh);
//- Destructor
virtual ~buoyancy();
// Member Functions
//- Update gh and ghf following mesh-motion
void moveMesh();
// Member Operators
//- Disallow default bitwise assignment
void operator=(const buoyancy&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,262 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "isothermalFluid.H"
#include "constrainHbyA.H"
#include "constrainPressure.H"
#include "adjustPhi.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::correctBuoyantPressure()
{
// Local references to the buoyancy paramaters
const volScalarField& gh = buoyancy->gh;
const surfaceScalarField& ghf = buoyancy->ghf;
const uniformDimensionedScalarField pRef = buoyancy->pRef;
const volScalarField& psi = thermo.psi();
rho = thermo.rho();
rho.relax();
fvVectorMatrix& UEqn = tUEqn.ref();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
const volScalarField rAU("rAU", 1.0/UEqn.A());
const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
tmp<volScalarField> rAtU
(
pimple.consistent()
? volScalarField::New("rAtU", 1.0/(1.0/rAU - UEqn.H1()))
: tmp<volScalarField>(nullptr)
);
tmp<surfaceScalarField> rhorAtUf
(
pimple.consistent()
? surfaceScalarField::New("rhoRAtUf", fvc::interpolate(rho*rAtU()))
: tmp<surfaceScalarField>(nullptr)
);
const volScalarField& rAAtU = pimple.consistent() ? rAtU() : rAU;
const surfaceScalarField& rhorAAtUf =
pimple.consistent() ? rhorAtUf() : rhorAUf;
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
if (pimple.nCorrPiso() <= 1)
{
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
const bool adjustMass =
mesh.schemes().steady() && adjustPhi(phiHbyA, U, p_rgh);
const surfaceScalarField ghGradRhof(-ghf*fvc::snGrad(rho)*mesh.magSf());
phiHbyA += rhorAUf*ghGradRhof;
tmp<fvScalarMatrix> tp_rghEqn;
if (pimple.transonic())
{
const surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))
*fvc::relative(phiHbyA, rho, U)
);
const fvScalarMatrix divPhidp(fvm::div(phid, p));
phiHbyA -= divPhidp.flux();
if (pimple.consistent())
{
const surfaceScalarField gradpf(fvc::snGrad(p_rgh)*mesh.magSf());
phiHbyA += (rhorAAtUf - rhorAUf)*gradpf;
HbyA += (rAAtU - rAU)*fvc::reconstruct(gradpf - ghGradRhof);
}
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAAtUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA) + divPhidp
==
fvModels().source(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
tp_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAAtUf, p);
fvScalarMatrix& p_rghEqn = tp_rghEqn.ref();
// Relax the pressure equation to ensure diagonal-dominance
p_rghEqn.relax();
p_rghEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
p_rghEqn.solve();
}
}
else
{
if (pimple.consistent())
{
const surfaceScalarField gradpf(fvc::snGrad(p_rgh)*mesh.magSf());
phiHbyA += (rhorAAtUf - rhorAUf)*gradpf;
HbyA += (rAAtU - rAU)*fvc::reconstruct(gradpf - ghGradRhof);
}
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAAtUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
==
fvModels().source(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
tp_rghEqn = p_rghDDtEqn - fvm::laplacian(rhorAAtUf, p_rgh);
fvScalarMatrix& p_rghEqn = tp_rghEqn.ref();
p_rghEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
p_rghEqn.solve();
}
}
const fvScalarMatrix& p_rghEqn = tp_rghEqn();
phi = phiHbyA + p_rghEqn.flux();
if (!mesh.schemes().steady())
{
p = p_rgh + rho*gh + pRef;
const bool constrained = fvConstraints().constrain(p);
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (constrained)
{
rho = thermo.rho();
}
correctDensity();
}
continuityErrors();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
p = p_rgh + rho*gh + pRef;
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAAtU*fvc::reconstruct((ghGradRhof + p_rghEqn.flux()/rhorAAtUf));
U.correctBoundaryConditions();
fvConstraints().constrain(U);
K = 0.5*magSqr(U);
if (mesh.schemes().steady())
{
if (fvConstraints().constrain(p))
{
p_rgh = p - rho*gh - pRef;
p_rgh.correctBoundaryConditions();
}
}
// For steady compressible closed-volume cases adjust the pressure level
// to obey overall mass continuity
if (adjustMass && !thermo.incompressible())
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/fvc::domainIntegrate(psi);
p_rgh = p - rho*gh - pRef;
p_rgh.correctBoundaryConditions();
}
if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass)
{
rho = thermo.rho();
}
// Correct rhoUf if the mesh is moving
fvc::correctRhoUf(rhoUf, rho, U, phi, MRF);
if (mesh.schemes().steady() || pimple.simpleRho())
{
rho.relax();
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,28 +21,27 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Calculates and outputs the mean and maximum Diffusion Numbers for the solid
regions
\*---------------------------------------------------------------------------*/
#ifndef solidRegionDiffNo_H
#define solidRegionDiffNo_H
#include "isothermalFluid.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
namespace Foam
void Foam::solvers::isothermalFluid::correctDensity()
{
scalar solidRegionDiffNo
fvScalarMatrix rhoEqn
(
const fvMesh& mesh,
const Time& runTime,
const volScalarField& Cprho,
const volScalarField& kappa
fvm::ddt(rho) + fvc::div(phi)
==
fvModels().source(rho)
);
fvConstraints().constrain(rhoEqn);
rhoEqn.solve();
fvConstraints().constrain(rho);
}
#endif
// ************************************************************************* //

View File

@ -0,0 +1,245 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "isothermalFluid.H"
#include "constrainHbyA.H"
#include "constrainPressure.H"
#include "adjustPhi.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::correctPressure()
{
const volScalarField& psi = thermo.psi();
rho = thermo.rho();
rho.relax();
fvVectorMatrix& UEqn = tUEqn.ref();
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
const volScalarField rAU("rAU", 1.0/UEqn.A());
const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
tmp<volScalarField> rAtU
(
pimple.consistent()
? volScalarField::New("rAtU", 1.0/(1.0/rAU - UEqn.H1()))
: tmp<volScalarField>(nullptr)
);
tmp<surfaceScalarField> rhorAtUf
(
pimple.consistent()
? surfaceScalarField::New("rhoRAtUf", fvc::interpolate(rho*rAtU()))
: tmp<surfaceScalarField>(nullptr)
);
const volScalarField& rAAtU = pimple.consistent() ? rAtU() : rAU;
const surfaceScalarField& rhorAAtUf =
pimple.consistent() ? rhorAtUf() : rhorAUf;
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
if (pimple.nCorrPiso() <= 1)
{
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf))
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
bool adjustMass = false;
if (pimple.transonic())
{
const surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/fvc::interpolate(rho))
*fvc::relative(phiHbyA, rho, U)
);
const fvScalarMatrix divPhidp(fvm::div(phid, p));
phiHbyA -= divPhidp.flux();
if (pimple.consistent())
{
phiHbyA += (rhorAAtUf - rhorAUf)*fvc::snGrad(p)*mesh.magSf();
HbyA += (rAAtU - rAU)*fvc::grad(p);
}
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA) + divPhidp
==
fvModels().source(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAAtUf, p));
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
if (pimple.consistent())
{
phiHbyA += (rhorAAtUf - rhorAUf)*fvc::snGrad(p)*mesh.magSf();
HbyA += (rAAtU - rAU)*fvc::grad(p);
}
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
fvc::makeRelative(phiHbyA, rho, U);
if (mesh.schemes().steady())
{
adjustMass = adjustPhi(phiHbyA, U, p);
}
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
fvModels().source(psi, p, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAAtUf, p));
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
if (!mesh.schemes().steady())
{
const bool constrained = fvConstraints().constrain(p);
// Thermodynamic density update
thermo.correctRho(psi*p - psip0);
if (constrained)
{
rho = thermo.rho();
}
correctDensity();
}
continuityErrors();
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAAtU*fvc::grad(p);
U.correctBoundaryConditions();
fvConstraints().constrain(U);
K = 0.5*magSqr(U);
if (mesh.schemes().steady())
{
fvConstraints().constrain(p);
}
// For steady compressible closed-volume cases adjust the pressure level
// to obey overall mass continuity
if (adjustMass && !thermo.incompressible())
{
p += (initialMass - fvc::domainIntegrate(thermo.rho()))
/fvc::domainIntegrate(psi);
p.correctBoundaryConditions();
}
if (mesh.schemes().steady() || pimple.simpleRho() || adjustMass)
{
rho = thermo.rho();
}
// Correct rhoUf if the mesh is moving
fvc::correctRhoUf(rhoUf, rho, U, phi, MRF);
if (mesh.schemes().steady() || pimple.simpleRho())
{
rho.relax();
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,412 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "isothermalFluid.H"
#include "localEulerDdtScheme.H"
#include "hydrostaticInitialisation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
defineTypeNameAndDebug(isothermalFluid, 0);
addToRunTimeSelectionTable(solver, isothermalFluid, fvMesh);
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::read()
{
maxCo =
runTime.controlDict().lookupOrDefault<scalar>("maxCo", 1.0);
maxDeltaT_ =
runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", great);
correctPhi = pimple.dict().lookupOrDefault
(
"correctPhi",
correctPhi
);
checkMeshCourantNo = pimple.dict().lookupOrDefault
(
"checkMeshCourantNo",
checkMeshCourantNo
);
moveMeshOuterCorrectors = pimple.dict().lookupOrDefault
(
"moveMeshOuterCorrectors",
moveMeshOuterCorrectors
);
}
void Foam::solvers::isothermalFluid::correctCoNum()
{
const scalarField sumPhi
(
fvc::surfaceSum(mag(phi))().primitiveField()/rho.primitiveField()
);
CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
const scalar meanCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
Info<< "Courant Number mean: " << meanCoNum
<< " max: " << CoNum << endl;
}
void Foam::solvers::isothermalFluid::continuityErrors()
{
scalar sumLocalContErr = 0;
scalar globalContErr = 0;
if (mesh.schemes().steady())
{
const volScalarField contErr(fvc::div(phi));
sumLocalContErr =
runTime.deltaTValue()
*mag(contErr)().weightedAverage(mesh.V()).value();
globalContErr =
runTime.deltaTValue()
*contErr.weightedAverage(mesh.V()).value();
}
else
{
const dimensionedScalar totalMass = fvc::domainIntegrate(rho);
sumLocalContErr =
(fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass).value();
globalContErr =
(fvc::domainIntegrate(rho - thermo.rho())/totalMass).value();
}
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solvers::isothermalFluid::isothermalFluid
(
fvMesh& mesh,
autoPtr<fluidThermo> thermoPtr
)
:
solver(mesh),
thermo_(thermoPtr),
thermo(thermo_()),
p(thermo.p()),
rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
thermo.rho()
),
dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar(p.dimensions()/dimTime, 0)
),
buoyancy(buoyancy::New(mesh)),
p_rgh(buoyancy.valid() ? buoyancy->p_rgh : p),
pressureReference
(
p,
p_rgh,
pimple.dict(),
thermo.incompressible()
),
U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rho*U) & mesh.Sf()
),
K("K", 0.5*magSqr(U)),
turbulence
(
compressible::momentumTransportModel::New
(
rho,
U,
phi,
thermo
)
),
initialMass(fvc::domainIntegrate(rho)),
cumulativeContErr(0),
MRF(mesh),
CoNum(0)
{
// Read the controls
read();
thermo.validate("isothermalFluid", "h", "e");
mesh.schemes().setFluxRequired(p.name());
turbulence->validate();
if (buoyancy.valid())
{
hydrostaticInitialisation
(
p_rgh,
p,
rho,
U,
buoyancy->gh,
buoyancy->ghf,
buoyancy->pRef,
thermo,
pimple.dict()
);
}
if (mesh.dynamic())
{
Info<< "Constructing face momentum rhoUf" << endl;
rhoUf = new surfaceVectorField
(
IOobject
(
"rhoUf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(rho*U)
);
}
if (transient())
{
correctCoNum();
}
else if (LTS)
{
Info<< "Using LTS" << endl;
trDeltaT = tmp<volScalarField>
(
new volScalarField
(
IOobject
(
fv::localEulerDdt::rDeltaTName,
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless/dimTime, 1),
extrapolatedCalculatedFvPatchScalarField::typeName
)
);
}
}
Foam::solvers::isothermalFluid::isothermalFluid(fvMesh& mesh)
:
isothermalFluid(mesh, fluidThermo::New(mesh))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::solvers::isothermalFluid::~isothermalFluid()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar Foam::solvers::isothermalFluid::maxDeltaT() const
{
if (CoNum > small)
{
const scalar deltaT = maxCo*runTime.deltaTValue()/CoNum;
return min(min(deltaT, fvModels().maxDeltaT()), maxDeltaT_);
}
else
{
return runTime.deltaTValue();
}
}
void Foam::solvers::isothermalFluid::preSolve()
{
// Read the controls
read();
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
if (correctPhi)
{
divrhoU = new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
}
fvModels().preUpdateMesh();
// Store momentum to set rhoUf for introduced faces
if (rhoUf.valid())
{
rhoU = new volVectorField("rhoU", rho*U);
}
// Update the mesh for topology change, mesh to mesh mapping
mesh.update();
if (transient())
{
correctCoNum();
}
else if (LTS)
{
setRDeltaT();
}
}
void Foam::solvers::isothermalFluid::thermophysicalPredictor()
{
thermo.correct();
}
void Foam::solvers::isothermalFluid::pressureCorrector()
{
while (pimple.correct())
{
if (buoyancy.valid())
{
correctBuoyantPressure();
}
else
{
correctPressure();
}
}
tUEqn.clear();
}
void Foam::solvers::isothermalFluid::momentumTransportCorrector()
{
if (pimple.transportCorr())
{
turbulence->correct();
}
}
void Foam::solvers::isothermalFluid::thermophysicalTransportCorrector()
{}
void Foam::solvers::isothermalFluid::postSolve()
{
rhoU.clear();
divrhoU.clear();
if (!mesh.schemes().steady())
{
rho = thermo.rho();
}
}
// ************************************************************************* //

View File

@ -0,0 +1,308 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::solvers::isothermalFluid
Description
Solver module for steady or transient turbulent flow of compressible
isothermal fluids with optional mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
Reference:
\verbatim
Greenshields, C. J., & Weller, H. G. (2022).
Notes on Computational Fluid Dynamics: General Principles.
CFD Direct Ltd.: Reading, UK.
\endverbatim
SourceFiles
isothermalFluid.C
\*---------------------------------------------------------------------------*/
#ifndef isothermalFluid_H
#define isothermalFluid_H
#include "solver.H"
#include "fluidThermo.H"
#include "compressibleMomentumTransportModels.H"
#include "buoyancy.H"
#include "pressureReference.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
/*---------------------------------------------------------------------------*\
Class isothermalFluid Declaration
\*---------------------------------------------------------------------------*/
class isothermalFluid
:
public solver
{
protected:
// Control parameters
//- Maximum allowed Courant number
scalar maxCo;
//- Maximum time-step
scalar maxDeltaT_;
//- Switch to correct the flux after mesh change
bool correctPhi;
//- Switch to check the mesh Courant number after mesh change
bool checkMeshCourantNo;
//- Switch to move the mesh at the start of every PIMPLE outer corrected
// rather than the first corrector only which is the default
bool moveMeshOuterCorrectors;
// Thermophysical properties
//- Pointer to the fluid thermophysical properties
autoPtr<fluidThermo> thermo_;
//- Reference to the fluid thermophysical properties
fluidThermo& thermo;
//- Reference to the pressure field
volScalarField& p;
//- The continuity density field
volScalarField rho;
//- Rate of change of the pressure
// Used in the enthalpy equation
volScalarField dpdt;
// Optional buoyancy
//- Pointer to the optional buoyancy force
// Case is considered buoyant if the p_rgh field exists
autoPtr<solvers::buoyancy> buoyancy;
//- Reference to the buoyant pressure for buoyant cases
// otherwise to the pressure
volScalarField& p_rgh;
// Pressure reference
//- Pressure reference
Foam::pressureReference pressureReference;
// Kinematic properties
//- Velocity field
volVectorField U;
//- Mass-flux field
surfaceScalarField phi;
//- Kinetic energy field
// Used in the energy equation
volScalarField K;
// Momentum transport
//- Pointer to the momentum transport model
autoPtr<compressible::momentumTransportModel> turbulence;
// Continuity properties
//- Initial mass in the region
dimensionedScalar initialMass;
//- Current cumulative continuity error
scalar cumulativeContErr;
// Optional models
//- MRF zone list
IOMRFZoneList MRF;
// Time-step controls
//- Current maximum Courant number
scalar CoNum;
//- Current mean Courant number
scalar meanCoNum;
// Cached temporary fields
//- Pointer to the surface momentum field
// used to recreate the flux after mesh-change
autoPtr<surfaceVectorField> rhoUf;
//- Pointer to the vol momentum field
// used for mesh-change to set rhoUf for introduced faces
autoPtr<volVectorField> rhoU;
//- Pointer to the vol momentum divergence field
// used in correctPhi to ensure the corrected phi has the
// same divergence
autoPtr<volScalarField> divrhoU;
//- Optional LTS reciprocal time-step field
tmp<volScalarField> trDeltaT;
//- Cached momentum matrix
// shared between the momentum predictor and pressure corrector
tmp<fvVectorMatrix> tUEqn;
private:
// Private Member Data
//- Pointer to the demand driven fvModels MeshObject
mutable Foam::fvModels* fvModelsPtr;
//- Pointer to the demand driven fvConstraints MeshObject
mutable Foam::fvConstraints* fvConstraintsPtr;
// Private Member Functions
//- Read controls
void read();
//- Set rDeltaT for LTS
virtual void setRDeltaT();
//- Correct the cached Courant numbers
void correctCoNum();
//- Check mesh Courant numbers for moving mesh cases
void meshCourantNo() const;
//- Calculate and print the continuity errors
void continuityErrors();
//- Construct the continuity equation and correct the density
void correctDensity();
//- Construct the pressure equation
// and correct the pressure and velocity
void correctPressure();
//- Construct the buoyant pressure equation
// and correct the pressure and velocity
void correctBuoyantPressure();
public:
//- Runtime type information
TypeName("isothermalFluid");
// Constructors
//- Construct from region mesh and thermophysical properties
isothermalFluid(fvMesh& mesh, autoPtr<fluidThermo>);
//- Construct from region mesh
isothermalFluid(fvMesh& mesh);
//- Disallow default bitwise copy construction
isothermalFluid(const isothermalFluid&) = delete;
//- Destructor
virtual ~isothermalFluid();
// Member Functions
//- Return the current maximum time-step for stable solution
virtual scalar maxDeltaT() const;
//- Called at the start of the time-step, before the PIMPLE loop
virtual void preSolve();
//- Called at the start of the PIMPLE loop to move the mesh
virtual bool moveMesh();
//- Called at the start of the PIMPLE loop
virtual void prePredictor();
//- Construct and optionally solve the momentum equation
virtual void momentumPredictor();
//- Construct and solve the energy equation,
// convert to temperature
// and update thermophysical and transport properties
virtual void thermophysicalPredictor();
//- Construct and solve the pressure equation in the PISO loop
virtual void pressureCorrector();
//- Correct the momentum transport modelling
// Newtonian, non-Newtonian or turbulent
virtual void momentumTransportCorrector();
//- Correct the thermophysical transport modelling
virtual void thermophysicalTransportCorrector();
//- Called after the PIMPLE loop at the end of the time-step
virtual void postSolve();
// Member Operators
//- Disallow default bitwise assignment
void operator=(const isothermalFluid&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,38 +21,56 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
setInitialMultiRegionDeltaT
Description
Set the initial timestep for the CHT MultiRegion solver.
\*---------------------------------------------------------------------------*/
if (adjustTimeStep)
#include "isothermalFluid.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::momentumPredictor()
{
if ((runTime.timeIndex() == 0) && ((CoNum > small) || (DiNum > small)))
MRF.correctBoundaryVelocity(U);
tUEqn =
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevTau(U)
==
fvModels().source(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvConstraints().constrain(UEqn);
if (pimple.momentumPredictor())
{
if (CoNum < small)
if (buoyancy.valid())
{
CoNum = small;
}
if (DiNum < small)
{
DiNum = small;
}
runTime.setDeltaT
(
min
solve
(
min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
min(runTime.deltaTValue(), maxDeltaT)
)
);
Info<< "deltaT = " << runTime.deltaT().value() << endl;
UEqn
==
fvc::reconstruct
(
(
- buoyancy->ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
)*mesh.magSf()
)
);
}
else
{
solve(UEqn == -fvc::grad(p));
}
fvConstraints().constrain(U);
K = 0.5*magSqr(U);
}
}
// ************************************************************************* //

View File

@ -0,0 +1,105 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "isothermalFluid.H"
#include "CorrectPhi.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::meshCourantNo() const
{
if (checkMeshCourantNo)
{
const scalarField sumPhi
(
fvc::surfaceSum(mag(mesh.phi()))().primitiveField()
);
const scalar meshCoNum
(
0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue()
);
const scalar meanMeshCoNum
(
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue()
);
Info<< "Mesh Courant Number mean: " << meanMeshCoNum
<< " max: " << meshCoNum << endl;
}
}
bool Foam::solvers::isothermalFluid::moveMesh()
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
// Move the mesh
mesh.move();
if (mesh.changing())
{
if (buoyancy.valid())
{
buoyancy->moveMesh();
}
MRF.update();
if (correctPhi)
{
// Calculate absolute flux
// from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
correctUphiBCs(rho, U, phi, true);
CorrectPhi
(
phi,
p,
rho,
thermo.psi(),
dimensionedScalar("rAUf", dimTime, 1),
divrhoU(),
pimple
);
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
}
meshCourantNo();
return true;
}
}
return false;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,14 +21,27 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
readSolidTimeControls
Description
Read the control parameters used in the solid
\*---------------------------------------------------------------------------*/
scalar maxDi = runTime.controlDict().lookupOrDefault<scalar>("maxDi", 10.0);
#include "isothermalFluid.H"
#include "CorrectPhi.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::prePredictor()
{
if
(
!mesh.schemes().steady()
&& !pimple.simpleRho()
&& pimple.firstIter()
)
{
correctDensity();
}
fvModels().correct();
}
// ************************************************************************* //

View File

@ -1,4 +1,37 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "isothermalFluid.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::solvers::isothermalFluid::setRDeltaT()
{
const volScalarField& psi = thermo.psi();
volScalarField& rDeltaT = trDeltaT.ref();
const dictionary& pimpleDict = pimple.dict();
@ -94,3 +127,6 @@
<< ", " << gMax(1/rDeltaT.primitiveField()) << endl;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,5 @@
multicomponentFluid.C
setRDeltaT.C
thermophysicalPredictor.C
LIB = $(FOAM_LIBBIN)/libmulticomponentFluid

View File

@ -1,6 +1,5 @@
EXE_INC = \
-I. \
-I$(FOAM_SOLVERS)/combustion/reactingFoam \
-I../isothermalFluid/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
-I$(LIB_SRC)/ThermophysicalTransportModels/lnInclude \
@ -12,22 +11,23 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude \
-I$(LIB_SRC)/combustionModels/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
LIB_LIBS = \
-lisothermalFluid \
-lfluidThermophysicalModels \
-lspecie \
-lchemistryModel \
-lODE \
-lcombustionModels \
-lmulticomponentThermophysicalModels \
-lmomentumTransportModels \
-lcompressibleMomentumTransportModels \
-lthermophysicalTransportModels \
-lfluidMulticomponentThermophysicalTransportModels \
-lmulticomponentThermophysicalModels \
-lspecie \
-lfluidThermophysicalModels \
-lchemistryModel \
-lODE \
-lcombustionModels \
-lfiniteVolume \
-lfvModels \
-lfvConstraints \

View File

@ -0,0 +1,94 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "multicomponentFluid.H"
#include "localEulerDdtScheme.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
defineTypeNameAndDebug(multicomponentFluid, 0);
addToRunTimeSelectionTable(solver, multicomponentFluid, fvMesh);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solvers::multicomponentFluid::multicomponentFluid(fvMesh& mesh)
:
isothermalFluid
(
mesh,
autoPtr<fluidThermo>(fluidMulticomponentThermo::New(mesh).ptr())
),
thermo(refCast<fluidMulticomponentThermo>(isothermalFluid::thermo)),
composition(thermo.composition()),
Y(composition.Y()),
reaction(combustionModel::New(thermo, turbulence())),
thermophysicalTransport
(
fluidMulticomponentThermophysicalTransportModel::New
(
turbulence(),
thermo
)
)
{
forAll(Y, i)
{
fields.add(Y[i]);
}
fields.add(thermo.he());
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::solvers::multicomponentFluid::~multicomponentFluid()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::multicomponentFluid::thermophysicalTransportCorrector()
{
if (pimple.transportCorr())
{
thermophysicalTransport->correct();
}
}
// ************************************************************************* //

View File

@ -0,0 +1,143 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::solvers::multicomponentFluid
Description
Solver module for steady or transient turbulent flow of compressible
reacting fluids with optional mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
Reference:
\verbatim
Greenshields, C. J., & Weller, H. G. (2022).
Notes on Computational Fluid Dynamics: General Principles.
CFD Direct Ltd.: Reading, UK.
\endverbatim
SourceFiles
multicomponentFluid.C
\*---------------------------------------------------------------------------*/
#ifndef multicomponentFluid_H
#define multicomponentFluid_H
#include "isothermalFluid.H"
#include "fluidMulticomponentThermo.H"
#include "combustionModel.H"
#include "fluidMulticomponentThermophysicalTransportModel.H"
#include "multivariateScheme.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
/*---------------------------------------------------------------------------*\
Class multicomponentFluid Declaration
\*---------------------------------------------------------------------------*/
class multicomponentFluid
:
public isothermalFluid
{
protected:
// Compositon
fluidMulticomponentThermo& thermo;
basicSpecieMixture& composition;
PtrList<volScalarField>& Y;
autoPtr<combustionModel> reaction;
multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
// Thermophysical transport
autoPtr<fluidMulticomponentThermophysicalTransportModel>
thermophysicalTransport;
private:
// Private Member Functions
//- Set rDeltaT for LTS
virtual void setRDeltaT();
public:
//- Runtime type information
TypeName("multicomponentFluid");
// Constructors
//- Construct from region mesh
multicomponentFluid(fvMesh& mesh);
//- Disallow default bitwise copy construction
multicomponentFluid(const multicomponentFluid&) = delete;
//- Destructor
virtual ~multicomponentFluid();
// Member Functions
//- Construct and solve the energy equation,
// convert to temperature
// and update thermophysical and transport properties
virtual void thermophysicalPredictor();
//- Correct the thermophysical transport modelling
virtual void thermophysicalTransportCorrector();
// Member Operators
//- Disallow default bitwise assignment
void operator=(const multicomponentFluid&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2013-2022 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,6 +23,12 @@ License
\*---------------------------------------------------------------------------*/
#include "multicomponentFluid.H"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::solvers::multicomponentFluid::setRDeltaT()
{
volScalarField& rDeltaT = trDeltaT.ref();
@ -70,7 +76,7 @@ License
{
volScalarField::Internal rDeltaTT
(
mag(reaction->Qdot())/(alphaTemp*rho*thermo.Cp()*T)
mag(reaction->Qdot())/(alphaTemp*rho*thermo.Cp()*thermo.T())
);
Info<< " Temperature = "
@ -137,7 +143,7 @@ License
}
else
{
IOWarningIn(args.executable().c_str(), Yref)
WarningInFunction
<< "Cannot find any active species in Yref " << Yref
<< endl;
}

View File

@ -0,0 +1,107 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "multicomponentFluid.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::multicomponentFluid::thermophysicalPredictor()
{
tmp<fv::convectionScheme<scalar>> mvConvection
(
fv::convectionScheme<scalar>::New
(
mesh,
fields,
phi,
mesh.schemes().div("div(phi,Yi_h)")
)
);
reaction->correct();
forAll(Y, i)
{
if (composition.solve(i))
{
volScalarField& Yi = Y[i];
fvScalarMatrix YiEqn
(
fvm::ddt(rho, Yi)
+ mvConvection->fvmDiv(phi, Yi)
+ thermophysicalTransport->divj(Yi)
==
reaction->R(Yi)
+ fvModels().source(rho, Yi)
);
YiEqn.relax();
fvConstraints().constrain(YiEqn);
YiEqn.solve("Yi");
fvConstraints().constrain(Yi);
}
}
composition.normalise();
volScalarField& he = thermo.he();
fvScalarMatrix EEqn
(
fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+ fvc::ddt(rho, K) + fvc::div(phi, K)
+ (
he.name() == "e"
? mvConvection->fvcDiv(fvc::absolute(phi, rho, U), p/rho)
: -dpdt
)
+ thermophysicalTransport->divq(he)
==
reaction->Qdot()
+ (
buoyancy.valid()
? fvModels().source(rho, he) + rho*(U & buoyancy->g)
: fvModels().source(rho, he)
)
);
EEqn.relax();
fvConstraints().constrain(EEqn);
EEqn.solve();
fvConstraints().constrain(he);
thermo.correct();
}
// ************************************************************************* //

View File

@ -0,0 +1,5 @@
regionSolvers/regionSolvers.C
setDeltaT.C
foamMultiRun.C
EXE = $(FOAM_APPBIN)/foamMultiRun

View File

@ -0,0 +1,6 @@
EXE_INC = \
-IregionSolvers \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume

View File

@ -0,0 +1,209 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
foamMultiRun
Description
Loads and executes an OpenFOAM solver modules for each region of a
multiregion simulation e.g. for conjugate heat transfer.
The region solvers are specified in the \c regionSolvers dictionary entry in
\c controlDict, containing a list of pairs of region and solver names,
e.g. for a two region case with one fluid region named
liquid and one solid region named tubeWall:
\verbatim
regionSolvers
{
liquid fluid;
tubeWall solid;
}
\endverbatim
The \c regionSolvers entry is a dictionary to support name substitutions to
simplify the specification of a single solver type for a set of
regions, e.g.
\verbatim
fluidSolver fluid;
solidSolver solid;
regionSolvers
{
tube1 $fluidSolver;
tubeWall1 solid;
tube2 $fluidSolver;
tubeWall2 solid;
tube3 $fluidSolver;
tubeWall3 solid;
}
\endverbatim
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
Usage
\b foamMultiRun [OPTION]
- \par -libs '(\"lib1.so\" ... \"libN.so\")'
Specify the additional libraries loaded
Example usage:
- To update and run a \c chtMultiRegion case add the following entries to
the controlDict:
\verbatim
application foamMultiRun;
regionSolvers
{
fluid fluid;
solid solid;
}
\endverbatim
then execute \c foamMultiRun
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "solver.H"
#include "regionSolvers.H"
#include "pimpleMultiRegionControl.H"
#include "setDeltaT.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
regionSolvers regionSolvers(runTime);
PtrList<fvMesh> regions(regionSolvers.size());
PtrList<solver> solvers(regionSolvers.size());
forAll(regionSolvers, i)
{
regions.set
(
i,
new fvMesh
(
IOobject
(
regionSolvers[i].first(),
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
solvers.set
(
i,
solver::New(regionSolvers[i].second(), regions[i])
);
}
// Create the outer PIMPLE loop and control structure
pimpleMultiRegionControl pimple(runTime, solvers);
// Set the initial time-step
setDeltaT(runTime, solvers);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
forAll(solvers, i)
{
solvers[i].preSolve();
}
// Adjust the time-step according to the solver maxDeltaT
adjustDeltaT(runTime, solvers);
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// Multi-region PIMPLE corrector loop
while (pimple.loop())
{
forAll(solvers, i)
{
solvers[i].moveMesh();
}
forAll(solvers, i)
{
solvers[i].prePredictor();
}
forAll(solvers, i)
{
solvers[i].momentumPredictor();
}
while (pimple.correctEnergy())
{
forAll(solvers, i)
{
solvers[i].thermophysicalPredictor();
}
}
forAll(solvers, i)
{
solvers[i].pressureCorrector();
}
forAll(solvers, i)
{
solvers[i].momentumTransportCorrector();
solvers[i].thermophysicalTransportCorrector();
}
}
forAll(solvers, i)
{
solvers[i].postSolve();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,85 +23,80 @@ License
\*---------------------------------------------------------------------------*/
#include "regionProperties.H"
#include "IOdictionary.H"
#include "argList.H"
#include "polyMesh.H"
#include "regionSolvers.H"
#include "Time.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::regionProperties::regionProperties(const Time& runTime)
:
HashTable<wordList>
(
IOdictionary
Foam::regionSolvers::regionSolvers(const Time& runTime)
{
if (runTime.controlDict().found("regionSolvers"))
{
const dictionary& regionSolversDict =
runTime.controlDict().subDict("regionSolvers");
forAllConstIter(dictionary, regionSolversDict, iter)
{
append(Pair<word>(iter().keyword(), iter().stream()));
}
}
else
{
// Partial backward-compatibility
// Converts the regions entry in the regionProperties dictionary into
// the regionSolvers list
typeIOobject<IOdictionary> regionPropertiesHeader
(
IOobject
(
"regionProperties",
runTime.time().constant(),
runTime.db(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
IOobject::MUST_READ
)
).lookup("regions")
)
{}
);
if (regionPropertiesHeader.headerOk())
{
HashTable<wordList> regions
(
IOdictionary(regionPropertiesHeader).lookup("regions")
);
if (regions.found("solid"))
{
const wordList& fluidRegions = regions["solid"];
forAll(fluidRegions, i)
{
append(Pair<word>(fluidRegions[i], "solid"));
}
}
if (regions.found("fluid"))
{
const wordList& fluidRegions = regions["fluid"];
forAll(fluidRegions, i)
{
append(Pair<word>(fluidRegions[i], "fluid"));
}
}
}
else
{
FatalIOErrorInFunction(runTime.controlDict())
<< "regionSolvers list missing from "
<< runTime.controlDict().name()
<< exit(FatalIOError);
}
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::regionProperties::~regionProperties()
Foam::regionSolvers::~regionSolvers()
{}
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
const Foam::word& Foam::regionDir(const word& regionName)
{
return
regionName == polyMesh::defaultRegion
? word::null
: regionName;
}
Foam::wordList Foam::selectRegionNames(const argList& args, const Time& runTime)
{
const bool allRegions = args.optionFound("allRegions");
wordList regionNames;
if (allRegions)
{
const regionProperties rp(runTime);
forAllConstIter(HashTable<wordList>, rp, iter)
{
const wordList& regions = iter();
forAll(regions, i)
{
if (findIndex(regionNames, regions[i]) == -1)
{
regionNames.append(regions[i]);
}
}
}
}
else
{
word regionName;
if (args.optionReadIfPresent("region", regionName))
{
regionNames = wordList(1, regionName);
}
else
{
regionNames = wordList(1, polyMesh::defaultRegion);
}
}
return regionNames;
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,66 +22,59 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::regionProperties
Foam::regionSolvers
Description
Simple class to hold region information for coupled region simulations.
Gives per physics ('fluid', 'solid') the names of the regions. There
is no assumption on this level that one region should only have one
set of physics.
Simple class to hold the list of region solvers
SourceFiles
regionProperties.C
regionSolvers.C
\*---------------------------------------------------------------------------*/
#ifndef regionProperties_H
#define regionProperties_H
#ifndef regionSolvers_H
#define regionSolvers_H
#include "Time.H"
#include "HashTable.H"
#include "word.H"
#include "Pair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class Time;
/*---------------------------------------------------------------------------*\
Class regionProperties Declaration
Class regionSolvers Declaration
\*---------------------------------------------------------------------------*/
class regionProperties
class regionSolvers
:
public HashTable<wordList>
public List<Pair<word>>
{
public:
// Constructors
//- Construct from components
regionProperties(const Time& runTime);
regionSolvers(const Time& runTime);
//- Disallow default bitwise copy construction
regionProperties(const regionProperties&) = delete;
regionSolvers(const regionSolvers&) = delete;
//- Destructor
~regionProperties();
~regionSolvers();
// Member Operators
//- Disallow default bitwise assignment
void operator=(const regionProperties&) = delete;
void operator=(const regionSolvers&) = delete;
};
const word& regionDir(const word& regionName);
wordList selectRegionNames(const argList& args, const Time& runTime);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -0,0 +1,93 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "setDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::setDeltaT(Time& runTime, const PtrList<solver>& solvers)
{
if
(
runTime.timeIndex() == 0
&& runTime.controlDict().lookupOrDefault("adjustTimeStep", false)
)
{
scalar deltaT = great;
forAll(solvers, i)
{
if (solvers[i].transient())
{
deltaT = min(deltaT, solvers[i].maxDeltaT());
}
}
if (deltaT != great)
{
runTime.setDeltaT(deltaT);
}
}
}
void Foam::adjustDeltaT(Time& runTime, const PtrList<solver>& solvers)
{
if (runTime.controlDict().lookupOrDefault("adjustTimeStep", false))
{
scalar deltaT = great;
forAll(solvers, i)
{
if (solvers[i].transient())
{
const scalar maxDeltaTi = solvers[i].maxDeltaT();
deltaT = min
(
deltaT,
min
(
min
(
maxDeltaTi,
runTime.deltaTValue() + 0.1*maxDeltaTi
),
1.2*runTime.deltaTValue()
)
);
}
}
if (deltaT != great)
{
runTime.setDeltaT(deltaT);
Info<< "deltaT = " << runTime.deltaTValue() << endl;
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,63 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
Foam::setDeltaT
Description
Foam::setDeltaT sets the initial time-step according to the solver maxDeltaT
Foam::adjustDeltaT adjust the time-step according to the solver maxDeltaT
SourceFiles
setDeltaT.C
\*---------------------------------------------------------------------------*/
#ifndef setDeltaT_H
#define setDeltaT_H
#include "solver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//- Set the initial time-step according to the solver maxDeltaT
void setDeltaT(Time& runTime, const PtrList<solver>& solvers);
//- Adjust the time-step according to the solver maxDeltaT
void adjustDeltaT(Time& runTime, const PtrList<solver>& solvers);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,4 @@
setDeltaT.C
foamRun.C
EXE = $(FOAM_APPBIN)/foamRun

Some files were not shown because too many files have changed in this diff Show More