ENH: Updates to rotor momentum source

src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,33 +47,34 @@ void Foam::bladeModel::interpolateWeights
     scalar& ddx
 ) const
 {
-    scalar x = -GREAT;
+    i2 = 0;
     label nElem = values.size();
 
-    i2 = 0;
+    if (nElem == 1)
-    while ((x < xIn) && (i2 < nElem))
     {
-        x = values[i2];
-        i2++;
-    }
-
-    if (i2 == 0)
-    {
-        i1 = i2;
-        ddx = 0.0;
-        return;
-    }
-    else if (i2 == values.size())
-    {
-        i2 = values.size() - 1;
         i1 = i2;
         ddx = 0.0;
         return;
     }
     else
     {
-        i1 = i2 - 1;
+        while ((values[i2] < xIn) && (i2 < nElem))
-        ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+        {
+            i2++;
+        }
+
+        if (i2 == nElem)
+        {
+            i2 = nElem - 1;
+            i1 = i2;
+            ddx = 0.0;
+            return;
+        }
+        else
+        {
+            i1 = i2 - 1;
+            ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+        }
     }
 }
 
@@ -120,14 +121,8 @@ Foam::bladeModel::bladeModel(const dictionary& dict)
     }
     else
     {
-        FatalErrorIn
+        FatalErrorIn("Foam::bladeModel::bladeModel(const dictionary&)")
-        (
+            << "No blade data specified" << exit(FatalError);
-            "Foam::bladeModel::bladeModel"
-            "("
-                "const dictionary&, "
-                "const word&"
-            ")"
-        )   << "No blade data specified" << exit(FatalError);
     }
 }
 

src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C

@@ -49,33 +49,34 @@ void Foam::lookupProfile::interpolateWeights
     scalar& ddx
 ) const
 {
-    scalar x = -GREAT;
+    i2 = 0;
     label nElem = values.size();
 
-    i2 = 0;
+    if (nElem == 1)
-    while ((x < xIn) && (i2 < nElem))
     {
-        x = values[i2];
-        i2++;
-    }
-
-    if (i2 == 0)
-    {
-        i1 = i2;
-        ddx = 0.0;
-        return;
-    }
-    else if (i2 == nElem)
-    {
-        i2 = nElem - 1;
         i1 = i2;
         ddx = 0.0;
         return;
     }
     else
     {
-        i1 = i2 - 1;
+        while ((values[i2] < xIn) && (i2 < nElem))
-        ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+        {
+            i2++;
+        }
+
+        if (i2 == nElem)
+        {
+            i2 = nElem - 1;
+            i1 = i2;
+            ddx = 0.0;
+            return;
+        }
+        else
+        {
+            i1 = i2 - 1;
+            ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+        }
     }
 }
 

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C

@@ -67,6 +67,7 @@ namespace Foam
 
 void Foam::rotorDiskSource::checkData()
 {
+    // set inflow type
     switch (selectionMode())
     {
         case smCellSet:
@@ -129,11 +130,10 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
 
     const label nInternalFaces = mesh_.nInternalFaces();
     const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
-    const vectorField& Sf = mesh_.Sf().internalField();
+    const vectorField& Sf = mesh_.Sf();
-    const scalarField& magSf = mesh_.magSf().internalField();
+    const scalarField& magSf = mesh_.magSf();
 
     vector n = vector::zero;
-    label nFace = 0;
 
     // calculate cell addressing for selected cells
     labelList cellAddr(mesh_.nCells(), -1);
@@ -158,7 +158,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
     // correct for parallel running
     syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
 
-
     // add internal field contributions
     for (label faceI = 0; faceI < nInternalFaces; faceI++)
     {
@@ -173,7 +172,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
             {
                 area_[own] += magSf[faceI];
                 n += Sf[faceI];
-                nFace++;
             }
         }
         else if ((own == -1) && (nbr != -1))
@@ -184,7 +182,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
             {
                 area_[nbr] += magSf[faceI];
                 n -= Sf[faceI];
-                nFace++;
             }
         }
     }
