ENH: Updated rotor source trim sub-model

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

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);
+        // create normalisation coefficient
-        f[0] += force[cellI] & yawAxis;
+        scalar radius = x[i].x();
-        f[1] += moment & pitchAxis;
+        scalar coeff2 = coeff1*pow4(radius);
-        f[2] += moment & rollAxis;
+        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_[0] = readScalar(targetDict.lookup("thrustCoeff"));
-    target_[1] = readScalar(targetDict.lookup("mPitch"));
+    target_[1] = readScalar(targetDict.lookup("pitchCoeff"));
-    target_[2] = readScalar(targetDict.lookup("mRoll"));
+    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);
-
+if (iter == 0) Info<< "old=" << old << endl;
             // 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
+            Info<< "    solution not converged in " << iter
-            (
+                << " iterations, final residual = " << err
-                "void Foam::targetForceTrim::correct"
+                << "(" << tol_ << ")" << endl;
-                "("
-                    "const vectorField&, "
-                    "vectorField&"
-                ")"
-            )   << "Trim routine not converged in " << iter
-                << " iterations, max residual = " << err << endl;
         }
         else
         {
-            Info<< type() << ": converged in " << iter
+            Info<< "    final residual = " << err << "(" << tol_
-                << " iterations" << endl;
+                << "), 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

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
+#ifndef targetCoeffTrim_H
-#define targetForceTrim_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
+        //- Calculate the rotor force and moment coefficients vector
-        vector calcForce
+        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