Commit cbe5974c authored by Andrew Heather's avatar Andrew Heather
Browse files

adding tabulated U wall function BC and utility to generate tables

parent 6a62baae
#!/bin/sh
wclean libso tabulatedWallFunction
wclean
#!/bin/sh
wmake libso tabulatedWallFunction
wmake
wallFunctionTableApp.C
EXE = $(FOAM_APPBIN)/wallFunctionTable
EXE_INC = \
-I./tabulatedWallFunction/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-ltabulatedWallFunctions \
-lfiniteVolume
tabulatedWallFunction/tabulatedWallFunction.C
tabulatedWallFunction/newTabulatedWallFunction.C
SpaldingsLaw/SpaldingsLaw.C
general/general.C
LIB = $(FOAM_LIBBIN)/libtabulatedWallFunctions
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
-I$(LIB_SRC)/transportModels
LIB_LIBS = \
-lfiniteVolume
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-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
\*---------------------------------------------------------------------------*/
#include "SpaldingsLaw.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace tabulatedWallFunctions
{
defineTypeNameAndDebug(SpaldingsLaw, 0);
addToRunTimeSelectionTable
(
tabulatedWallFunction,
SpaldingsLaw,
dictionary
);
}
}
const Foam::label Foam::tabulatedWallFunctions::SpaldingsLaw::maxIters_ = 1000;
const Foam::scalar
Foam::tabulatedWallFunctions::SpaldingsLaw::tolerance_ = 1e-4;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::tabulatedWallFunctions::SpaldingsLaw::invertFunction()
{
// Initialise u+ and R
scalar Re = 0.0;
scalar uPlus = 1;
// Populate the table
for (label i=0; i<invertedTable_.size(); i++)
{
if (invertedTable_.log10())
{
Re = pow(10, (i*invertedTable_.dx() + invertedTable_.x0()));
}
else
{
Re = i*invertedTable_.dx() + invertedTable_.x0();
}
// Use latest available u+ estimate
if (i > 0)
{
uPlus = invertedTable_[i-1];
}
// Newton iterations to determine u+
label iter = 0;
scalar error = GREAT;
do
{
scalar kUPlus = min(kappa_*uPlus, 50);
scalar A =
E_*sqr(uPlus)
+ uPlus
*(exp(kUPlus) - pow3(kUPlus)/6 - 0.5*sqr(kUPlus) - kUPlus - 1);
scalar f = - Re + A/E_;
scalar df =
(
2*E_*uPlus
+ exp(kUPlus)*(kUPlus + 1)
- 2/3*pow3(kUPlus)
- 1.5*sqr(kUPlus)
- 2*kUPlus
- 1
)/E_;
scalar uPlusNew = uPlus - f/(df + ROOTVSMALL);
error = mag((uPlus - uPlusNew)/uPlusNew);
uPlus = uPlusNew;
} while (error > tolerance_ && ++iter < maxIters_);
if (iter == maxIters_)
{
WarningIn("SpaldingsLaw::invertFunction()")
<< "Newton iterations not converged:" << nl
<< " iters = " << iter << ", error = " << error << endl;
}
// Set new values - constrain u+ >= 0
invertedTable_[i] = max(0, uPlus);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::tabulatedWallFunctions::SpaldingsLaw::SpaldingsLaw
(
const dictionary& dict,
const polyMesh& mesh
)
:
tabulatedWallFunction(dict, mesh, typeName),
kappa_(readScalar(coeffDict_.lookup("kappa"))),
E_(readScalar(coeffDict_.lookup("E")))
{
invertFunction();
if (debug)
{
writeData(Info);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::tabulatedWallFunctions::SpaldingsLaw::~SpaldingsLaw()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
Foam::scalar Foam::tabulatedWallFunctions::SpaldingsLaw::yPlus
(
const scalar uPlus
) const
{
scalar kUPlus = min(kappa_*uPlus, 50);
return
uPlus
+ 1/E_*(exp(kUPlus) - pow3(kUPlus)/6 - 0.5*sqr(kUPlus) - kUPlus - 1);
}
Foam::scalar Foam::tabulatedWallFunctions::SpaldingsLaw::Re
(
const scalar uPlus
) const
{
return uPlus*yPlus(uPlus);
}
void Foam::tabulatedWallFunctions::SpaldingsLaw::writeData(Ostream& os) const
{
if (invertedTable_.log10())
{
os << "log10(Re), y+, u+:" << endl;
forAll(invertedTable_, i)
{
scalar uPlus = invertedTable_[i];
scalar Re = ::log10(this->Re(uPlus));
scalar yPlus = this->yPlus(uPlus);
os << Re << ", " << yPlus << ", " << uPlus << endl;
}
}
else
{
os << "Re, y+, u+:" << endl;
forAll(invertedTable_, i)
{
scalar uPlus = invertedTable_[i];
scalar Re = this->Re(uPlus);
scalar yPlus = this->yPlus(uPlus);
os << Re << ", " << yPlus << ", " << uPlus << endl;
}
}
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-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
Class
Foam::tabulatedWallFunctions::SpaldingsLaw
Description
Computes U+ as a function of Reynolds number by inverting Spaldings law.
Example dictionary specification:
tabulatedWallFunction SpaldingsLaw;
// Output table info
tableName uPlusWallFunctionData; // Output table name
log10 yes; // Rey interpreted as log10(Rey)
dx 0.2; // Interval log10(Rey)
x0 -3; // Minimum log10(Rey)
xMax 7; // Maximum log10(Rey)
SpaldingsLawCoeffs
{
kappa 0.41; // Von Karman constant
E 9.8; // Law-of-the-wall E coefficient
}
SourceFiles
SpaldingsLaw.C
\*---------------------------------------------------------------------------*/
#ifndef SpaldingsLaw_H
#define SpaldingsLaw_H
#include "tabulatedWallFunction.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace tabulatedWallFunctions
{
/*---------------------------------------------------------------------------*\
Class SpaldingsLaw Declaration
\*---------------------------------------------------------------------------*/
class SpaldingsLaw
:
public tabulatedWallFunction
{
protected:
// Protected data
//- Von Karman constant
scalar kappa_;
//- Law-of-the-wall E coefficient
scalar E_;
// Newton iteration solution properties
//- Maximum number of iterations
static const label maxIters_;
//- Tolerance
static const scalar tolerance_;
// Protected member functions
//- Invert the function
virtual void invertFunction();
public:
//- Run-time type information
TypeName("SpaldingsLaw");
// Constructors
SpaldingsLaw(const dictionary& dict, const polyMesh& mesh);
//- Destructor
virtual ~SpaldingsLaw();
// Member Functions
// Access
//- Return y+ as a function of u+
virtual scalar yPlus(const scalar uPlus) const;
//- Return Reynolds number as a function of u+
virtual scalar Re(const scalar uPlus) const;
// I-O
//- Write to Ostream
virtual void writeData(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace tabulatedWallFunctions
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-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
\*---------------------------------------------------------------------------*/
#include "general.H"
#include "addToRunTimeSelectionTable.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace tabulatedWallFunctions
{
defineTypeNameAndDebug(general, 0);
addToRunTimeSelectionTable
(
tabulatedWallFunction,
general,
dictionary
);
}
}
template<>
const char*
Foam::NamedEnum<Foam::tabulatedWallFunctions::general::interpolationType, 1>::
names[] =
{
"linear"
};
const
Foam::NamedEnum<Foam::tabulatedWallFunctions::general::interpolationType, 1>
Foam::tabulatedWallFunctions::general::interpolationTypeNames_;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::tabulatedWallFunctions::general::invertTable()
{
scalarList Rey(uPlus_.size(), 0.0);
// Calculate Reynolds number
forAll(uPlus_, i)
{
Rey[i] = yPlus_[i]*uPlus_[i];
if (invertedTable_.log10())
{
Rey[i] = ::log10(max(ROOTVSMALL, Rey[i]));
}
}
// Populate the U+ vs Re table
forAll(invertedTable_, i)
{
scalar Re = i*invertedTable_.dx() + invertedTable_.x0();
invertedTable_[i] = max(0, interpolate(Re, Rey, uPlus_));
}
}
Foam::scalar Foam::tabulatedWallFunctions::general::interpolate
(
const scalar xi,
const scalarList& x,
const scalarList& fx
) const
{
switch (interpType_)
{
case itLinear:
{
if (xi < x[0])
{
return fx[0];
}
else if (xi > x[x.size()-1])
{
return fx[x.size()-1];
}
else
{
label i2 = 0;
while (x[i2] < xi)
{
i2++;
}
label i1 = i2 - 1;
return (xi - x[i1])/(x[i2] - x[i1])*(fx[i2] - fx[i1]) + fx[i1];
}
break;
}
default:
{
FatalErrorIn
(
"tabulatedWallFunctions::general::interpolate"
"("
"const scalar, "
"const scalarList&, "
"const scalarList&"
")"
) << "Unknown interpolation method" << nl
<< abort(FatalError);
}
}
return 0.0;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::tabulatedWallFunctions::general::general
(
const dictionary& dict,
const polyMesh& mesh
)
:
tabulatedWallFunction(dict, mesh, typeName),
interpType_(interpolationTypeNames_[coeffDict_.lookup("interpType")]),
yPlus_(),
uPlus_(),
log10YPlus_(coeffDict_.lookup("log10YPlus")),
log10UPlus_(coeffDict_.lookup("log10UPlus"))
{
List<Tuple2<scalar, scalar> > inputTable = coeffDict_.lookup("inputTable");
if (inputTable.size() < 2)
{
FatalErrorIn
(
"tabulatedWallFunctions::general::general"
"("
"const dictionary&, "
"const polyMesh&"
")"
) << "Input table must have at least 2 values" << nl
<< exit(FatalError);
}
yPlus_.setSize(inputTable.size());
uPlus_.setSize(inputTable.size());
forAll(inputTable, i)
{
if (log10YPlus_)
{
yPlus_[i] = pow(10, inputTable[i].first());
}
else
{
yPlus_[i] = inputTable[i].first();
}
if (log10UPlus_)
{
uPlus_[i] = pow(10, inputTable[i].second());
}
else
{
uPlus_[i] = inputTable[i].second();
}
}