@@ -210,7 +207,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
                 {
                     area_[own] += magSfp[j];
                     n += Sfp[j];
-                    nFace++;
                 }
             }
         }
@@ -222,11 +218,10 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
                 const label own = cellAddr[mesh_.faceOwner()[faceI]];
                 const vector nf = Sfp[j]/magSfp[j];
 
-                if ((own != -1) && ((nf & axis) > 0.8))
+                if ((own != -1) && ((nf & axis) > tol))
                 {
                     area_[own] += magSfp[j];
                     n += Sfp[j];
-                    nFace++;
                 }
             }
         }
@@ -359,32 +354,30 @@ void Foam::rotorDiskSource::constructGeometry()
 
     forAll(cells_, i)
     {
-        const label cellI = cells_[i];
+        if (area_[i] > ROOTVSMALL)
-
-        // position in (planar) rotor co-ordinate system
-        x_[i] = coordSys_.localPosition(C[cellI]);
-
-        // cache max radius
-        rMax_ = max(rMax_, x_[i].x());
-
-        // swept angle relative to rDir axis [radians] in range 0 -> 2*pi
-        scalar psi = x_[i].y();
-        if (psi < 0)
         {
-            psi += mathematical::twoPi;
+            const label cellI = cells_[i];
+
+            // position in (planar) rotor co-ordinate system
+            x_[i] = coordSys_.localPosition(C[cellI]);
+
+            // cache max radius
+            rMax_ = max(rMax_, x_[i].x());
+
+            // swept angle relative to rDir axis [radians] in range 0 -> 2*pi
+            scalar psi = x_[i].y();
+
+            // blade flap angle [radians]
+            scalar beta =
+                flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi);
+
+            // determine rotation tensor to convert from planar system into the
+            // rotor cone system
+            scalar c = cos(beta);
+            scalar s = sin(beta);
+            R_[i] = tensor(c, 0, -s, 0, 1, 0, s, 0, c);
+            invR_[i] = R_[i].T();
         }
-
-        // blade flap angle [radians]
-        scalar beta = flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi);
-
-        // determine rotation tensor to convert from planar system into the
-        // rotor cone system
-        scalar cNeg = cos(-beta);
-        scalar sNeg = sin(-beta);
-        R_[i] = tensor(cNeg, 0.0, -sNeg, 0.0, 1.0, 0.0, sNeg, 0.0, cNeg);
-        scalar cPos = cos(beta);
-        scalar sPos = sin(beta);
-        invR_[i] = tensor(cPos, 0.0, -sPos, 0.0, 1.0, 0.0, sPos, 0.0, cPos);
     }
 }
 
@@ -440,6 +433,7 @@ Foam::rotorDiskSource::rotorDiskSource
 :
     basicSource(name, modelType, dict, mesh),
     rhoName_("none"),
+    rhoRef_(1.0),
     omega_(0.0),
     nBlades_(0),
     inletFlow_(ifLocal),
@@ -477,13 +471,15 @@ void Foam::rotorDiskSource::calculate
     const bool output
 ) const
 {
-    tmp<volScalarField> trho;
-    if (rhoName_ != "none")
-    {
-        trho = mesh_.lookupObject<volScalarField>(rhoName_);
-    }
-
     const scalarField& V = mesh_.V();
+    const bool compressible = rhoName_ != "none";
+
+    tmp<volScalarField> trho
+    (
+        compressible
+      ? mesh_.lookupObject<volScalarField>(rhoName_)
+      : volScalarField::null()
+    );
 
 
     // logging info
@@ -503,17 +499,17 @@ void Foam::rotorDiskSource::calculate
             // velocity in local cylindrical reference frame
             vector Uc = coordSys_.localVector(U[cellI]);
 
-            // apply correction in local system due to coning
+            // transform from rotor cylindrical into local coning system
             Uc = R_[i] & Uc;
 
             // set radial component of velocity to zero
             Uc.x() = 0.0;
 
-            // remove blade linear velocity from blade normal component
+            // set blade normal component of velocity
-            Uc.y() -= radius*omega_;
+            Uc.y() = radius*omega_ - Uc.y();
 
             // determine blade data for this radius
-            // i1 = index of upper bound data point in blade list
+            // i2 = index of upper radius bound data point in blade list
             scalar twist = 0.0;
             scalar chord = 0.0;
             label i1 = -1;
@@ -529,17 +525,7 @@ void Foam::rotorDiskSource::calculate
             }
 
             // effective angle of attack
