Commit 5ce9333c authored by mattijs's avatar mattijs
Browse files

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

parents 0685f59b 68965204
buoyantBoussinesqFoam.C
EXE = $(FOAM_APPBIN)/buoyantBoussinesqFoam
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lincompressibleRASModels \
-lincompressibleTransportModels
{
volScalarField kappaEff
(
"kappaEff",
turbulence->nu() + turbulence->nut()/Prt
);
fvScalarMatrix TEqn
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(kappaEff, T)
);
TEqn.relax();
TEqn.solve();
}
// Solve the momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
UEqn().relax();
solve
(
UEqn()
==
-fvc::reconstruct
(
(
fvc::snGrad(pd)
- betaghf*fvc::snGrad(T)
) * mesh.magSf()
)
);
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
buoyantBoussinesqSimpleFoam
Description
Steady-state solver for buoyant, turbulent flow of incompressible fluids
Uses the Boussinesq approximation:
\f[
rho_{eff} = 1 - beta(T - T_{ref})
\f]
where:
\f$ rho_{eff} \f$ = the effective (driving) density
beta = thermal expansion coefficient [1/K]
T = temperature [K]
\f$ T_{ref} \f$ = reference temperature [K]
Valid when:
\f[
rho_{eff} << 1
\f]
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "createFields.H"
# include "initContinuityErrs.H"
# include "readTimeControls.H"
# include "CourantNo.H"
# include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readTimeControls.H"
# include "readPISOControls.H"
# include "CourantNo.H"
# include "setDeltaT.H"
# include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
# include "TEqn.H"
# include "pdEqn.H"
}
turbulence->correct();
if (runTime.write())
{
# include "writeAdditionalFields.H"
}
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<basicThermo> thermo
(
basicThermo::New(mesh)
);
volScalarField rho
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"rho",
"T",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
thermo->rho()
mesh
);
volScalarField& p = thermo->p();
volScalarField& h = thermo->h();
const volScalarField& T = thermo->T();
Info<< "Reading field U\n" << endl;
volVectorField U
// kinematic pd
Info<< "Reading field pd\n" << endl;
volScalarField pd
(
IOobject
(
"U",
"pd",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
......@@ -37,33 +29,12 @@
mesh
);
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::RASModel> turbulence
(
compressible::RASModel::New
(
rho,
U,
phi,
thermo()
)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));
Info<< "Creating field pd\n" << endl;
volScalarField pd
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"pd",
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
......@@ -72,9 +43,18 @@
mesh
);
p = pd + rho*gh + pRef;
thermo->correct();
# include "createPhi.H"
# include "readTransportProperties.H"
Info<< "Creating turbulence model\n" << endl;
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
Info<< "Calculating field beta*(g.h)\n" << endl;
surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
label pdRefCell = 0;
scalar pdRefValue = 0.0;
......@@ -85,13 +65,3 @@
pdRefCell,
pdRefValue
);
Info<< "Creating radiation model\n" << endl;
autoPtr<radiation::radiationModel> radiation
(
radiation::radiationModel::New(T)
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
{
volScalarField rUA("rUA", 1.0/UEqn().A());
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
U = rUA*UEqn().H();
UEqn.clear();
surfaceScalarField phiU
(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi)
);
phi = phiU + betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rUAf, pd) == fvc::div(phi)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi += pdEqn.flux();
}
}
U -= rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
#include "continuityErrs.H"
}
singlePhaseTransportModel laminarTransport(U, phi);
// thermal expansion coefficient [1/K]
dimensionedScalar beta(laminarTransport.lookup("beta"));
// reference temperature [K]
dimensionedScalar TRef(laminarTransport.lookup("TRef"));
// reference kinematic pressure [m2/s2]
dimensionedScalar pRef(laminarTransport.lookup("pRef"));
// turbulent Prandtl number
dimensionedScalar Prt(laminarTransport.lookup("Prt"));
{
volScalarField rhoEff
(
IOobject
(
"rhoEff",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
1.0 - beta*(T - TRef)
);
rhoEff.write();
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pd + rhoEff*(g & mesh.C()) + pRef
);
p.write();
}
......@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "createFields.H"
# include "createRadiationModel.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
// check convergence
if (maxResidual < convergenceCriterion)
{
Info<< "reached convergence criterion: " << convergenceCriterion << endl;
runTime.writeAndEnd();
Info<< "latestTime = " << runTime.timeName() << endl;
}
Info<< "Creating radiation model\n" << endl;
autoPtr<radiation::radiationModel> radiation
(
radiation::radiationModel::New(thermo->T())
);
// initialize values for convergence checks
scalar eqnResidual = 1, maxResidual = 0;
scalar convergenceCriterion = 0;
simple.readIfPresent("convergence", convergenceCriterion);
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......
......@@ -12,6 +12,6 @@
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T())
<< endl;
Info<< "Min/max T:" << min(thermo.T()).value() << ' '
<< max(thermo.T()).value() << endl;
}
......@@ -52,27 +52,47 @@ int main(int argc, char *argv[])
list1 = -1;
list1.print(Info);
Info<< "\ntest zero\n";
list1 = 0;
list1.print(Info);
Info<< "\ntest set() with default argument (max_value)\n";
list1.set(3);
list1.print(Info);
Info<< "\ntest assign between references\n";
list1[2] = 3;
list1[4] = list1[2];
list1.print(Info);
Info<< "\ntest assign between references, with chaining\n";
list1[4] = list1[2] = 1;
list1[0] = list1[4] = 1;
list1.print(Info);
Info<< "\ntest assign between references, with chaining and auto-vivify\n";
list1[1] = list1[8] = list1[10] = list1[14] = 2;
list1.print(Info);
{
const PackedList<3>& constLst = list1;
Info<< "\ntest operator[] const with out-of-range index\n";
constLst.print(Info);
if (!constLst[20])
if (constLst[20])
{
Info<< "[20] is true (unexpected)\n";
}
else
{
Info<< "[20] is false (expected) list size should be unchanged (const)\n";
}
constLst.print(Info);
Info<< "\ntest operator[] non-const with out-of-range index\n";
if (!list1[20])
if (list1[20])
{
Info<< "[20] is true (unexpected)\n";
}
else
{
Info<< "[20] is false (expected) but list was resized?? (non-const)\n";
}
......@@ -91,6 +111,14 @@ int main(int argc, char *argv[])
list1.resize(8, list1.max_value());
list1.print(Info);
Info<< "\ntest flip() function\n";
list1.flip();
list1.print(Info);
Info<< "\nre-flip()\n";
list1.flip();
list1.print(Info);
Info<< "\ntest set() function\n";
list1.set(1, 5);
list1.print(Info);
......@@ -188,15 +216,23 @@ int main(int argc, char *argv[])
{
Info<< "\ntest assignment of iterator\n";
list1.print(Info);
PackedList<3>::iterator cit = list1[25];
cit.print(Info);
Info<< "cend()\n";
list1.end().print(Info);
PackedList<3>::iterator cit = list1[100];
Info<< "out-of-range: ";
cit.print(Info);
cit = list1[15];
Info<< "in-range: ";
cit.print(Info);
Info<< "out-of-range: ";
cit = list1[1000];
cit.print(Info);
}
for
(
PackedList<3>::iterator cit = list1[5];
PackedList<3>::iterator cit = list1[30];
cit != list1.end();
++cit
)
......@@ -204,14 +240,19 @@ int main(int argc, char *argv[])
cit.print(Info);
}
// Info<< "\ntest operator[] auto-vivify\n";
// const unsigned int val = list1[45];
//
// Info<< "list[45]:" << val << "\n";
// list1[45] = list1.max_value();
// Info<< "list[45]:" << list1[45] << "\n";
// list1[49] = list1.max_value();
// list1.print(Info);
Info<< "\ntest operator[] auto-vivify\n";
Info<< "size:" << list1.size() << "\n";
const unsigned int val = list1[45];
Info<< "list[45]:" << val << "\n";
Info<< "size after read:" << list1.size() << "\n";
list1[45] = list1.max_value();
Info<< "size after write:" << list1.size() << "\n";
Info<< "list[45]:" << list1[45] << "\n";
list1[49] = list1[100];
list1.print(Info);
Info<< "\ntest copy constructor + append\n";
......@@ -235,6 +276,23 @@ int main(int argc, char *argv[])
Info<< "removed final value: " << list3.remove() << endl;
list3.print(Info);
List<bool> list4(4, true);
{
const List<bool>& constLst = list4;
Info<< "\ntest operator[] const with out-of-range index\n";
Info<< constLst << endl;
if (constLst[20])
{
Info<< "[20] is true (unexpected)\n";
}
else
{
Info<< "[20] is false (expected) list size should be unchanged (const)\n";
}
Info<< constLst << endl;
}
Info<< "\n\nDone.\n";
return 0;
......
......@@ -32,6 +32,7 @@ Description
#include "fileName.H"
#include "SubList.H"
#include "IOobject.H"
#include "IOstreams.H"
#include "OSspecific.H"
......@@ -61,7 +62,8 @@ int main()
<< endl;