rotorDiskSourceTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "rotorDiskSource.H"
27 #include "volFields.H"
28 #include "unitConversion.H"
29 
30 using namespace Foam::constant;
31 
32 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
33 
34 template<class RhoFieldType>
36 (
37  const RhoFieldType& rho,
38  const vectorField& U,
39  const scalarField& thetag,
40  vectorField& force,
41  const bool divideVolume,
42  const bool output
43 ) const
44 {
45  const scalarField& V = mesh_.V();
46 
47  // Logging info
48  scalar dragEff = 0.0;
49  scalar liftEff = 0.0;
50  scalar AOAmin = GREAT;
51  scalar AOAmax = -GREAT;
52 
53  forAll(cells_, i)
54  {
55  if (area_[i] > ROOTVSMALL)
56  {
57  const label cellI = cells_[i];
58 
59  const scalar radius = x_[i].x();
60 
61  // Transform velocity into local cylindrical reference frame
62  vector Uc = cylindrical_->invTransform(U[cellI], i);
63 
64  // Transform velocity into local coning system
65  Uc = R_[i] & Uc;
66 
67  // Set radial component of velocity to zero
68  Uc.x() = 0.0;
69 
70  // Set blade normal component of velocity
71  Uc.y() = radius*omega_ - Uc.y();
72 
73  // Determine blade data for this radius
74  // i2 = index of upper radius bound data point in blade list
75  scalar twist = 0.0;
76  scalar chord = 0.0;
77  label i1 = -1;
78  label i2 = -1;
79  scalar invDr = 0.0;
80  blade_.interpolate(radius, twist, chord, i1, i2, invDr);
81 
82  // Flip geometric angle if blade is spinning in reverse (clockwise)
83  scalar alphaGeom = thetag[i] + twist;
84  if (omega_ < 0)
85  {
86  alphaGeom = mathematical::pi - alphaGeom;
87  }
88 
89  // Effective angle of attack
90  scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
91  if (alphaEff > mathematical::pi)
92  {
93  alphaEff -= mathematical::twoPi;
94  }
95  if (alphaEff < -mathematical::pi)
96  {
97  alphaEff += mathematical::twoPi;
98  }
99 
100  AOAmin = min(AOAmin, alphaEff);
101  AOAmax = max(AOAmax, alphaEff);
102 
103  // Determine profile data for this radius and angle of attack
104  const label profile1 = blade_.profileID()[i1];
105  const label profile2 = blade_.profileID()[i2];
106 
107  scalar Cd1 = 0.0;
108  scalar Cl1 = 0.0;
109  profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);
110 
111  scalar Cd2 = 0.0;
112  scalar Cl2 = 0.0;
113  profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);
114 
115  scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
116  scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
117 
118  // Apply tip effect for blade lift
119  scalar tipFactor = neg(radius/rMax_ - tipEffect_);
120 
121  // Calculate forces perpendicular to blade
122  scalar pDyn = 0.5*rho[cellI]*magSqr(Uc);
123 
124  scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
125  vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
126 
127  // Accumulate forces
128  dragEff += rhoRef_*localForce.y();
129  liftEff += rhoRef_*localForce.z();
130 
131  // Transform force from local coning system into rotor cylindrical
132  localForce = invR_[i] & localForce;
133 
134  // Transform force into global Cartesian co-ordinate system
135  force[cellI] = cylindrical_->transform(localForce, i);
136 
137  if (divideVolume)
138  {
139  force[cellI] /= V[cellI];
140  }
141  }
142  }
143 
144  if (output)
145  {
146  reduce(AOAmin, minOp<scalar>());
147  reduce(AOAmax, maxOp<scalar>());
148  reduce(dragEff, sumOp<scalar>());
149  reduce(liftEff, sumOp<scalar>());
150 
151  Info<< type() << " output:" << nl
152  << " min/max(AOA) = " << radToDeg(AOAmin) << ", "
153  << radToDeg(AOAmax) << nl
154  << " Effective drag = " << dragEff << nl
155  << " Effective lift = " << liftEff << endl;
156  }
157 }
158 
159 
160 template<class Type>
162 (
163  const word& name,
164  const List<Type>& values,
165  const bool writeNow
166 ) const
167 {
169 
170  if (mesh_.time().outputTime() || writeNow)
171  {
172  tmp<fieldType> tfield
173  (
174  new fieldType
175  (
176  IOobject
177  (
178  name,
179  mesh_.time().timeName(),
180  mesh_,
183  ),
184  mesh_,
186  )
187  );
188 
189  Field<Type>& field = tfield().internalField();
190 
191  if (cells_.size() != values.size())
192  {
193  FatalErrorIn("") << "cells_.size() != values_.size()"
194  << abort(FatalError);
195  }
196 
197  forAll(cells_, i)
198  {
199  const label cellI = cells_[i];
200  field[cellI] = values[i];
201  }
202 
203  tfield().write();
204  }
205 }
206 
207 
208 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar neg(const dimensionedScalar &ds)
Collection of constants.
labelList f(nPoints)
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
volScalarField pDyn(IOobject( "pDyn", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedScalar("zero", dimPressure, 0.0))
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void writeField(const word &name, const List< Type > &values, const bool writeNow=false) const
Helper function to write rotor values.
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
messageStream Info
void calculate(const RhoFieldType &rho, const vectorField &U, const scalarField &thetag, vectorField &force, const bool divideVolume=true, const bool output=true) const
Calculate forces.
Generic dimensioned Type class.
static const char nl
Definition: Ostream.H:260
const Cmpt & y() const
Definition: VectorI.H:71
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
const Cmpt & x() const
Definition: VectorI.H:65
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
Unit conversion functions.
const Cmpt & z() const
Definition: VectorI.H:77
errorManip< error > abort(error &err)
Definition: errorManip.H:131
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
Generic GeometricField class.
Traits class for primitives.
Definition: pTraits.H:50
error FatalError
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const scalar twoPi(2 *pi)
A class for managing temporary objects.
Definition: PtrList.H:118
const volScalarField & alphaEff
Definition: setAlphaEff.H:49