-            scalar alphaEff =
+            scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
-                mathematical::pi + atan2(Uc.z(), Uc.y()) - alphaGeom;
-
-            if (alphaEff > mathematical::pi)
-            {
-                alphaEff -= mathematical::twoPi;
-            }
-            if (alphaEff < -mathematical::pi)
-            {
-                alphaEff += mathematical::twoPi;
-            }
 
             AOAmin = min(AOAmin, alphaEff);
             AOAmax = max(AOAmax, alphaEff);
@@ -564,21 +550,21 @@ void Foam::rotorDiskSource::calculate
 
             // calculate forces perpendicular to blade
             scalar pDyn = 0.5*magSqr(Uc);
-            if (trho.valid())
+            if (compressible)
             {
                 pDyn *= trho()[cellI];
             }
 
-            scalar f = pDyn*chord*nBlades_*area_[i]/mathematical::twoPi;
+            scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
-            vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
+            vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
+
+            // accumulate forces
+            dragEff += rhoRef_*localForce.y();
+            liftEff += rhoRef_*localForce.z();
 
             // convert force from local coning system into rotor cylindrical
             localForce = invR_[i] & localForce;
 
-            // accumulate forces
-            dragEff += localForce.y();
-            liftEff += localForce.z();
-
             // convert force to global cartesian co-ordinate system
             force[cellI] = coordSys_.globalVector(localForce);
 
@@ -616,6 +602,7 @@ void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
     }
     else
     {
+        coeffs_.lookup("rhoRef") >> rhoRef_;
         dims.reset(dimForce/dimVolume/dimDensity);
     }
 
@@ -666,6 +653,8 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
         coeffs_.lookup("fieldNames") >> fieldNames_;
         applied_.setSize(fieldNames_.size(), false);
 
+
+        // read co-ordinate system/geometry invariant properties
         scalar rpm(readScalar(coeffs_.lookup("rpm")));
         omega_ = rpm/60.0*mathematical::twoPi;
 
@@ -683,14 +672,17 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
         flap_.beta1 = degToRad(flap_.beta1);
         flap_.beta2 = degToRad(flap_.beta2);
 
-        trim_->read(coeffs_);        
-
-        checkData();
 
+        // create co-ordinate system
         createCoordinateSystem();
 
+        // read co-odinate system dependent properties
+        checkData();
+
         constructGeometry();
 
+        trim_->read(coeffs_);
+
         if (debug)
         {
             writeField("alphag", trim_->thetag()(), true);

src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H

@@ -148,6 +148,9 @@ protected:
         //- Name of density field
         word rhoName_;
 
+        //- Reference density if rhoName = 'none'
+        scalar rhoRef_;
+
         //- Rotational speed [rad/s]
         //  Positive anti-clockwise when looking along -ve lift direction
         scalar omega_;

src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C

@@ -71,11 +71,6 @@ void Foam::fixedTrim::read(const dictionary& dict)
     forAll(thetag_, i)
     {
         scalar psi = x[i].y();
-        if (psi < 0)
-        {
-            psi += mathematical::twoPi;
-        }
-
         thetag_[i] = theta0 + theta1c*cos(psi) + theta1s*sin(psi);
     }
 }

src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C

@@ -142,11 +142,6 @@ Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
     forAll(t, i)
     {
         scalar psi = x[i].y();
-        if (psi < 0)
-        {
-            psi += mathematical::twoPi;
-        }
-
         t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi);
     }