diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
index 06e93e4901ffd259ff8f6f770f17062fd3faad56..407de61b3988746ec85d797ea8651d0fdb6e2b19 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
@@ -32,7 +32,7 @@ using namespace Foam::vectorTools;
void Foam::conformalVoronoiMesh::insertBoundingPoints()
{
- pointField farPts = geometryToConformTo_.globalBounds().points();
+ pointField farPts(geometryToConformTo_.globalBounds().points());
// Shift corners of bounds relative to origin
farPts -= geometryToConformTo_.globalBounds().midpoint();
diff --git a/src/fieldSources/Make/files b/src/fieldSources/Make/files
index f642dd2b65a2e50f0760a1af10b807395346dab6..eced302bfded50b2988169b216e37a6ac27d068a 100644
--- a/src/fieldSources/Make/files
+++ b/src/fieldSources/Make/files
@@ -18,7 +18,7 @@ basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C
basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C
basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
-basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
+basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
basicSource/actuationDiskSource/actuationDiskSource.C
basicSource/radialActuationDiskSource/radialActuationDiskSource.C
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
index 6014b623fb90cb587292edc40a7772bf78f4a4c6..c3f042fe0160277a4773b46da6f4b334259e2e41 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
@@ -465,22 +465,15 @@ Foam::rotorDiskSource::~rotorDiskSource()
void Foam::rotorDiskSource::calculate
(
const vectorField& U,
- const scalarField& alphag,
+ const scalarField& thetag,
vectorField& force,
const bool divideVolume,
const bool output
) const
{
const scalarField& V = mesh_.V();
- const bool compressible = rhoName_ != "none";
-
- tmp<volScalarField> trho
- (
- compressible
- ? mesh_.lookupObject<volScalarField>(rhoName_)
- : volScalarField::null()
- );
-
+ const bool compressible = this->compressible();
+ tmp<volScalarField> trho(rho());
// logging info
scalar dragEff = 0.0;
@@ -518,7 +511,7 @@ void Foam::rotorDiskSource::calculate
blade_.interpolate(radius, twist, chord, i1, i2, invDr);
// flip geometric angle if blade is spinning in reverse (clockwise)
- scalar alphaGeom = alphag[i] + twist;
+ scalar alphaGeom = thetag[i] + twist;
if (omega_ < 0)
{
alphaGeom = mathematical::pi - alphaGeom;
@@ -685,7 +678,7 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
if (debug)
{
- writeField("alphag", trim_->thetag()(), true);
+ writeField("thetag", trim_->thetag()(), true);
writeField("faceArea", area_, true);
}
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
index a2f856988da5ba9418dd264486941ce5473676b8..7fd7d33964602b8574da194ff178eb882b546445 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
+++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
@@ -253,7 +253,11 @@ public:
// Member Functions
// Access
-
+
+ //- Return the rotational speed [rad/s]
+ // Positive anti-clockwise when looking along -ve lift direction
+ inline scalar omega() const;
+
//- Return the cell centre positions in local rotor frame
// (Cylindrical r, theta, z)
inline const List<point>& x() const;
@@ -261,6 +265,12 @@ public:
//- Return the rotor co-ordinate system (r, theta, z)
inline const cylindricalCS& coordSys() const;
+ //- Return true if solving a compressible case
+ inline bool compressible() const;
+
+ //- Return the density field [kg/m3]
+ inline tmp<volScalarField> rho() const;
+
// Evaluation
@@ -268,7 +278,7 @@ public:
void calculate
(
const vectorField& U,
- const scalarField& alphag,
+ const scalarField& thetag,
vectorField& force,
const bool divideVolume = true,
const bool output = true
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H
index 23ca3b0e80d6aa258c6db2d998385e824ec21b84..7c194e8713e4beb538292f9e2fd3a6f1c44a9ed5 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H
+++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceI.H
@@ -25,6 +25,13 @@ License
#include "rotorDiskSource.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::rotorDiskSource::omega() const
+{
+ return omega_;
+}
+
const Foam::List<Foam::point>& Foam::rotorDiskSource::x() const
{
return x_;
@@ -37,5 +44,23 @@ const Foam::cylindricalCS& Foam::rotorDiskSource::coordSys() const
}
+bool Foam::rotorDiskSource::compressible() const
+{
+ return rhoName_ != "none";
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::rotorDiskSource::rho() const
+{
+ if (compressible())
+ {
+ return mesh_.lookupObject<volScalarField>(rhoName_);
+ }
+ else
+ {
+ return volScalarField::null();
+ }
+}
+
// ************************************************************************* //
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
similarity index 67%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
rename to src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
index f22226f84c1cab0adfa0d0c5562ccdcf59a91249..d53d73f3ac454741e6c81556f8fc29bbdcb97df4 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.C
@@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
-#include "targetForceTrim.H"
+#include "targetCoeffTrim.H"
#include "addToRunTimeSelectionTable.H"
#include "unitConversion.H"
#include "mathematicalConstants.H"
@@ -34,15 +34,15 @@ using namespace Foam::constant;
namespace Foam
{
- defineTypeNameAndDebug(targetForceTrim, 0);
+ defineTypeNameAndDebug(targetCoeffTrim, 0);
- addToRunTimeSelectionTable(trimModel, targetForceTrim, dictionary);
+ addToRunTimeSelectionTable(trimModel, targetCoeffTrim, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
-Foam::vector Foam::targetForceTrim::calcForce
+Foam::vector Foam::targetCoeffTrim::calcCoeffs
(
const vectorField& U,
const scalarField& thetag,
@@ -51,34 +51,51 @@ Foam::vector Foam::targetForceTrim::calcForce
{
rotor_.calculate(U, thetag, force, false, false);
+ bool compressible = rotor_.compressible();
+ tmp<volScalarField> trho = rotor_.rho();
+
const labelList& cells = rotor_.cells();
const vectorField& C = rotor_.mesh().C();
+ const List<point>& x = rotor_.x();
const vector& origin = rotor_.coordSys().origin();
const vector& rollAxis = rotor_.coordSys().e1();
const vector& pitchAxis = rotor_.coordSys().e2();
const vector& yawAxis = rotor_.coordSys().e3();
- vector f(vector::zero);
+ scalar coeff1 = alpha_*sqr(rotor_.omega())*mathematical::pi;
+
+ vector cf(vector::zero);
forAll(cells, i)
{
label cellI = cells[i];
- vector moment = force[cellI]^(C[cellI] - origin);
- f[0] += force[cellI] & yawAxis;
- f[1] += moment & pitchAxis;
- f[2] += moment & rollAxis;
+ // create normalisation coefficient
+ scalar radius = x[i].x();
+ scalar coeff2 = coeff1*pow4(radius);
+ if (compressible)
+ {
+ coeff2 *= trho()[cellI];
+ }
+
+ vector fc = force[cellI];
+ vector mc = fc^(C[cellI] - origin);
+
+ // add to coefficient vector
+ cf[0] += (fc & yawAxis)/(coeff2 + ROOTVSMALL);
+ cf[1] += (mc & pitchAxis)/(coeff2*radius + ROOTVSMALL);
+ cf[2] += (mc & rollAxis)/(coeff2*radius + ROOTVSMALL);
}
- reduce(f, sumOp<vector>());
+ reduce(cf, sumOp<vector>());
- return f;
+ return cf;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-Foam::targetForceTrim::targetForceTrim
+Foam::targetCoeffTrim::targetCoeffTrim
(
const rotorDiskSource& rotor,
const dictionary& dict
@@ -91,7 +108,8 @@ Foam::targetForceTrim::targetForceTrim
nIter_(50),
tol_(1e-8),
relax_(1.0),
- dTheta_(degToRad(0.1))
+ dTheta_(degToRad(0.1)),
+ alpha_(1.0)
{
read(dict);
}
@@ -99,20 +117,20 @@ Foam::targetForceTrim::targetForceTrim
// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * //
-Foam::targetForceTrim::~targetForceTrim()
+Foam::targetCoeffTrim::~targetCoeffTrim()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-void Foam::targetForceTrim::read(const dictionary& dict)
+void Foam::targetCoeffTrim::read(const dictionary& dict)
{
trimModel::read(dict);
const dictionary& targetDict(coeffs_.subDict("target"));
- target_[0] = readScalar(targetDict.lookup("fThrust"));
- target_[1] = readScalar(targetDict.lookup("mPitch"));
- target_[2] = readScalar(targetDict.lookup("mRoll"));
+ target_[0] = readScalar(targetDict.lookup("thrustCoeff"));
+ target_[1] = readScalar(targetDict.lookup("pitchCoeff"));
+ target_[2] = readScalar(targetDict.lookup("rollCoeff"));
const dictionary& pitchAngleDict(coeffs_.subDict("pitchAngles"));
theta_[0] = degToRad(readScalar(pitchAngleDict.lookup("theta0Ini")));
@@ -129,10 +147,12 @@ void Foam::targetForceTrim::read(const dictionary& dict)
{
dTheta_ = degToRad(dTheta_);
}
+
+ alpha_ = readScalar(coeffs_.lookup("alpha"));
}
-Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
+Foam::tmp<Foam::scalarField> Foam::targetCoeffTrim::thetag() const
{
const List<vector>& x = rotor_.x();
@@ -149,34 +169,38 @@ Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
}
-void Foam::targetForceTrim::correct(const vectorField& U, vectorField& force)
+void Foam::targetCoeffTrim::correct(const vectorField& U, vectorField& force)
{
if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
{
+ Info<< type() << ":" << nl
+ << " solving for target trim coefficients" << nl;
+
// iterate to find new pitch angles to achieve target force
scalar err = GREAT;
label iter = 0;
tensor J(tensor::zero);
+ vector old = vector::zero;
while ((err > tol_) && (iter < nIter_))
{
// cache initial theta vector
vector theta0(theta_);
// set initial values
- vector old = calcForce(U, thetag(), force);
+ old = calcCoeffs(U, thetag(), force);
// construct Jacobian by perturbing the pitch angles
// by +/-(dTheta_/2)
for (label pitchI = 0; pitchI < 3; pitchI++)
{
theta_[pitchI] -= dTheta_/2.0;
- vector f0 = calcForce(U, thetag(), force);
+ vector cf0 = calcCoeffs(U, thetag(), force);
theta_[pitchI] += dTheta_;
- vector f1 = calcForce(U, thetag(), force);
+ vector cf1 = calcCoeffs(U, thetag(), force);
- vector ddTheta = (f1 - f0)/dTheta_;
+ vector ddTheta = (cf1 - cf0)/dTheta_;
J[pitchI + 0] = ddTheta[0];
J[pitchI + 3] = ddTheta[1];
@@ -201,23 +225,21 @@ void Foam::targetForceTrim::correct(const vectorField& U, vectorField& force)
if (iter == nIter_)
{
- WarningIn
- (
- "void Foam::targetForceTrim::correct"
- "("
- "const vectorField&, "
- "vectorField&"
- ")"
- ) << "Trim routine not converged in " << iter
- << " iterations, max residual = " << err << endl;
+ Info<< " solution not converged in " << iter
+ << " iterations, final residual = " << err
+ << "(" << tol_ << ")" << endl;
}
else
{
- Info<< type() << ": converged in " << iter
- << " iterations" << endl;
+ Info<< " final residual = " << err << "(" << tol_
+ << "), iterations = " << iter << endl;
}
- Info<< " new pitch angles:" << nl
+ Info<< " current and target coefficients:" << nl
+ << " thrust = " << old[0] << ", " << target_[0] << nl
+ << " pitch = " << old[1] << ", " << target_[1] << nl
+ << " roll = " << old[2] << ", " << target_[2] << nl
+ << " new pitch angles [deg]:" << nl
<< " theta0 = " << radToDeg(theta_[0]) << nl
<< " theta1c = " << radToDeg(theta_[1]) << nl
<< " theta1s = " << radToDeg(theta_[2]) << nl
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.H b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H
similarity index 65%
rename from src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.H
rename to src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H
index a030cfb4a915f3b3dfcd0c5268a4e0849b4d26b6..e387a4c1630fa36b14cdd80e76d59eb45bb28391 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.H
+++ b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetCoeff/targetCoeffTrim.H
@@ -22,18 +22,55 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
- Foam::targetForceTrim
+ Foam::targetCoeffTrim
Description
- Target force trim coefficients
+ Target trim coefficients
+
+ Solves:
+
+ c^old + J.d(theta) = c^target
+
+ Where:
+
+ n = time level
+ c = coefficient vector (thrust force, roll moment, pitch moment)
+ theta = pitch angle vector (collective, roll, pitch)
+ J = Jacobian [3x3] matrix
+
+
+ The trimmed pitch angles are found via solving the above with a
+ Newton-Raphson iterative method. The solver tolerance can be user-input,
+ using the 'tol' entry.
+
+
+ The coefficients are determined using:
+
+ force
+ c = ---------------------------------
+ alpha*pi*rho*(omega^2)*(radius^4)
+
+ and
+
+ moment
+ c = ---------------------------------
+ alpha*pi*rho*(omega^2)*(radius^5)
+
+ Where:
+
+ alpha = user-input conversion coefficient (default = 1)
+ rho = desity
+ omega = rotor angulr velocity
+ pi = mathematical pi
+
SourceFiles
- targetForceTrim.C
+ targetCoeffTrim.C
\*---------------------------------------------------------------------------*/
-#ifndef targetForceTrim_H
-#define targetForceTrim_H
+#ifndef targetCoeffTrim_H
+#define targetCoeffTrim_H
#include "trimModel.H"
#include "tensor.H"
@@ -45,10 +82,10 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
- Class targetForceTrim Declaration
+ Class targetCoeffTrim Declaration
\*---------------------------------------------------------------------------*/
-class targetForceTrim
+class targetCoeffTrim
:
public trimModel
{
@@ -60,7 +97,7 @@ protected:
//- Number of iterations between calls to 'correct'
label calcFrequency_;
- //- Target force [N]
+ //- Target coefficient vector (thrust force, roll moment, pitch moment)
vector target_;
//- Pitch angles (collective, roll, pitch) [rad]
@@ -78,11 +115,14 @@ protected:
//- Perturbation angle used to determine jacobian
scalar dTheta_;
+ //- Coefficient to allow for conversion between US and EU definitions
+ scalar alpha_;
+
// Protected member functions
- //- Calculate the rotor forces
- vector calcForce
+ //- Calculate the rotor force and moment coefficients vector
+ vector calcCoeffs
(
const vectorField& U,
const scalarField& alphag,
@@ -93,13 +133,13 @@ protected:
public:
//- Run-time type information
- TypeName("targetForceTrim");
+ TypeName("targetCoeffTrim");
//- Constructor
- targetForceTrim(const rotorDiskSource& rotor, const dictionary& dict);
+ targetCoeffTrim(const rotorDiskSource& rotor, const dictionary& dict);
//- Destructor
- virtual ~targetForceTrim();
+ virtual ~targetCoeffTrim();
// Member functions