zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: turbulentDFSEMInlet: refactoring by PatchFunction1

Browse Source 
- ENH: turbulentDFSEMInlet: add normalisation factors for
input Reynolds stresses, mean velocity and integral-length
scales as entries `Uref` and `Lref`.
- ENH: turbulentDFSEMInlet: add scaling factor entries, `scale`
and `m`, to enable users to tune C1 normalisation coefficient,
if need be.
- BUG: turbulentDFSEM: (fixes #1004 #1744 #2089)
- see #2090 for theoretical issues related to the DFSEM method.
This commit is contained in:
Kutalmis Bercin
2021-05-11 10:23:01 +01:00
committed by Andrew Heather
parent c6759692ba
commit b9c174312b
6 changed files with 450 additions and 678 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2020 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -49,37 +49,37 @@ bool Foam::eddy::setScales
vector& alpha
) const
{
// Static array of gamma^2 vs c2 coefficient
// Static array of gamma^2 vs c2 coefficient (PCR:Table 1)
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));
const scalar gamma = Foam::sqrt(scalar(gamma2));
// c2 coefficient retrieved from array
scalar c2 = gamma2VsC2[gamma2 - 1];
const scalar c2 = gamma2VsC2[gamma2 - 1];
// Length scale in largest eigenvalue direction
label d1 = dir1_;
label d2 = (d1 + 1) % 3;
label d3 = (d1 + 2) % 3;
// Length scale in the largest eigenvalue direction
const label d1 = dir1_;
const label d2 = (d1 + 1) % 3;
const 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));
// (PCR:Eq. 13)
const vector sigma2(cmptMultiply(sigma, sigma));
const scalar slos2 = cmptSum(cmptDivide(lambda, sigma2));
bool ok = true;
for (label beta = 0; beta < 3; ++beta)
for (label beta = 0; beta < vector::nComponents; ++beta)
{
scalar x = slos2 - 2*lambda[beta]/sigma2[beta];
const scalar x = slos2 - 2*lambda[beta]/sigma2[beta];
if (x < 0)
{
@ -88,6 +88,7 @@ bool Foam::eddy::setScales
}
else
{
// (SST:Eq. 23)
alpha[beta] = e[beta]*sqrt(x/(2*c2));
}
}
@ -145,7 +146,7 @@ Foam::eddy::eddy
dir1_(0)
{
// Principal stresses - eigenvalues returned in ascending order
vector lambda(eigenValues(R));
const vector lambda(eigenValues(R));
// Eddy rotation from principal-to-global axes
// - given by the 3 eigenvectors of the Reynold stress tensor as rows in
@ -174,17 +175,17 @@ Foam::eddy::eddy
{
// Random length scale ratio, gamma = sigmax/sigmay = sigmax/sigmaz
// - using gamma^2 to ease lookup of c2 coefficient
label g2 = Gamma2[i];
const label gamma2 = Gamma2[i];
if (setScales(sigmaX, g2, e, lambda, sigma_, alpha_))
if (setScales(sigmaX, gamma2, e, lambda, sigma_, alpha_))
{
found = true;
break;
}
}
// Normalisation coefficient (eq. 11)
// Note: sqrt(10*V)/sqrt(nEddy) applied outside when computing uDash
// Normalisation coefficient (PCR:Eq. 11)
// Note: sqrt(10*V)/sqrt(nEddy) applied outside when computing uPrime
c1_ = cmptAv(sigma_)/cmptProduct(sigma_)*cmptMin(sigma_);
if (found)
@ -226,27 +227,27 @@ Foam::eddy::eddy(const eddy& e)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::vector Foam::eddy::uDash(const point& xp, const vector& n) const
Foam::vector Foam::eddy::uPrime(const point& xp, const vector& n) const
{
// Relative position inside eddy (global system)
const vector r = cmptDivide(xp - position(n), sigma_);
// Relative position inside eddy (global system) (PCR:p. 524)
const vector r(cmptDivide(xp - position(n), sigma_));
if (mag(r) > 1)
if (mag(r) >= scalar(1))
{
return vector::zero;
}
// Relative position inside eddy (eddy principal system)
const vector rp = Rpg_.T() & r;
const vector rp(Rpg_.T() & r);
// Shape function (eddy principal system)
const vector q = cmptMultiply(sigma_, vector::one - cmptMultiply(rp, rp));
const vector q(cmptMultiply(sigma_, vector::one - cmptMultiply(rp, rp)));
// Fluctuating velocity (eddy principal system) (eq. 8)
const vector uDashp = cmptMultiply(q, rp^alpha_);
// Fluctuating velocity (eddy principal system) (PCR:Eq. 8)
const vector uPrimep(cmptMultiply(q, rp^alpha_));
// Convert into global system (eq. 10)
return c1_*(Rpg_ & uDashp);
// Convert into global system (PCR:Eq. 10)
return c1_*(Rpg_ & uPrimep);
}
@ -256,7 +257,7 @@ void Foam::eddy::writeCentreOBJ
Ostream& os
) const
{
point p = position(n);
const point p(position(n));
os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
}
@ -295,14 +296,14 @@ Foam::label Foam::eddy::writeSurfaceOBJ
x[nEddyPoints - 1] = - axisDir*s[dir1_];
label eddyPtI = 1;
for (label axisI = 1; axisI < nFaceAxis; axisI++)
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_])));
const scalar z = s[dir1_]*cos(axisI*dPhi);
const scalar r = sqrt(sqr(s[dir2])*(1 - sqr(z)/sqr(s[dir1_])));
for (label thetaI = 0; thetaI < nFaceTheta; thetaI++)
for (label thetaI = 0; thetaI < nFaceTheta; ++thetaI)
{
scalar theta = thetaI*dTheta;
const scalar theta = thetaI*dTheta;
point& p = x[eddyPtI++];
p[dir1_] = z;
p[dir2] = r*sin(theta);
@ -313,33 +314,33 @@ Foam::label Foam::eddy::writeSurfaceOBJ
// Write points
forAll(x, i)
{
point p = position(n) + (Rpg_ & x[i]);
const 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++)
for (label faceI = 0; faceI < nFaceTheta; ++faceI)
{
label p1 = pointI + 1;
label p2 = p1 + faceI + 1;
const label p1 = pointI + 1;
const 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;
const label q1 = pointI + nEddyPoints;
const 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 axisI = 1; axisI < nFaceAxis - 1; ++axisI)
{
for (label thetaI = 0; thetaI < nFaceTheta; thetaI++)
for (label thetaI = 0; thetaI < nFaceTheta; ++thetaI)
{
label p1 = pointI + 1 + (axisI - 1)*nFaceTheta + thetaI + 1;
label p2 = p1 + nFaceTheta;
const label p1 = pointI + 1 + (axisI - 1)*nFaceTheta + thetaI + 1;
const label p2 = p1 + nFaceTheta;
label p3 = p2 + 1;
label p4 = p1 + 1;

38
src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddy.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -69,7 +69,7 @@ Ostream& operator<<(Ostream& os, const eddy& e);
class eddy
{
// Private data
// Private Data
static label Gamma2Values[8];
static UList<label> Gamma2;
@ -83,7 +83,7 @@ class eddy
//- Distance from reference position in normal direction
scalar x_;
//- Length scales in 3-D space
//- Integral-length scales in 3-D space
vector sigma_;
//- Time-averaged intensity
@ -133,14 +133,18 @@ public:
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
const scalar sigmaX, // integral-length scale
const symmTensor& R, // Reynolds stress tensor
Random& rndGen
);
//- Construct copy
//- Copy construct
eddy(const eddy& e);
// Public Data
//- Flag to activate debug statements
static int debug;
@ -149,26 +153,26 @@ public:
// Access
//- Return the patch face index that spawned the eddy
inline label patchFaceI() const;
inline label patchFaceI() const noexcept;
//- Return the reference position
inline const point& position0() const;
inline const point& position0() const noexcept;
//- Return the distance from the reference position
inline scalar x() const;
inline scalar x() const noexcept;
//- Return the lLength scales in 3-D space
inline const vector& sigma() const;
//- Return the length scales in 3-D space
inline const vector& sigma() const noexcept;
//- Return the time-averaged intensity
inline const vector& alpha() const;
inline const vector& alpha() const noexcept;
//- Return the coordinate system transformation from local
// principal to global axes
inline const tensor& Rpg() const;
//- principal to global axes
inline const tensor& Rpg() const noexcept;
//- Return the model coefficient c1
inline scalar c1() const;
inline scalar c1() const noexcept;
//- Return the eddy position
inline point position(const vector& n) const;
@ -192,10 +196,10 @@ public:
inline boundBox bounds(const bool global = true) const;
// Evaluate
// Evaluation
//- Return the fluctuating velocity contribution at local point xp
vector uDash(const point& xp, const vector& n) const;
vector uPrime(const point& xp, const vector& n) const;
// Writing

26
src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/eddy/eddyI.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -39,54 +39,54 @@ Foam::scalar Foam::eddy::epsi(Random& rndGen) const
}
inline Foam::label Foam::eddy::patchFaceI() const
inline Foam::label Foam::eddy::patchFaceI() const noexcept
{
return patchFaceI_;
}
inline const Foam::point& Foam::eddy::position0() const
inline const Foam::point& Foam::eddy::position0() const noexcept
{
return position0_;
}
inline Foam::scalar Foam::eddy::x() const
inline Foam::scalar Foam::eddy::x() const noexcept
{
return x_;
}
inline const Foam::vector& Foam::eddy::sigma() const
inline const Foam::vector& Foam::eddy::sigma() const noexcept
{
return sigma_;
}
inline const Foam::vector& Foam::eddy::alpha() const
inline const Foam::vector& Foam::eddy::alpha() const noexcept
{
return alpha_;
}
inline const Foam::tensor& Foam::eddy::Rpg() const
inline const Foam::tensor& Foam::eddy::Rpg() const noexcept
{
return Rpg_;
}
inline Foam::scalar Foam::eddy::c1() const noexcept
{
return c1_;
}
inline Foam::point Foam::eddy::position(const vector& n) const
{
return position0_ + n*x_;
}
inline Foam::scalar Foam::eddy::c1() const
{
return c1_;
}
Foam::vector Foam::eddy::epsilon(Random& rndGen) const
{
return vector(epsi(rndGen), epsi(rndGen), epsi(rndGen));

570
src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2020 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,13 +27,9 @@ License
\*---------------------------------------------------------------------------*/
#include "turbulentDFSEMInletFvPatchVectorField.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "momentOfInertia.H"
#include "OFstream.H"
#include "globalIndex.H"
#include "rawIOField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -51,7 +47,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ() const
const labelList& boundaryPoints = pp.boundaryPoints();
const pointField& localPoints = pp.localPoints();
vector offset = patchNormal_*maxSigmaX_;
const vector offset(patchNormal_*maxSigmaX_);
forAll(boundaryPoints, i)
{
point p = localPoints[boundaryPoints[i]];
@ -65,24 +61,6 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ() const
p -= offset;
os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
}
// Draw lines between points
// Note: need to order to avoid crossing patch
//const label nPoint = boundaryPoints.size();
//
//forAll(boundaryPoints, i)
//{
// label i1 = i;
// label i2 = (i + 1) % nPoint;
// os << "l " << i1 << " " << i2 << nl;
//}
//
//forAll(boundaryPoints, i)
//{
// label i1 = i + nPoint;
// label i2 = ((i + 1) % nPoint) + nPoint;
// os << "l " << i1 << " " << i2 << nl;
//}
}
{
@ -106,141 +84,31 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::writeLumleyCoeffs() const
{
// Output list of xi vs eta
// Before interpolation/raw data
if (interpolateR_)
OFstream os(db().time().path()/"lumley_interpolated.out");
os << "# xi" << token::TAB << "eta" << endl;
const scalar t = db().time().timeOutputValue();
const symmTensorField R(R_->value(t)/sqr(Uref_));
forAll(R, faceI)
{
const fileName valsFile
(
fileName
(
this->db().time().globalPath()
/this->db().time().constant()
/"boundaryData"
/this->patch().name()
/"0"
/"R"
)
);
// Normalised anisotropy tensor
const symmTensor devR(dev(R[faceI]/(tr(R[faceI]))));
IOobject io
(
valsFile, // absolute path
this->db().time(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false, // no need to register
true // is global object (currently not used)
);
// Second tensor invariant
const scalar ii = min(0, invariantII(devR));
const rawIOField<symmTensor> Rexp(io, false);
// Third tensor invariant
const scalar iii = invariantIII(devR);
OFstream os(db().time().path()/"lumley_input.out");
os << "# xi" << token::TAB << "eta" << endl;
forAll(Rexp, faceI)
{
// Normalised anisotropy tensor
symmTensor devR = dev(Rexp[faceI]/(tr(Rexp[faceI])));
// Second tensor invariant
scalar ii = min(0, invariantII(devR));
// Third tensor invariant
scalar iii = invariantIII(devR);
// xi, eta
// See Pope - characterization of Reynolds-stress anisotropy
scalar xi = cbrt(0.5*iii);
scalar eta = sqrt(-ii/3.0);
os << xi << token::TAB << eta << token::TAB
<< ii << token::TAB << iii << endl;
}
// xi, eta
// See Pope - characterization of Reynolds-stress anisotropy
const scalar xi = cbrt(0.5*iii);
const scalar eta = sqrt(-ii/3.0);
os << xi << token::TAB << eta << token::TAB
<< ii << token::TAB << iii << endl;
}
// After interpolation
{
OFstream os(db().time().path()/"lumley_interpolated.out");
os << "# xi" << token::TAB << "eta" << endl;
forAll(R_, faceI)
{
// Normalised anisotropy tensor
symmTensor devR = dev(R_[faceI]/(tr(R_[faceI])));
// Second tensor invariant
scalar ii = min(0, invariantII(devR));
// Third tensor invariant
scalar iii = invariantIII(devR);
// xi, eta
// See Pope - characterization of Reynolds-stress anisotropy
scalar xi = cbrt(0.5*iii);
scalar eta = sqrt(-ii/3.0);
os << xi << token::TAB << eta << token::TAB
<< ii << token::TAB << iii << endl;
}
}
}
const Foam::pointToPointPlanarInterpolation&
Foam::turbulentDFSEMInletFvPatchVectorField::patchMapper() const
{
// Initialise interpolation (2D planar interpolation by triangulation)
if (!mapperPtr_)
{
const fileName samplePointsFile
(
this->db().time().globalPath()
/this->db().time().constant()
/"boundaryData"
/this->patch().name()
/"points"
);
IOobject io
(
samplePointsFile, // absolute path
this->db().time(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false, // no need to register
true // is global object (currently not used)
);
// Read data
const rawIOField<point> samplePoints(io, false);
DebugInFunction
<< " Read " << samplePoints.size() << " sample points from "
<< samplePointsFile << endl;
// tbd: run-time selection
bool nearestOnly =
(
!mapMethod_.empty()
&& mapMethod_ != "planarInterpolation"
);
// Allocate the interpolator
mapperPtr_.reset
(
new pointToPointPlanarInterpolation
(
samplePoints,
this->patch().patch().faceCentres(),
perturb_,
nearestOnly
)
);
}
return *mapperPtr_;
}
@ -252,7 +120,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
patchNormal_ = -gAverage(nf);
// Check that patch is planar
scalar error = max(magSqr(patchNormal_ + nf));
const scalar error = max(magSqr(patchNormal_ + nf));
if (error > SMALL)
{
@ -293,9 +161,9 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
}
}
forAll(sumTriMagSf_, i)
for (auto& s : sumTriMagSf_)
{
sumTriMagSf_[i] = 0.0;
s = 0.0;
}
sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
@ -303,7 +171,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
Pstream::listCombineScatter(sumTriMagSf_);
for (label i = 1; i < triMagSf.size(); i++)
for (label i = 1; i < triMagSf.size(); ++i)
{
triMagSf[i] += triMagSf[i-1];
}
@ -314,7 +182,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
triCumulativeMagSf_.transfer(triMagSf);
// Convert sumTriMagSf_ into cumulative sum of areas per proc
for (label i = 1; i < sumTriMagSf_.size(); i++)
for (label i = 1; i < sumTriMagSf_.size(); ++i)
{
sumTriMagSf_[i] += sumTriMagSf_[i-1];
}
@ -336,7 +204,9 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
{
const scalarField& magSf = patch().magSf();
//const scalarField cellDx(Foam::sqrt(magSf));
const scalarField L(L_->value(db().time().timeOutputValue())/Lref_);
// (PCF:Eq. 14)
const scalarField cellDx(max(Foam::sqrt(magSf), 2/patch().deltaCoeffs()));
// Inialise eddy box extents
@ -344,13 +214,10 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
{
scalar& s = sigmax_[faceI];
// Length scale in x direction (based on eq. 14)
s = mag(L_[faceI]);
s = min(s, kappa_*delta_);
// Allow eddies to be smaller than the mesh scale as suggested by
// the reference?
// s = min(s, nCellPerEddy_*cellDx[faceI]);
// Average length scale (SST:Eq. 24)
// Personal communication regarding (PCR:Eq. 14)
// - the min operator in Eq. 14 is a typo, and should be a max operator
s = min(mag(L[faceI]), kappa_*delta_);
s = max(s, nCellPerEddy_*cellDx[faceI]);
}
@ -383,7 +250,8 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
if (global)
{
scalar areaFraction = rndGen_.globalPosition<scalar>(0, patchArea_);
const scalar areaFraction =
rndGen_.globalPosition<scalar>(0, patchArea_);
// Determine which processor to use
label procI = 0;
@ -400,7 +268,7 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
{
// Find corresponding decomposed face triangle
label triI = 0;
scalar offset = sumTriMagSf_[procI];
const scalar offset = sumTriMagSf_[procI];
forAllReverse(triCumulativeMagSf_, i)
{
if (areaFraction > triCumulativeMagSf_[i] + offset)
@ -423,8 +291,8 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
{
// Find corresponding decomposed face triangle on local processor
label triI = 0;
scalar maxAreaLimit = triCumulativeMagSf_.last();
scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
const scalar maxAreaLimit = triCumulativeMagSf_.last();
const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
forAllReverse(triCumulativeMagSf_, i)
{
@ -450,6 +318,9 @@ Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
{
const scalar t = db().time().timeOutputValue();
const symmTensorField R(R_->value(t)/sqr(Uref_));
DynamicList<eddy> eddies(size());
// Initialise eddy properties
@ -465,18 +336,18 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
{
// Get new parallel consistent position
pointIndexHit pos(setNewPosition(true));
label faceI = pos.index();
const label patchFaceI = pos.index();
// Note: only 1 processor will pick up this face
if (faceI != -1)
if (patchFaceI != -1)
{
eddy e
(
faceI,
patchFaceI,
pos.hitPoint(),
rndGen_.position<scalar>(-maxSigmaX_, maxSigmaX_),
sigmax_[faceI],
R_[faceI],
sigmax_[patchFaceI],
R[patchFaceI],
rndGen_
);
@ -526,16 +397,21 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
WarningInFunction
<< "Patch: " << patch().patch().name()
<< " on field " << internalField().name()
<< ": No eddies seeded - please check your set-up" << endl;
<< ": No eddies seeded - please check your set-up"
<< endl;
}
}
void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
(
const vector& UBulk,
const scalar deltaT
)
{
const scalar t = db().time().timeOutputValue();
const symmTensorField R(R_->value(t)/sqr(Uref_));
// Note: all operations applied to local processor only
label nRecycled = 0;
@ -543,7 +419,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
forAll(eddies_, eddyI)
{
eddy& e = eddies_[eddyI];
e.move(deltaT*(UMean_ & patchNormal_));
e.move(deltaT*(UBulk & patchNormal_));
const scalar position0 = e.x();
@ -555,17 +431,17 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
while (search && iter++ < seedIterMax_)
{
// Spawn new eddy with new random properties (intensity etc)
pointIndexHit pos(setNewPosition(false));
label faceI = pos.index();
// Spawn new eddy with new random properties (intensity etc)
pointIndexHit pos(setNewPosition(false));
const label patchFaceI = pos.index();
e = eddy
e = eddy
(
faceI,
patchFaceI,
pos.hitPoint(),
position0 - floor(position0/maxSigmaX_)*maxSigmaX_,
sigmax_[faceI],
R_[faceI],
sigmax_[patchFaceI],
R[patchFaceI],
rndGen_
);
@ -583,27 +459,28 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
if (debug && nRecycled > 0)
{
Info<< "Patch: " << patch().patch().name() << " recycled "
<< nRecycled << " eddies" << endl;
Info<< "Patch: " << patch().patch().name()
<< " recycled " << nRecycled << " eddies"
<< endl;
}
}
Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uDashEddy
Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uPrimeEddy
(
const List<eddy>& eddies,
const point& patchFaceCf
) const
{
vector uDash(Zero);
vector uPrime(Zero);
forAll(eddies, k)
{
const eddy& e = eddies[k];
uDash += e.uDash(patchFaceCf, patchNormal_);
uPrime += e.uPrime(patchFaceCf, patchNormal_);
}
return uDash;
return uPrime;
}
@ -629,8 +506,8 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
const eddy& e = eddies_[i];
// Eddy bounds
point x = e.position(patchNormal_);
boundBox ebb = e.bounds();
const point x(e.position(patchNormal_));
boundBox ebb(e.bounds());
ebb.min() += x;
ebb.max() += x;
@ -678,10 +555,10 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
if (procI != Pstream::myProcNo())
{
// What I need to receive is what other processor is sending to me
label nRecv = sendSizes[procI][Pstream::myProcNo()];
const label nRecv = sendSizes[procI][Pstream::myProcNo()];
constructMap[procI].setSize(nRecv);
for (label i = 0; i < nRecv; i++)
for (label i = 0; i < nRecv; ++i)
{
constructMap[procI][i] = segmentI++;
}
@ -738,37 +615,33 @@ turbulentDFSEMInletFvPatchVectorField
)
:
fixedValueFvPatchField<vector>(p, iF),
delta_(Zero),
d_(Zero),
kappa_(Zero),
perturb_(1e-5),
mapMethod_("nearestCell"),
mapperPtr_(nullptr),
interpolateR_(false),
interpolateL_(false),
interpolateU_(false),
R_(),
L_(),
U_(),
UMean_(Zero),
U_(nullptr),
R_(nullptr),
L_(nullptr),
delta_(1.0),
d_(1.0),
kappa_(0.41),
Uref_(1.0),
Lref_(1.0),
scale_(1.0),
m_(0.5),
nCellPerEddy_(5),
patchArea_(-1),
triFace_(),
triToFace_(),
triCumulativeMagSf_(),
sumTriMagSf_(Pstream::nProcs() + 1, Zero),
patchNormal_(Zero),
patchBounds_(boundBox::invertedBox),
eddies_(Zero),
nCellPerEddy_(5),
patchNormal_(Zero),
v0_(Zero),
rndGen_(Pstream::myProcNo()),
sigmax_(size(), Zero),
maxSigmaX_(Zero),
nEddy_(Zero),
nEddy_(0),
curTimeIndex_(-1),
patchBounds_(boundBox::invertedBox),
singleProc_(false),
writeEddies_(false)
{}
@ -784,37 +657,33 @@ turbulentDFSEMInletFvPatchVectorField
)
:
fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
U_(ptf.U_.clone(patch().patch())),
R_(ptf.R_.clone(patch().patch())),
L_(ptf.L_.clone(patch().patch())),
delta_(ptf.delta_),
d_(ptf.d_),
kappa_(ptf.kappa_),
perturb_(ptf.perturb_),
mapMethod_(ptf.mapMethod_),
mapperPtr_(nullptr),
interpolateR_(ptf.interpolateR_),
interpolateL_(ptf.interpolateL_),
interpolateU_(ptf.interpolateU_),
R_(ptf.R_, mapper),
L_(ptf.L_, mapper),
U_(ptf.U_, mapper),
UMean_(ptf.UMean_),
Uref_(ptf.Uref_),
Lref_(ptf.Lref_),
scale_(ptf.scale_),
m_(ptf.m_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchArea_(ptf.patchArea_),
triFace_(ptf.triFace_),
triToFace_(ptf.triToFace_),
triCumulativeMagSf_(ptf.triCumulativeMagSf_),
sumTriMagSf_(ptf.sumTriMagSf_),
patchNormal_(ptf.patchNormal_),
patchBounds_(ptf.patchBounds_),
eddies_(ptf.eddies_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchNormal_(ptf.patchNormal_),
v0_(ptf.v0_),
rndGen_(ptf.rndGen_),
sigmax_(ptf.sigmax_, mapper),
maxSigmaX_(ptf.maxSigmaX_),
nEddy_(Zero),
nEddy_(ptf.nEddy_),
curTimeIndex_(-1),
patchBounds_(ptf.patchBounds_),
singleProc_(ptf.singleProc_),
writeEddies_(ptf.writeEddies_)
{}
@ -829,46 +698,42 @@ turbulentDFSEMInletFvPatchVectorField
)
:
fixedValueFvPatchField<vector>(p, iF, dict),
delta_(dict.get<scalar>("delta")),
d_(dict.getOrDefault<scalar>("d", 1)),
kappa_(dict.getOrDefault<scalar>("kappa", 0.41)),
perturb_(dict.getOrDefault<scalar>("perturb", 1e-5)),
mapMethod_(dict.getOrDefault<word>("mapMethod", "nearestCell")),
mapperPtr_(nullptr),
interpolateR_(dict.getOrDefault("interpolateR", false)),
interpolateL_(dict.getOrDefault("interpolateL", false)),
interpolateU_(dict.getOrDefault("interpolateU", false)),
R_(interpolateOrRead<symmTensor>("R", dict, interpolateR_)),
L_(interpolateOrRead<scalar>("L", dict, interpolateL_)),
U_(interpolateOrRead<vector>("U", dict, interpolateU_)),
UMean_(Zero),
U_(PatchFunction1<vector>::New(patch().patch(), "U", dict)),
R_(PatchFunction1<symmTensor>::New(patch().patch(), "R", dict)),
L_(PatchFunction1<scalar>::New(patch().patch(), "L", dict)),
delta_(dict.getCheck<scalar>("delta", scalarMinMax::ge(0))),
d_(dict.getCheckOrDefault<scalar>("d", 1, scalarMinMax::ge(SMALL))),
kappa_(dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(0))),
Uref_(dict.getCheckOrDefault<scalar>("Uref", 1, scalarMinMax::ge(SMALL))),
Lref_(dict.getCheckOrDefault<scalar>("Lref", 1, scalarMinMax::ge(SMALL))),
scale_(dict.getCheckOrDefault<scalar>("scale", 1, scalarMinMax::ge(0))),
m_(dict.getCheckOrDefault<scalar>("m", 0.5, scalarMinMax::ge(0))),
nCellPerEddy_(dict.getOrDefault<label>("nCellPerEddy", 5)),
patchArea_(-1),
triFace_(),
triToFace_(),
triCumulativeMagSf_(),
sumTriMagSf_(Pstream::nProcs() + 1, Zero),
patchNormal_(Zero),
patchBounds_(boundBox::invertedBox),
eddies_(),
nCellPerEddy_(dict.getOrDefault<label>("nCellPerEddy", 5)),
patchNormal_(Zero),
v0_(Zero),
rndGen_(),
sigmax_(size(), Zero),
maxSigmaX_(Zero),
nEddy_(Zero),
nEddy_(0),
curTimeIndex_(-1),
patchBounds_(boundBox::invertedBox),
singleProc_(false),
writeEddies_(dict.getOrDefault("writeEddies", false))
{
eddy::debug = debug;
checkStresses(R_);
const scalar t = db().time().timeOutputValue();
const symmTensorField R(R_->value(t)/sqr(Uref_));
// Set UMean as patch area average value
UMean_ = gSum(U_*patch().magSf())/(gSum(patch().magSf()) + ROOTVSMALL);
checkStresses(R);
}
@ -879,37 +744,33 @@ turbulentDFSEMInletFvPatchVectorField
)
:
fixedValueFvPatchField<vector>(ptf),
U_(ptf.U_.clone(patch().patch())),
R_(ptf.R_.clone(patch().patch())),
L_(ptf.L_.clone(patch().patch())),
delta_(ptf.delta_),
d_(ptf.d_),
kappa_(ptf.kappa_),
perturb_(ptf.perturb_),
mapMethod_(ptf.mapMethod_),
mapperPtr_(nullptr),
interpolateR_(ptf.interpolateR_),
interpolateL_(ptf.interpolateL_),
interpolateU_(ptf.interpolateU_),
R_(ptf.R_),
L_(ptf.L_),
U_(ptf.U_),
UMean_(ptf.UMean_),
Uref_(ptf.Uref_),
Lref_(ptf.Lref_),
scale_(ptf.scale_),
m_(ptf.m_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchArea_(ptf.patchArea_),
triFace_(ptf.triFace_),
triToFace_(ptf.triToFace_),
triCumulativeMagSf_(ptf.triCumulativeMagSf_),
sumTriMagSf_(ptf.sumTriMagSf_),
patchNormal_(ptf.patchNormal_),
patchBounds_(ptf.patchBounds_),
eddies_(ptf.eddies_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchNormal_(ptf.patchNormal_),
v0_(ptf.v0_),
rndGen_(ptf.rndGen_),
sigmax_(ptf.sigmax_),
maxSigmaX_(ptf.maxSigmaX_),
nEddy_(Zero),
nEddy_(ptf.nEddy_),
curTimeIndex_(-1),
patchBounds_(ptf.patchBounds_),
singleProc_(ptf.singleProc_),
writeEddies_(ptf.writeEddies_)
{}
@ -923,37 +784,33 @@ turbulentDFSEMInletFvPatchVectorField
)
:
fixedValueFvPatchField<vector>(ptf, iF),
U_(ptf.U_.clone(patch().patch())),
R_(ptf.R_.clone(patch().patch())),
L_(ptf.L_.clone(patch().patch())),
delta_(ptf.delta_),
d_(ptf.d_),
kappa_(ptf.kappa_),
perturb_(ptf.perturb_),
mapMethod_(ptf.mapMethod_),
mapperPtr_(nullptr),
interpolateR_(ptf.interpolateR_),
interpolateL_(ptf.interpolateL_),
interpolateU_(ptf.interpolateU_),
R_(ptf.R_),
L_(ptf.L_),
U_(ptf.U_),
UMean_(ptf.UMean_),
Uref_(ptf.Uref_),
Lref_(ptf.Lref_),
scale_(ptf.scale_),
m_(ptf.m_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchArea_(ptf.patchArea_),
triFace_(ptf.triFace_),
triToFace_(ptf.triToFace_),
triCumulativeMagSf_(ptf.triCumulativeMagSf_),
sumTriMagSf_(ptf.sumTriMagSf_),
patchNormal_(ptf.patchNormal_),
patchBounds_(ptf.patchBounds_),
eddies_(ptf.eddies_),
nCellPerEddy_(ptf.nCellPerEddy_),
patchNormal_(ptf.patchNormal_),
v0_(ptf.v0_),
rndGen_(ptf.rndGen_),
sigmax_(ptf.sigmax_),
maxSigmaX_(ptf.maxSigmaX_),
nEddy_(Zero),
nEddy_(ptf.nEddy_),
curTimeIndex_(-1),
patchBounds_(ptf.patchBounds_),
singleProc_(ptf.singleProc_),
writeEddies_(ptf.writeEddies_)
{}
@ -981,11 +838,11 @@ bool Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
<< exit(FatalError);
}
scalar a_xx = sqrt(R.xx());
const scalar a_xx = sqrt(R.xx());
scalar a_xy = R.xy()/a_xx;
const scalar a_xy = R.xy()/a_xx;
scalar a_yy_2 = R.yy() - sqr(a_xy);
const scalar a_yy_2 = R.yy() - sqr(a_xy);
if (a_yy_2 < 0)
{
@ -996,13 +853,13 @@ bool Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
<< exit(FatalError);
}
scalar a_yy = Foam::sqrt(a_yy_2);
const scalar a_yy = Foam::sqrt(a_yy_2);
scalar a_xz = R.xz()/a_xx;
const scalar a_xz = R.xz()/a_xx;
scalar a_yz = (R.yz() - a_xy*a_xz)*a_yy;
const scalar a_yz = (R.yz() - a_xy*a_xz)/a_yy;
scalar a_zz_2 = R.zz() - sqr(a_xz) - sqr(a_yz);
const scalar a_zz_2 = R.zz() - sqr(a_xz) - sqr(a_yz);
if (a_zz_2 < 0)
{
@ -1013,7 +870,7 @@ bool Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
<< exit(FatalError);
}
scalar a_zz = Foam::sqrt(a_zz_2);
const scalar a_zz = Foam::sqrt(a_zz_2);
if (debug)
{
@ -1036,11 +893,18 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::autoMap
{
fixedValueFvPatchField<vector>::autoMap(m);
// Clear interpolator
mapperPtr_.clear();
R_.autoMap(m);
L_.autoMap(m);
U_.autoMap(m);
if (U_)
{
U_->autoMap(m);
}
if (R_)
{
R_->autoMap(m);
}
if (L_)
{
L_->autoMap(m);
}
sigmax_.autoMap(m);
}
@ -1054,15 +918,21 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::rmap
{
fixedValueFvPatchField<vector>::rmap(ptf, addr);
const turbulentDFSEMInletFvPatchVectorField& dfsemptf =
const auto& dfsemptf =
refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
R_.rmap(dfsemptf.R_, addr);
L_.rmap(dfsemptf.L_, addr);
U_.rmap(dfsemptf.U_, addr);
// Clear interpolator
mapperPtr_.clear();
if (U_)
{
U_->rmap(dfsemptf.U_(), addr);
}
if (R_)
{
R_->rmap(dfsemptf.R_(), addr);
}
if (L_)
{
L_->rmap(dfsemptf.L_(), addr);
}
sigmax_.rmap(dfsemptf.sigmax_, addr);
}
@ -1087,38 +957,38 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
if (curTimeIndex_ != db().time().timeIndex())
{
if (debug)
{
label n = eddies_.size();
Info<< "Number of eddies: " << returnReduce(n, sumOp<label>())
<< endl;
}
tmp<vectorField> UMean =
U_->value(db().time().timeOutputValue())/Uref_;
// (PCR:p. 522)
const vector UBulk
(
gSum(UMean()*patch().magSf())
/(gSum(patch().magSf()) + ROOTVSMALL)
);
// Move eddies using bulk velocity
const scalar deltaT = db().time().deltaTValue();
convectEddies(UBulk, deltaT);
// Move eddies using mean velocity
convectEddies(deltaT);
// Set velocity
// Set mean velocity
vectorField& U = *this;
//U = UMean_;
U = U_;
const pointField& Cf = patch().Cf();
U = UMean;
// Apply second part of normalisation coefficient
// Note: factor of 2 required to match reference stresses?
const scalar FACTOR = 2;
const scalar c = FACTOR*Foam::sqrt(10*v0_)/Foam::sqrt(scalar(nEddy_));
const scalar c =
scale_*Foam::pow(10*v0_, m_)/Foam::sqrt(scalar(nEddy_));
// In parallel, need to collect all eddies that will interact with
// local faces
const pointField& Cf = patch().Cf();
if (singleProc_ || !Pstream::parRun())
{
forAll(U, faceI)
{
U[faceI] += c*uDashEddy(eddies_, Cf[faceI]);
U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
}
}
else
@ -1126,7 +996,7 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
// Process local eddy contributions
forAll(U, faceI)
{
U[faceI] += c*uDashEddy(eddies_, Cf[faceI]);
U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
}
// Add contributions from overlapping eddies
@ -1139,30 +1009,21 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
if (eddies.size())
{
//Pout<< "Applying " << eddies.size()
// << " eddies from processor " << procI << endl;
forAll(U, faceI)
{
U[faceI] += c*uDashEddy(eddies, Cf[faceI]);
U[faceI] += c*uPrimeEddy(eddies, Cf[faceI]);
}
}
}
}
// Re-scale to ensure correct flow rate
scalar fCorr =
gSum((UMean_ & patchNormal_)*patch().magSf())
const scalar fCorr =
gSum((UBulk & patchNormal_)*patch().magSf())
/gSum(U & -patch().Sf());
U *= fCorr;
if (debug)
{
Info<< "Patch:" << patch().patch().name()
<< " min/max(U):" << gMin(U) << ", " << gMax(U) << endl;
}
curTimeIndex_ = db().time().timeIndex();
if (writeEddies_)
@ -1170,9 +1031,22 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
writeEddyOBJ();
}
if (debug && db().time().writeTime())
if (debug)
{
writeLumleyCoeffs();
Info<< "Magnitude of bulk velocity: " << UBulk << endl;
label n = eddies_.size();
Info<< "Number of eddies: " << returnReduce(n, sumOp<label>())
<< endl;
Info<< "Patch:" << patch().patch().name()
<< " min/max(U):" << gMin(U) << ", " << gMax(U)
<< endl;
if (db().time().writeTime())
{
writeLumleyCoeffs();
}
}
}
@ -1183,38 +1057,28 @@ void Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs()
void Foam::turbulentDFSEMInletFvPatchVectorField::write(Ostream& os) const
{
fvPatchField<vector>::write(os);
writeEntry("value", os);
os.writeEntry("delta", delta_);
os.writeEntryIfDifferent<scalar>("d", 1.0, d_);
os.writeEntryIfDifferent<scalar>("kappa", 0.41, kappa_);
os.writeEntryIfDifferent<scalar>("perturb", 1e-5, perturb_);
os.writeEntryIfDifferent<scalar>("Uref", 1.0, Uref_);
os.writeEntryIfDifferent<scalar>("Lref", 1.0, Lref_);
os.writeEntryIfDifferent<scalar>("scale", 1.0, scale_);
os.writeEntryIfDifferent<scalar>("m", 0.5, m_);
os.writeEntryIfDifferent<label>("nCellPerEddy", 5, nCellPerEddy_);
os.writeEntryIfDifferent("writeEddies", false, writeEddies_);
if (!interpolateR_)
if (U_)
{
R_.writeEntry("R", os);
U_->writeData(os);
}
if (!interpolateL_)
if (R_)
{
L_.writeEntry("L", os);
R_->writeData(os);
}
if (!interpolateU_)
if (L_)
{
U_.writeEntry("U", os);
}
if (!mapMethod_.empty())
{
os.writeEntryIfDifferent<word>
(
"mapMethod",
"nearestCell",
mapMethod_
);
L_->writeData(os);
}
writeEntry("value", os);
}

296
src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorField.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -31,56 +31,107 @@ Group
grpInletBoundaryConditions
Description
Velocity boundary condition including synthesised eddies for use with LES
and DES turbulent flows.
The \c turbulentDFSEMInlet is a synthesised-eddy based velocity inlet
boundary condition to generate synthetic turbulence-alike time-series
from a given set of turbulence statistics for LES and hybrid RANS-LES
computations.
Reference:
\verbatim
Poletto, R., Craft, T., & Revell, A. (2013).
A new divergence free synthetic eddy method for
the reproduction of inlet flow conditions for LES.
Flow, turbulence and combustion, 91(3), 519-539.
DOI:10.1007/s10494-013-9488-2
\endverbatim
Standard model (tag:PCR):
Poletto, R., Craft, T., & Revell, A. (2013).
A new divergence free synthetic eddy method for
the reproduction of inlet flow conditions for LES.
Flow, turbulence and combustion, 91(3), 519-539.
DOI:10.1007/s10494-013-9488-2
Reynolds stress, velocity and turbulence length scale values can either
be specified directly, or mapped. If mapping, the values should be
entered in the same form as the \c timeVaryingMappedFixedValue condition,
except that no interpolation in time is supported. These should be
located in the directory:
\verbatim
\$FOAM_CASE/constant/boundaryData/\<patchName\>/points
\$FOAM_CASE/constant/boundaryData/\<patchName\>/0/\{R|U|L\}
Expression for the average length scale (tag:SST):
Shur, M., Strelets, M., Travin, A.,
Probst, A., Probst, S., Schwamborn, D., ... & Revell, A. (2018).
Improved embedded approaches.
In: Mockett C., Haase W., Schwamborn D. (eds)
Go4Hybrid: Grey area mitigation for hybrid RANS-LES methods.
Notes on Numerical Fluid Mechanics and Multidisciplinary Design.
p. 51-87. Springer, Cham.
DOI:10.1007/978-3-319-52995-0_3
\endverbatim
Usage
Example of the boundary condition specification:
\verbatim
<patchName>
{
// Mandatory entries
type turbulentDFSEMInlet;
delta <scalar>;
R <PatchFunction1>;
U <PatchFunction1>;
L <PatchFunction1>;
// e.g.
// R uniform (<Rxx> <Rxy> <Rxz> <Ryy> <Ryz> <Rzz>);
// U uniform (<Ux> <Uy> <Uz>);
// L uniform <L>;
// Optional entries
d <scalar>;
nCellPerEddy <label>;
kappa <scalar>;
Uref <scalar>;
Lref <scalar>;
scale <scalar>;
m <scalar>;
writeEddies <bool>;
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Required | Default value
value | Restart value | yes |
delta | Local limiting length scale | yes |
R | Reynolds stress field | no |
U | Velocity field | no |
L | Turbulence length scale field | no |
d | Eddy density (fill fraction) | no | 1
kappa | Von Karman constant | no | 0.41
mapMethod | Method to map reference values | no | nearestCell
perturb | Point perturbation for interpolation | no | 1e-5
interpolateR | Flag to interpolate the R field | no | false
interpolateL | Flag to interpolate the L field | no | false
interpolateU | Flag to interpolate the U field | no | false
writeEddies | Flag to write eddies as OBJ file | no | no
Property | Description | Type | Reqd | Deflt
type | Type name: turbulentDFSEMInlet | word | yes | -
delta | Characteristic length scale | scalar | yes | -
R | Reynolds-stress tensor field <!--
--> | PatchFunction1\<symmTensor\> | yes | -
U | Mean velocity field <!--
--> | PatchFunction1<vector> | yes | -
L | Integral-length scale field <!--
--> | PatchFunction1<scalar> | yes | -
d | Ratio of sum of eddy volumes to eddy box volume <!--
--> i.e. eddy density (fill fraction) | scalar | no | 1.0
nCellPerEddy | Minimum eddy length in units of number of cells <!--
--> | label | no | 5
kappa | von Karman constant | scalar | no | 0.41
Uref | Normalisation factor for Reynolds-stress <!--
--> tensor and mean velocity | scalar | no | 1.0
Lref | Normalisation factor for integral-length scale <!--
--> | scalar | no | 1.0
scale | Heuristic scaling factor being applied <!--
--> on the normalisation factor C1 | scalar | no | 1.0
m | The power of V defined in C1 | scalar | no | 0.5
writeEddies | Flag to write eddies as OBJ file | bool | no | false
\endtable
Note
- The \c delta value typically represents the characteristic scale of flow
or flow domain, e.g. a channel half-height
- For \c R, \c U and \c L specification: if the entry is not user input,
it is assumed that the data will be mapped
The inherited entries are elaborated in:
- \link fixedValueFvPatchFields.H \endlink
- \link PatchFunction1.H \endlink
- \link MappedFile.H \endlink
SeeAlso
timeVaryingMappedFixedValueFvPatchField
turbulentDigitalFilterInlet
Note
- The \c delta value typically represents the characteristic scale of flow
or flow domain, e.g. a channel half height for plane channel flows.
- \c nCellPerEddy and \c scale entries are heuristic entries
which do not exist in the standard method, yet are necessary
to reproduce the results published in the original journal paper.
- In the original journal paper, \c C1 in Eq. 11 is not dimensionless.
It is not clear whether this dimensionality issue was intentional.
To alleviate this matter, users can alter the input entry \c m, which is
the power of the eddy-box volume defined in the \c C1, to obtain a
dimensionless \c C1 coefficient. The default value of \c m is 0.5,
which is the value stated in the original journal paper,
and \c m=1/3 leads to a dimensionless \c C1.
SourceFiles
turbulentDFSEMInletFvPatchVectorField.C
@ -94,15 +145,13 @@ SourceFiles
#include "Random.H"
#include "eddy.H"
#include "pointIndexHit.H"
#include "instantList.H"
#include "PatchFunction1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class pointToPointPlanarInterpolation;
/*---------------------------------------------------------------------------*\
Class turbulentDFSEMInletFvPatchVectorField Declaration
\*---------------------------------------------------------------------------*/
@ -116,50 +165,38 @@ class turbulentDFSEMInletFvPatchVectorField
//- Maximum number of attempts when seeding eddies
static label seedIterMax_;
//- Characteristic length scale, e.g. half channel height
//- Mean velocity field
autoPtr<PatchFunction1<vector>> U_;
//- Reynolds stress tensor field
autoPtr<PatchFunction1<symmTensor>> R_;
//- Integral-length scale field
autoPtr<PatchFunction1<scalar>> L_;
//- Characteristic length scale
const scalar delta_;
//- Ratio of sum of eddy volumes to eddy box volume; default = 1
//- Ratio of sum of eddy volumes to eddy box volume, i.e. eddy density
const scalar d_;
//- Von Karman constant
//- von Karman constant
const scalar kappa_;
//- Global numbering for faces
mutable autoPtr<globalIndex> globalFacesPtr_;
//- Normalisation factor for Reynolds-stress tensor and mean velocity
const scalar Uref_;
//- Normalisation factor for integral-length scale
const scalar Lref_;
// Table reading for patch inlet flow properties
//- Heuristic scaling factor being applied on the normalisation factor
const scalar scale_;
//- Fraction of perturbation (fraction of bounding box) to add
scalar perturb_;
//- The power of V defined in C1
const scalar m_;
//- Interpolation scheme to use (nearestCell | planarInterpolation)
word mapMethod_;
//- 2D interpolation (for 'planarInterpolation' mapMethod)
mutable autoPtr<pointToPointPlanarInterpolation> mapperPtr_;
//- Flag to identify to interpolate the R field
bool interpolateR_;
//- Flag to identify to interpolate the L field
bool interpolateL_;
//- Flag to identify to interpolate the U field
bool interpolateU_;
//- Reynolds stress tensor
symmTensorField R_;
//- Length scale
scalarField L_;
//- Inlet velocity
vectorField U_;
//- Mean inlet velocity
vector UMean_;
//- Minimum eddy length in units of number of cells
const label nCellPerEddy_;
// Patch information
@ -179,46 +216,51 @@ class turbulentDFSEMInletFvPatchVectorField
//- Cumulative area fractions per processor
scalarList sumTriMagSf_;
//- Patch normal into the domain
vector patchNormal_;
//- List of eddies
List<eddy> eddies_;
//- Patch bounds (local processor)
boundBox patchBounds_;
//- Minimum number of cells required to resolve an eddy
label nCellPerEddy_;
//- Patch normal into the domain
vector patchNormal_;
// Eddy information
//- Eddy box volume
scalar v0_;
//- List of eddies
List<eddy> eddies_;
//- Random number generator
Random rndGen_;
//- Eddy box volume
scalar v0_;
//- Length scale per patch face
scalarField sigmax_;
//- Random number generator
Random rndGen_;
//- Maximum length scale (across all processors)
scalar maxSigmaX_;
//- Integral-length scale per patch face
scalarField sigmax_;
//- Global number of eddies
label nEddy_;
//- Maximum integral-length scale (across all processors)
scalar maxSigmaX_;
//- Current time index (used for updating)
label curTimeIndex_;
//- Global number of eddies
label nEddy_;
//- Patch bounds (local processor)
boundBox patchBounds_;
//- Current time index (used for updating)
label curTimeIndex_;
//- Single processor contains all eddies (flag)
bool singleProc_;
//- Single processor contains all eddies (flag)
bool singleProc_;
//- Flag to write the eddies to file
bool writeEddies_;
//- Flag to write the eddies to file
bool writeEddies_;
// Private Member Functions
//- Write Lumley coefficients to file
void writeLumleyCoeffs() const;
//- Write eddy info in OBJ format
void writeEddyOBJ() const;
//- Initialise info for patch point search
void initialisePatch();
@ -231,37 +273,11 @@ class turbulentDFSEMInletFvPatchVectorField
//- Initialise eddies
void initialiseEddies();
//- Convect the eddies
void convectEddies(const scalar deltaT);
//- Convect the eddies with the bulk velocity
void convectEddies(const vector& UBulk, const scalar deltaT);
//- Calculate the velocity fluctuation at a point
vector uDashEddy(const List<eddy>& eddies, const point& globalX) const;
//- Helper function to interpolate values from the boundary data or
//- read from dictionary
template<class Type>
tmp<Field<Type>> interpolateOrRead
(
const word& fieldName,
const dictionary& dict,
bool& interpolateField
) const;
//- Helper function to interpolate values from the boundary data
template<class Type>
tmp<Field<Type>> interpolateBoundaryData
(
const word& fieldName
) const;
//- Write Lumley coefficients to file
void writeLumleyCoeffs() const;
//- Write eddy info in OBJ format
void writeEddyOBJ() const;
//- Return a reference to the patch mapper object
const pointToPointPlanarInterpolation& patchMapper() const;
//- Return velocity fluctuation vector at a given point
vector uPrimeEddy(const List<eddy>& eddies, const point& globalX) const;
//- Return eddies from remote processors that interact with local
//- processor
@ -345,11 +361,11 @@ public:
// Member Functions
//- Helper function to check that Reynold stresses are valid
//- Return true if input Reynold stresses are valid
static bool checkStresses(const symmTensorField& Rf);
// Mapping functions
// Mapping
//- Map (and resize as needed) from self given a mapping object
virtual void autoMap(const fvPatchFieldMapper& m);
@ -362,14 +378,16 @@ public:
);
// Evaluation functions
// Evaluation
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
// IO
//- Write
virtual void write(Ostream&) const;
};
@ -379,12 +397,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "turbulentDFSEMInletFvPatchVectorFieldTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

109
src/finiteVolume/fields/fvPatchFields/derived/turbulentDFSEMInlet/turbulentDFSEMInletFvPatchVectorFieldTemplates.C
View File

@ -1,109 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2016-2020 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 "pointToPointPlanarInterpolation.H"
#include "Time.H"
#include "rawIOField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::turbulentDFSEMInletFvPatchVectorField::interpolateOrRead
(
const word& fieldName,
const dictionary& dict,
bool& interpolateField
) const
{
if (dict.found(fieldName))
{
tmp<Field<Type>> tFld
(
new Field<Type>
(
fieldName,
dict,
this->patch().size()
)
);
interpolateField = false;
return tFld;
}
else
{
interpolateField = true;
return interpolateBoundaryData<Type>(fieldName);
}
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::turbulentDFSEMInletFvPatchVectorField::interpolateBoundaryData
(
const word& fieldName
) const
{
const word& patchName = this->patch().name();
const fileName valsFile
(
fileName
(
this->db().time().globalPath()
/this->db().time().constant()
/"boundaryData"
/patchName
/"0"
/fieldName
)
);
IOobject io
(
valsFile, // absolute path
this->db().time(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false, // no need to register
true // is global object (currently not used)
);
const rawIOField<Type> vals(io, false);
Info<< "Turbulent DFSEM patch " << patchName
<< ": interpolating field " << fieldName
<< " from " << valsFile << endl;
return patchMapper().interpolate(vals);
}
// ************************************************************************* //