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

ENH: Added new syntheic turbulence inlet boundary condition for LES/DES

Reference:
    Poletto, R., Craft, T., and Revell, A.,
    "A New Divergence Free Synthetic Eddy Method for the
    Reproduction of Inlet Flow Conditions for LES",
    Flow Turbulence Combust (2013) 91:519-539
parent c233552d
......@@ -184,6 +184,9 @@ $(derivedFvPatchFields)/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueF
$(derivedFvPatchFields)/totalPressure/totalPressureFvPatchScalarField.C
$(derivedFvPatchFields)/totalTemperature/totalTemperatureFvPatchScalarField.C
$(derivedFvPatchFields)/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C
$(derivedFvPatchFields)/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
$(derivedFvPatchFields)/turbulentDFSEMInlet/eddy/eddy.C
$(derivedFvPatchFields)/turbulentDFSEMInlet/eddy/eddyIO.C
$(derivedFvPatchFields)/turbulentInlet/turbulentInletFvPatchFields.C
$(derivedFvPatchFields)/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
$(derivedFvPatchFields)/uniformDensityHydrostaticPressure/uniformDensityHydrostaticPressureFvPatchScalarField.C
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "eddy.H"
#include "mathematicalConstants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
Foam::label Foam::eddy::Gamma2Values[] = {1, 2, 3, 4, 5, 6, 7, 8};
Foam::UList<Foam::label> Foam::eddy::Gamma2(&Gamma2Values[0], 8);
int Foam::eddy::debug = 0;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::eddy::setScales
(
const scalar sigmaX,
const label gamma2,
const vector& e,
const vector& lambda,
vector& sigma,
vector& alpha
) const
{
// Static array of gamma^2 vs c2 coefficient
static const scalar gamma2VsC2[8] =
{2, 1.875, 1.737, 1.75, 0.91, 0.825, 0.806, 1.5};
scalar gamma = Foam::sqrt(scalar(gamma2));
// c2 coefficient retrieved from array
scalar c2 = gamma2VsC2[gamma2 - 1];
// Length scale in largest eigenvalue direction
label d1 = dir1_;
label d2 = (d1 + 1) % 3;
label d3 = (d1 + 2) % 3;
sigma[d1] = sigmaX;
// Note: sigma_average = 1/3*(sigma_x + sigma_y + sigma_z)
// Substituting for sigma_y = sigma_x/gamma and sigma_z = sigma_y
//sigma[d1] = 3*sigmaX/(1 + 2/gamma);
// Other length scales equal, as function of major axis length and gamma
sigma[d2] = sigma[d1]/gamma;
sigma[d3] = sigma[d2];
vector sigma2 = cmptMultiply(sigma, sigma);
scalar slos2 = cmptSum(cmptDivide(lambda, sigma2));
bool ok = true;
for (label beta = 0; beta < 3; ++beta)
{
scalar x = slos2 - 2*lambda[beta]/sigma2[beta];
if (x < 0)
{
alpha[beta] = 0;
ok = false;
}
else
{
alpha[beta] = e[beta]*sqrt(x/(2*c2));
}
}
if (debug > 1)
{
Pout<< "c2:" << c2
<< ", gamma2:" << gamma2
<< ", gamma:" << gamma
<< ", lambda:" << lambda
<< ", sigma2: " << sigma2
<< ", slos2: " << slos2
<< ", sigmaX:" << sigmaX
<< ", sigma:" << sigma
<< ", alpha:" << alpha
<< endl;
}
return ok;
}
// * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
Foam::eddy::eddy()
:
patchFaceI_(-1),
position0_(vector::zero),
x_(0),
sigma_(vector::zero),
alpha_(vector::zero),
Rpg_(tensor::I),
c1_(-1),
dir1_(0)
{}
Foam::eddy::eddy
(
const label patchFaceI,
const point& position0,
const scalar x,
const scalar sigmaX,
const symmTensor& R,
cachedRandom& rndGen
)
:
patchFaceI_(patchFaceI),
position0_(position0),
x_(x),
sigma_(vector::zero),
alpha_(vector::zero),
Rpg_(tensor::I),
c1_(-1),
dir1_(0)
{
// Principal stresses - eigenvalues returned in ascending order
vector lambda = eigenValues(R);
// Eddy rotation from principal-to-global axes
// - given by the 3 eigenvectors of the Reynold stress tensor as rows in
// the result tensor (transposed transformation tensor)
// - returned in ascending eigenvalue order
Rpg_ = eigenVectors(R, lambda).T();
if (debug)
{
// Global->Principal transform = Rgp = Rpg.T()
// Rgp & R & Rgp.T() should have eigenvalues on its diagonal and
// zeros for all other components
Pout<< "Rpg.T() & R & Rpg: " << (Rpg_.T() & R & Rpg_) << endl;
}
// Set the eddy orientation to position of max eigenvalue
// (direction of eddy major axis, sigma_x in reference)
dir1_ = 2;
// Random vector of 1's and -1's
const vector e(epsilon(rndGen));
// Set intensities and length scales
bool found = false;
forAll(Gamma2, i)
{
// Random length scale ratio, gamma = sigmax/sigmay = sigmax/sigmaz
// - using gamma^2 to ease lookup of c2 coefficient
label g2 = Gamma2[i];
if (setScales(sigmaX, g2, e, lambda, sigma_, alpha_))
{
found = true;
break;
}
}
// Normalisation coefficient (eq. 11)
// Note: sqrt(10*V)/sqrt(nEddy) applied outside when computing uDash
c1_ = cmptAv(sigma_)/cmptProduct(sigma_)*cmptMin(sigma_);
if (found)
{
// Shuffle the gamma^2 values
rndGen.shuffle(Gamma2);
}
else
{
if (debug)
{
// If not found typically means that the stress has a repeated
// eigenvalue/not covered by the selection of Gamma values, e.g.
// as seen by range of applicability on Lumley diagram
WarningInFunction
<< "Unable to set eddy intensity for eddy: " << *this
<< endl;
}
// Remove the influence of this eddy/indicate that its initialisation
// failed
patchFaceI_ = -1;
}
}
Foam::eddy::eddy(const eddy& e)
:
patchFaceI_(e.patchFaceI_),
position0_(e.position0_),
x_(e.x_),
sigma_(e.sigma_),
alpha_(e.alpha_),
Rpg_(e.Rpg_),
c1_(e.c1_),
dir1_(e.dir1_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::vector Foam::eddy::uDash(const point& xp, const vector& n) const
{
// Relative position inside eddy (global system)
const vector r = cmptDivide(xp - position(n), sigma_);
if (mag(r) > 1)
{
return vector::zero;
}
// Relative position inside eddy (eddy principal system)
const vector rp = Rpg_.T() & r;
// Shape function (eddy principal system)
const vector q = cmptMultiply(sigma_, vector::one - cmptMultiply(rp, rp));
// Fluctuating velocity (eddy principal system) (eq. 8)
const vector uDashp = cmptMultiply(q, rp^alpha_);
// Convert into global system (eq. 10)
return c1_*(Rpg_ & uDashp);
}
void Foam::eddy::writeCentreOBJ
(
const vector& n,
Ostream& os
) const
{
point p = position(n);
os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
}
Foam::label Foam::eddy::writeSurfaceOBJ
(
const label pointOffset,
const vector& n,
Ostream& os
) const
{
if (patchFaceI_ < 0)
{
// Invalid eddy
return 0;
}
static const label nFaceAxis = 20;
static const label nFaceTheta = 22;
static const label nEddyPoints = (nFaceAxis - 1)*nFaceTheta + 2;
static FixedList<point, nEddyPoints> x;
static scalar dTheta = mathematical::twoPi/nFaceTheta;
static scalar dPhi = mathematical::pi/scalar(nFaceAxis);
label pointI = pointOffset;
const vector& s = sigma_;
const vector axisDir = tensor::I.vectorComponent(dir1_);
const label dir2 = (dir1_ + 1) % 3;
const label dir3 = (dir1_ + 2) % 3;
// Calculate the point positions
x[0] = axisDir*s[dir1_];
x[nEddyPoints - 1] = - axisDir*s[dir1_];
label eddyPtI = 1;
for (label axisI = 1; axisI < nFaceAxis; axisI++)
{
scalar z = s[dir1_]*cos(axisI*dPhi);
scalar r = sqrt(sqr(s[dir2])*(1 - sqr(z)/sqr(s[dir1_])));
for (label thetaI = 0; thetaI < nFaceTheta; thetaI++)
{
scalar theta = thetaI*dTheta;
point& p = x[eddyPtI++];
p[dir1_] = z;
p[dir2] = r*sin(theta);
p[dir3] = r*cos(theta);
}
}
// Write points
forAll(x, i)
{
point p = position(n) + (Rpg_ & x[i]);
os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
}
// Write the end cap tri faces
for (label faceI = 0; faceI < nFaceTheta; faceI++)
{
label p1 = pointI + 1;
label p2 = p1 + faceI + 1;
label p3 = p2 + 1;
if (faceI == nFaceTheta - 1) p3 -= nFaceTheta;
os << "f " << p1 << " " << p2 << " " << p3 << nl;
label q1 = pointI + nEddyPoints;
label q2 = q1 - faceI - 1;
label q3 = q2 - 1;
if (faceI == nFaceTheta - 1) q3 += nFaceTheta;
os << "f " << q1 << " " << q2 << " " << q3 << nl;
}
// Write quad faces
for (label axisI = 1; axisI < nFaceAxis - 1; axisI++)
{
for (label thetaI = 0; thetaI < nFaceTheta; thetaI++)
{
label p1 = pointI + 1 + (axisI - 1)*nFaceTheta + thetaI + 1;
label p2 = p1 + nFaceTheta;
label p3 = p2 + 1;
label p4 = p1 + 1;
if (thetaI == nFaceTheta - 1)
{
p3 -= nFaceTheta;
p4 -= nFaceTheta;
}
os << "f " << p1 << " " << p2 << " " << p3 << " " << p4 << nl;
}
}
return nEddyPoints;
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
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::eddy
Description
Class to describe eddies for the turbulentDFSEMInletFvPatchVectorField
boundary condition.
SourceFiles
eddy.C
eddyI.H
eddyIO.C
\*---------------------------------------------------------------------------*/
#ifndef turbulentDFSEMInletFvPatchVectorField_eddy_H
#define turbulentDFSEMInletFvPatchVectorField_eddy_H
#include "vector.H"
#include "point.H"
#include "tensor.H"
#include "cachedRandom.H"
#include "coordinateSystem.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class eddy Declaration
\*---------------------------------------------------------------------------*/
class eddy
{
// Private data
static label Gamma2Values[8];
static UList<label> Gamma2;
//- Patch face index that spawned the eddy
label patchFaceI_;
//- Reference position
point position0_;
//- Distance from reference position in normal direction
scalar x_;
//- Length scales in 3-D space
vector sigma_;
//- Time-averaged intensity
vector alpha_;
//- Co-ordinate system transformation from local to global axes
// X-direction aligned with max stress eigenvalue
tensor Rpg_;
//- Model coefficient c1
scalar c1_;
//- Index of streamwise direction
label dir1_;
// Private Member Functions
//- Set the eddy scales: length, intensity
bool setScales
(
const scalar sigmaX,
const label gamma2,
const vector& e,
const vector& lambda,
vector& sigma,
vector& alpha
) const;
//- Return a number with zero mean and unit variance
inline scalar epsi(cachedRandom& rndGen) const;
public:
// Constructors
//- Construct null
eddy();
//- Construct from Istream
eddy(Istream& is);
//- Construct from components
eddy
(
const label patchFaceI, // patch face index
const point& position0, // reference position
const scalar x, // distance from reference position
const scalar sigmaX, // length scale
const symmTensor& R, // Stress tensor
cachedRandom& rndGen
);
//- Construct copy
eddy(const eddy& e);
static int debug;
// Public Member Functions
// Access
//- Return the patch face index that spawned the eddy
inline label patchFaceI() const;
//- Return the reference position
inline const point& position0() const;
//- Return the distance from the reference position
inline scalar x() const;
//- Return the lLength scales in 3-D space
inline const vector& sigma() const;
//- Return the time-averaged intensity
inline const vector& alpha() const;
//- Return the co-ordinate system transformation from local
// principal to global axes
inline const tensor& Rpg() const;
//- Return the model coefficient c1
inline scalar c1() const;
//- Return the eddy position
inline point position(const vector& n) const;
//- Return the index of the streamwise direction
inline label dir1() const;
//- Return random vector of -1 and 1's
inline vector epsilon(cachedRandom& rndGen) const;
// Helper functions
//- Volume
inline scalar volume() const;
//- Move the eddy
inline void move(const scalar dx);
//- Eddy bounds
inline boundBox bounds(const bool global = true) const;
// Evaluate
//- Return the fluctuating velocity contribution at local point xp
vector uDash(const point& xp, const vector& n) const;
// Writing
//- Write the eddy centre in OBJ format
void writeCentreOBJ(const vector& n, Ostream& os) const;
//- Write the eddy surface in OBJ format
// Returns the number of points used to describe the eddy surface
label writeSurfaceOBJ
(
const label pointOffset,
const vector& n,
Ostream& os
) const;
// Member Operators
void operator=(const eddy& e);
// Friend Operators
friend bool operator==(const eddy& a, const eddy& b)
{
return
a.patchFaceI_ == b.patchFaceI_
&& a.position0_ == b.position0_
&& a.x_ == b.x_
&& a.sigma_ == b.sigma_
&& a.alpha_ == b.alpha_
&& a.Rpg_ == b.Rpg_
&& a.c1_ == b.c1_
&& a.dir1_ == b.dir1_;
}