diff --git a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C b/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C
index 68f0c2a5bf926e1812886b95de3a6ea9f5bfb42e..878398220ef20ce73d72a8a7c823a6d71b2c525b 100644
--- a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C
+++ b/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C
@@ -66,6 +66,12 @@ void Foam::actuationDiskSource::checkData() const
<< "disk direction vector is approximately zero"
<< exit(FatalIOError);
}
+ if (returnReduce(upstreamCellId_, maxOp<label>()) == -1)
+ {
+ FatalErrorIn("Foam::actuationDiskSource::checkData()")
+ << "upstream location " << upstreamPoint_ << " not found in mesh"
+ << exit(FatalIOError);
+ }
}
@@ -83,7 +89,9 @@ Foam::actuationDiskSource::actuationDiskSource
diskDir_(coeffs_.lookup("diskDir")),
Cp_(readScalar(coeffs_.lookup("Cp"))),
Ct_(readScalar(coeffs_.lookup("Ct"))),
- diskArea_(readScalar(coeffs_.lookup("diskArea")))
+ diskArea_(readScalar(coeffs_.lookup("diskArea"))),
+ upstreamPoint_(coeffs_.lookup("upstreamPoint")),
+ upstreamCellId_(-1)
{
coeffs_.lookup("fieldNames") >> fieldNames_;
applied_.setSize(fieldNames_.size(), false);
@@ -91,6 +99,8 @@ Foam::actuationDiskSource::actuationDiskSource
Info<< " - creating actuation disk zone: "
<< this->name() << endl;
+ upstreamCellId_ = mesh.findCell(upstreamPoint_);
+
checkData();
}
diff --git a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H b/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H
index 077442589258e6d73bdf45a31d511af6225abb57..80cd7bd54a88c2da40b4580de14b90308cdf3b48 100644
--- a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H
+++ b/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H
@@ -48,6 +48,7 @@ Description
Cp 0.1; // power coefficient
Ct 0.5; // thrust coefficient
diskArea 5.0; // disk area
+ upstreamPoint (0 0 0); // upstream point
}
@@ -92,6 +93,12 @@ protected:
//- Disk area
scalar diskArea_;
+ //- Upstream point sample
+ point upstreamPoint_;
+
+ //- Upstream cell ID
+ label upstreamCellId_;
+
private:
diff --git a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C b/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
index dbaf95c4b94d77bf14ae3e547462dbc59b40cd2e..959e86f5b3fedf05088ef1ee7bcb739e13605add 100644
--- a/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
+++ b/src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C
@@ -42,20 +42,27 @@ void Foam::actuationDiskSource::addActuationDiskAxialInertialResistance
) const
{
scalar a = 1.0 - Cp_/Ct_;
- scalarField T(cells.size());
vector uniDiskDir = diskDir_/mag(diskDir_);
tensor E(tensor::zero);
E.xx() = uniDiskDir.x();
E.yy() = uniDiskDir.y();
E.zz() = uniDiskDir.z();
- forAll(cells, i)
+ vector upU = vector(VGREAT, VGREAT, VGREAT);
+ scalar upRho = VGREAT;
+ if (upstreamCellId_ != -1)
{
- T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a);
+ upU = U[upstreamCellId_];
+ upRho = rho[upstreamCellId_];
}
+ reduce(upU, minOp<vector>());
+ reduce(upRho, minOp<scalar>());
+
+ scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1 - a);
+
forAll(cells, i)
{
- Usource[cells[i]] -= ((Vcells[cells[i]]/V())*T[i]*E) & U[cells[i]];
+ Usource[cells[i]] += ((Vcells[cells[i]]/V())*T*E) & upU;
}
}
diff --git a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H b/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H
index 4c464bcfd36754571d947eca2af534ca456a8969..134e93b5eae64116d38e0e8048a29980598a8f27 100644
--- a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H
+++ b/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H
@@ -53,6 +53,7 @@ Description
Ct 0.5; // thrust coefficient
diskArea 5.0; // disk area
coeffs (0.1 0.5 0.01); // radial distribution coefficients
+ upstreamPoint (0 0 0); // upstream point
}
diff --git a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C b/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
index 127263c94f2ff5dbc937020e6d8e45c862678b76..adea6154550cf1642d52834d7af47b26efe18d64 100644
--- a/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
+++ b/src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C
@@ -44,11 +44,9 @@ addRadialActuationDiskAxialInertialResistance
) const
{
scalar a = 1.0 - Cp_/Ct_;
- scalarField T(cells.size());
scalarField Tr(cells.size());
const vector uniDiskDir = diskDir_/mag(diskDir_);
-
tensor E(tensor::zero);
E.xx() = uniDiskDir.x();
E.yy() = uniDiskDir.y();
@@ -65,21 +63,27 @@ addRadialActuationDiskAxialInertialResistance
+ radialCoeffs_[1]*sqr(maxR)/2.0
+ radialCoeffs_[2]*pow4(maxR)/3.0;
- forAll(cells, i)
+ vector upU = vector(VGREAT, VGREAT, VGREAT);
+ scalar upRho = VGREAT;
+ if (upstreamCellId_ != -1)
{
- T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a);
+ upU = U[upstreamCellId_];
+ upRho = rho[upstreamCellId_];
+ }
+ reduce(upU, minOp<vector>());
+ reduce(upRho, minOp<scalar>());
+ scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1.0 - a);
+ forAll(cells, i)
+ {
scalar r2 = magSqr(mesh().cellCentres()[cells[i]] - avgCentre);
Tr[i] =
- T[i]
+ T
*(radialCoeffs_[0] + radialCoeffs_[1]*r2 + radialCoeffs_[2]*sqr(r2))
/intCoeffs;
- }
- forAll(cells, i)
- {
- Usource[cells[i]] -= ((Vcells[cells[i]]/V_)*Tr[i]*E) & U[cells[i]];
+ Usource[cells[i]] += ((Vcells[cells[i]]/V_)*Tr[i]*E) & upU;
}
if (debug)
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
index a3d313e30285ffb9eaf10ac5c4f658518f6abec0..4b4273ff32ffcbb39accd680c50c2ac391ec8496 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C
@@ -51,6 +51,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
z_(pTraits<vector>::zero),
z0_(0),
kappa_(0.41),
+ Uref_(0),
+ Href_(0),
zGround_(0)
{}
@@ -69,6 +71,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
z_(ptf.z_),
z0_(ptf.z0_),
kappa_(ptf.kappa_),
+ Uref_(ptf.Uref_),
+ Href_(ptf.Href_),
zGround_(ptf.zGround_)
{}
@@ -82,10 +86,12 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
)
:
fixedValueFvPatchScalarField(p, iF),
- Ustar_(readScalar(dict.lookup("Ustar"))),
+ Ustar_(p.size()),
z_(dict.lookup("z")),
z0_(readScalar(dict.lookup("z0"))),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+ Uref_(readScalar(dict.lookup("Uref"))),
+ Href_(readScalar(dict.lookup("Href"))),
zGround_("zGround", dict, p.size())
{
if (mag(z_) < SMALL)
@@ -103,6 +109,11 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
<< abort(FatalError);
}
+ forAll (Ustar_, i)
+ {
+ Ustar_[i] = kappa_*Uref_/(log((Href_ + z0_[i])/max(z0_[i] , 0.001)));
+ }
+
z_ /= mag(z_);
evaluate();
@@ -121,6 +132,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
z_(blpsf.z_),
z0_(blpsf.z0_),
kappa_(blpsf.kappa_),
+ Uref_(blpsf.Uref_),
+ Href_(blpsf.Href_),
zGround_(blpsf.zGround_)
{}
@@ -138,14 +151,16 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::updateCoeffs()
void atmBoundaryLayerInletEpsilonFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
- os.writeKeyword("Ustar")
- << Ustar_ << token::END_STATEMENT << nl;
os.writeKeyword("z")
<< z_ << token::END_STATEMENT << nl;
os.writeKeyword("z0")
<< z0_ << token::END_STATEMENT << nl;
os.writeKeyword("kappa")
<< kappa_ << token::END_STATEMENT << nl;
+ os.writeKeyword("Uref")
+ << Uref_ << token::END_STATEMENT << nl;
+ os.writeKeyword("Href")
+ << Href_ << token::END_STATEMENT << nl;
os.writeKeyword("zGround")
<< zGround_ << token::END_STATEMENT << nl;
writeEntry("value", os);
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
index ab08c621c95da1aa092c751da6192598e066a7d8..a1ff656334b9c5f17a9b0eefafc35615008fb6bc 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H
@@ -41,6 +41,16 @@ Description
z0 is the surface roughness lenght
zGround minium vlaue in z direction
+
+ and:
+
+ Ustar = K Uref/ln((Zref + z0)/z0)
+
+ where:
+
+ Uref is the reference velocity at Zref
+ Zref is the reference height.
+
\endverbatim
Reference:
@@ -78,17 +88,23 @@ class atmBoundaryLayerInletEpsilonFvPatchScalarField
// Private data
//- Frictional velocity
- const scalar Ustar_;
+ scalarField Ustar_;
//- Direction of the z-coordinate
vector z_;
//- Surface roughness lenght
- const scalar z0_;
+ scalarField z0_;
//- Von Karman constant
const scalar kappa_;
+ //- Reference velocity
+ const scalar Uref_;
+
+ //- Reference hight
+ const scalar Href_;
+
//- Minimum corrdinate value in z direction
const scalarField zGround_;
@@ -158,7 +174,7 @@ public:
// Member functions
//- Return max value
- scalar Ustar() const
+ scalarField Ustar() const
{
return Ustar_;
}
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
index d2c3dfde62d47a8570b63b1ed07482463a381e92..628b300afd389db3bda931085dae6cbcc68b6089 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C
@@ -88,16 +88,16 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
)
:
fixedValueFvPatchVectorField(p, iF),
- Ustar_(0),
+ Ustar_(p.size()),
n_(dict.lookup("n")),
z_(dict.lookup("z")),
- z0_(readScalar(dict.lookup("z0"))),
+ z0_("z0", dict, p.size()),
kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
Uref_(readScalar(dict.lookup("Uref"))),
Href_(readScalar(dict.lookup("Href"))),
zGround_("zGround", dict, p.size())
{
- if (mag(n_) < SMALL || mag(z_) < SMALL || mag(z0_) < SMALL)
+ if (mag(n_) < SMALL || mag(z_) < SMALL)
{
FatalErrorIn
(
@@ -108,14 +108,17 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
"onst dictionary&"
")"
)
- << "magnitude of n, z and z0 vectors must be greater than zero"
+ << "magnitude of n or z must be greater than zero"
<< abort(FatalError);
}
n_ /= mag(n_);
z_ /= mag(z_);
- Ustar_ = kappa_*Uref_/(log((Href_ + z0_)/max(z0_ , 0.001)));
+ forAll (Ustar_, i)
+ {
+ Ustar_[i] = kappa_*Uref_/(log((Href_ + z0_[i])/max(z0_[i] , 0.001)));
+ }
evaluate();
}
@@ -153,12 +156,12 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
if ((coord[i] - zGround_[i]) < Href_)
{
Un[i] =
- (Ustar_/kappa_)
- * log((coord[i] - zGround_[i] + z0_)/max(z0_, 0.001));
+ (Ustar_[i]/kappa_)
+ * log((coord[i] - zGround_[i] + z0_[i])/max(z0_[i], 0.001));
}
else
{
- Un[i] = (Uref_);
+ Un[i] = Uref_;
}
}
@@ -171,8 +174,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
{
fvPatchVectorField::write(os);
- os.writeKeyword("z0")
- << z0_ << token::END_STATEMENT << nl;
+ zGround_.writeEntry("z0", os) ;
os.writeKeyword("n")
<< n_ << token::END_STATEMENT << nl;
os.writeKeyword("z")
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
index a3a87ab9d2931cea7413a3df5bda2793e67466ce..1cfd6616c73f79942318a4c74907c00363728bbd 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H
@@ -92,7 +92,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField
// Private data
//- Frictional velocity
- scalar Ustar_;
+ scalarField Ustar_;
//- Flow direction
vector n_;
@@ -101,7 +101,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField
vector z_;
//- Surface roughness lenght
- const scalar z0_;
+ scalarField z0_;
//- Von Karman constant
const scalar kappa_;
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
index 2ec45c7f60bdbfb3b1dd9832fbfc4bc5110a820d..162de9255a0d560f93176252f8c25468dbd59ac8 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C
@@ -65,22 +65,10 @@ tmp<scalarField> nutkAtmRoughWallFunctionFvPatchScalarField::calcNut() const
scalar uStar = Cmu25*sqrt(k[faceCellI]);
scalar yPlus = uStar*y[faceI]/nuw[faceI];
- scalar Edash = (y[faceI] + z0_[faceI] - zGround_[faceI])/z0_[faceI];
+ scalar Edash = (y[faceI] + z0_[faceI])/z0_[faceI];
- scalar limitingNutw = max(nutw[faceI], nuw[faceI]);
-
- // To avoid oscillations limit the change in the wall viscosity
- // which is particularly important if it temporarily becomes zero
nutw[faceI] =
- max
- (
- min
- (
- nuw[faceI]
- *(yPlus*kappa_/log(max(Edash, 1+1e-4)) - 1),
- 2*limitingNutw
- ), 0.5*limitingNutw
- );
+ nuw[faceI]*(yPlus*kappa_/log(max(Edash, 1+1e-4)) - 1);
if (debug)
{
@@ -105,8 +93,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
)
:
nutkWallFunctionFvPatchScalarField(p, iF),
- z0_(p.size(), 0.0),
- zGround_(p.size(), 0.0)
+ z0_(p.size(), 0.0)
{}
@@ -120,8 +107,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
)
:
nutkWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
- z0_(ptf.z0_, mapper),
- zGround_(ptf.zGround_, mapper)
+ z0_(ptf.z0_, mapper)
{}
@@ -134,8 +120,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
)
:
nutkWallFunctionFvPatchScalarField(p, iF, dict),
- z0_("z0", dict, p.size()),
- zGround_("zGround", dict, p.size())
+ z0_("z0", dict, p.size())
{}
@@ -146,8 +131,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
)
:
nutkWallFunctionFvPatchScalarField(rwfpsf),
- z0_(rwfpsf.z0_),
- zGround_(rwfpsf.zGround_)
+ z0_(rwfpsf.z0_)
{}
@@ -159,8 +143,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
)
:
nutkWallFunctionFvPatchScalarField(rwfpsf, iF),
- z0_(rwfpsf.z0_),
- zGround_(rwfpsf.zGround_)
+ z0_(rwfpsf.z0_)
{}
@@ -173,7 +156,6 @@ void nutkAtmRoughWallFunctionFvPatchScalarField::autoMap
{
nutkWallFunctionFvPatchScalarField::autoMap(m);
z0_.autoMap(m);
- zGround_.autoMap(m);
}
@@ -189,7 +171,6 @@ void nutkAtmRoughWallFunctionFvPatchScalarField::rmap
refCast<const nutkAtmRoughWallFunctionFvPatchScalarField>(ptf);
z0_.rmap(nrwfpsf.z0_, addr);
- zGround_.rmap(nrwfpsf.zGround_, addr);
}
@@ -198,7 +179,6 @@ void nutkAtmRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
fvPatchField<scalar>::write(os);
writeLocalEntries(os);
z0_.writeEntry("z0", os);
- zGround_.writeEntry("zGround", os);
writeEntry("value", os);
}
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.H
index 803985ced7e879c4d06f21a64f7be01f98e8a332..aa045de4234d019479a1cbb1ebd2ccc4eca8e575 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.H
@@ -74,9 +74,6 @@ protected:
//- Surface roughness lenght
scalarField z0_;
- //- Minimum corrdinate value in z direction
- scalarField zGround_;
-
// Protected Member Functions