rotorDiskSource.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 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 "trimModel.H"
28 #include "fvMatrices.H"
29 #include "geometricOneField.H"
30 #include "syncTools.H"
31 #include "axesRotation.H"
33 
34 using namespace Foam::constant;
35 
36 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40  namespace fv
41  {
44  }
45 
46  template<> const char*
48  {
49  "auto",
50  "specified"
51  };
52 
55 
56  template<> const char*
58  {
59  "fixed",
60  "surfaceNormal",
61  "local"
62  };
63 
66 }
67 
68 
69 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
70 
71 void Foam::fv::rotorDiskSource::readCoeffs()
72 {
73  UName_ = coeffs().lookupOrDefault<word>("U", "U");
74 
75  // Read co-ordinate system/geometry invariant properties
76  scalar rpm(coeffs().lookup<scalar>("rpm"));
77  omega_ = rpm/60.0*mathematical::twoPi;
78 
79  coeffs().lookup("nBlades") >> nBlades_;
80 
81  inletFlow_ = inletFlowTypeNames_.read(coeffs().lookup("inletFlowType"));
82 
83  coeffs().lookup("tipEffect") >> tipEffect_;
84 
85  const dictionary& flapCoeffs(coeffs().subDict("flapCoeffs"));
86  flapCoeffs.lookup("beta0") >> flap_.beta0;
87  flapCoeffs.lookup("beta1c") >> flap_.beta1c;
88  flapCoeffs.lookup("beta2s") >> flap_.beta2s;
89  flap_.beta0 = degToRad(flap_.beta0);
90  flap_.beta1c = degToRad(flap_.beta1c);
91  flap_.beta2s = degToRad(flap_.beta2s);
92 
93  // Create co-ordinate system
94  createCoordinateSystem();
95 
96  // Read co-odinate system dependent properties
97  checkData();
98 
99  constructGeometry();
100 
101  trim_->read(coeffs());
102 
103  if (debug)
104  {
105  writeField("thetag", trim_->thetag()(), true);
106  writeField("faceArea", area_, true);
107  }
108 }
109 
110 
111 void Foam::fv::rotorDiskSource::checkData()
112 {
113  // Set inflow type
114  switch (set_.selectionType())
115  {
119  {
120  // Set the profile ID for each blade section
121  profiles_.connectBlades(blade_.profileName(), blade_.profileID());
122  switch (inletFlow_)
123  {
125  {
126  coeffs().lookup("inletVelocity") >> inletVelocity_;
127  break;
128  }
129  case inletFlowType::surfaceNormal:
130  {
131  scalar UIn
132  (
133  coeffs().lookup<scalar>("inletNormalVelocity")
134  );
135  inletVelocity_ = -coordSys_.R().e3()*UIn;
136  break;
137  }
138  case inletFlowType::local:
139  {
140  break;
141  }
142  default:
143  {
145  << "Unknown inlet velocity type" << abort(FatalError);
146  }
147  }
148  break;
149  }
150  default:
151  {
153  << "Source cannot be used with '"
154  << fvCellSet::selectionTypeNames[set_.selectionType()]
155  << "' mode. Please use one of: " << nl
162  << exit(FatalError);
163  }
164  }
165 }
166 
167 
168 void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
169 {
170  area_ = 0.0;
171 
172  static const scalar tol = 0.8;
173 
174  const label nInternalFaces = mesh().nInternalFaces();
175  const polyBoundaryMesh& pbm = mesh().boundaryMesh();
176  const vectorField& Sf = mesh().Sf();
177  const scalarField& magSf = mesh().magSf();
178 
179  vector n = Zero;
180 
181  // Calculate cell addressing for selected cells
182  labelList cellAddr(mesh().nCells(), -1);
183  UIndirectList<label>(cellAddr, set_.cells()) =
184  identityMap(set_.nCells());
185  labelList nbrFaceCellAddr(mesh().nFaces() - nInternalFaces, -1);
186  forAll(pbm, patchi)
187  {
188  const polyPatch& pp = pbm[patchi];
189 
190  if (pp.coupled())
191  {
192  forAll(pp, i)
193  {
194  label facei = pp.start() + i;
195  label nbrFacei = facei - nInternalFaces;
196  label own = mesh().faceOwner()[facei];
197  nbrFaceCellAddr[nbrFacei] = cellAddr[own];
198  }
199  }
200  }
201 
202  // Correct for parallel running
203  syncTools::swapBoundaryFaceList(mesh(), nbrFaceCellAddr);
204 
205  // Add internal field contributions
206  for (label facei = 0; facei < nInternalFaces; facei++)
207  {
208  const label own = cellAddr[mesh().faceOwner()[facei]];
209  const label nbr = cellAddr[mesh().faceNeighbour()[facei]];
210 
211  if ((own != -1) && (nbr == -1))
212  {
213  vector nf = Sf[facei]/magSf[facei];
214 
215  if ((nf & axis) > tol)
216  {
217  area_[own] += magSf[facei];
218  n += Sf[facei];
219  }
220  }
221  else if ((own == -1) && (nbr != -1))
222  {
223  vector nf = Sf[facei]/magSf[facei];
224 
225  if ((-nf & axis) > tol)
226  {
227  area_[nbr] += magSf[facei];
228  n -= Sf[facei];
229  }
230  }
231  }
232 
233 
234  // Add boundary contributions
235  forAll(pbm, patchi)
236  {
237  const polyPatch& pp = pbm[patchi];
238  const vectorField& Sfp = mesh().Sf().boundaryField()[patchi];
239  const scalarField& magSfp = mesh().magSf().boundaryField()[patchi];
240 
241  if (pp.coupled())
242  {
243  forAll(pp, j)
244  {
245  const label facei = pp.start() + j;
246  const label own = cellAddr[mesh().faceOwner()[facei]];
247  const label nbr = nbrFaceCellAddr[facei - nInternalFaces];
248  const vector nf = Sfp[j]/magSfp[j];
249 
250  if ((own != -1) && (nbr == -1) && ((nf & axis) > tol))
251  {
252  area_[own] += magSfp[j];
253  n += Sfp[j];
254  }
255  }
256  }
257  else
258  {
259  forAll(pp, j)
260  {
261  const label facei = pp.start() + j;
262  const label own = cellAddr[mesh().faceOwner()[facei]];
263  const vector nf = Sfp[j]/magSfp[j];
264 
265  if ((own != -1) && ((nf & axis) > tol))
266  {
267  area_[own] += magSfp[j];
268  n += Sfp[j];
269  }
270  }
271  }
272  }
273 
274  if (correct)
275  {
276  reduce(n, sumOp<vector>());
277  axis = n/mag(n);
278  }
279 
280  if (debug)
281  {
282  volScalarField area
283  (
284  IOobject
285  (
286  name() + ":area",
287  mesh().time().name(),
288  mesh(),
291  ),
292  mesh(),
294  );
295  UIndirectList<scalar>(area.primitiveField(), set_.cells()) = area_;
296 
297  Info<< type() << ": " << name() << " writing field " << area.name()
298  << endl;
299 
300  area.write();
301  }
302 }
303 
304 
305 void Foam::fv::rotorDiskSource::createCoordinateSystem()
306 {
307  // Construct the local rotor co-prdinate system
308  vector origin(Zero);
309  vector axis(Zero);
310  vector refDir(Zero);
311 
312  geometryModeType gm =
313  geometryModeTypeNames_.read(coeffs().lookup("geometryMode"));
314 
315  switch (gm)
316  {
317  case geometryModeType::automatic:
318  {
319  // Determine rotation origin (cell volume weighted)
320  scalar sumV = 0.0;
321  const scalarField& V = mesh().V();
322  const vectorField& C = mesh().C();
323 
324  const labelUList cells = set_.cells();
325 
326  forAll(cells, i)
327  {
328  const label celli = cells[i];
329  sumV += V[celli];
330  origin += V[celli]*C[celli];
331  }
332  reduce(origin, sumOp<vector>());
333  reduce(sumV, sumOp<scalar>());
334  origin /= sumV;
335 
336  // Determine first radial vector
337  vector dx1(Zero);
338  scalar magR = -great;
339  forAll(cells, i)
340  {
341  const label celli = cells[i];
342  vector test = C[celli] - origin;
343  if (mag(test) > magR)
344  {
345  dx1 = test;
346  magR = mag(test);
347  }
348  }
349  reduce(dx1, maxMagSqrOp<vector>());
350  magR = mag(dx1);
351 
352  // Determine second radial vector and cross to determine axis
353  forAll(cells, i)
354  {
355  const label celli = cells[i];
356  vector dx2 = C[celli] - origin;
357  if (mag(dx2) > 0.5*magR)
358  {
359  axis = dx1 ^ dx2;
360  if (mag(axis) > small)
361  {
362  break;
363  }
364  }
365  }
366  reduce(axis, maxMagSqrOp<vector>());
367  axis /= mag(axis);
368 
369  // Correct the axis direction using a point above the rotor
370  {
371  vector pointAbove(coeffs().lookup("pointAbove"));
372  vector dir = pointAbove - origin;
373  dir /= mag(dir);
374  if ((dir & axis) < 0)
375  {
376  axis *= -1.0;
377  }
378  }
379 
380  coeffs().lookup("refDirection") >> refDir;
381 
382  cylindrical_.reset
383  (
384  new cylindrical(axis, origin, UIndirectList<vector>(C, cells)())
385  );
386 
387  // Set the face areas and apply correction to calculated axis
388  // e.g. if cellZone is more than a single layer in thickness
389  setFaceArea(axis, true);
390 
391  break;
392  }
393  case geometryModeType::specified:
394  {
395  coeffs().lookup("origin") >> origin;
396  coeffs().lookup("axis") >> axis;
397  coeffs().lookup("refDirection") >> refDir;
398 
399  cylindrical_.reset
400  (
401  new cylindrical
402  (
403  axis,
404  origin,
405  UIndirectList<vector>(mesh().C(), set_.cells())()
406  )
407  );
408 
409  setFaceArea(axis, false);
410 
411  break;
412  }
413  }
414 
415  coordSys_ = coordinateSystems::cylindrical
416  (
417  "rotorCoordSys",
418  origin,
419  axis,
420  refDir,
421  false
422  );
423 
424  const scalar sumArea = gSum(area_);
425  const scalar diameter = Foam::sqrt(4.0*sumArea/mathematical::pi);
426  Info<< " Rotor geometry:" << nl
427  << " - disk diameter = " << diameter << nl
428  << " - disk area = " << sumArea << nl
429  << " - origin = " << coordSys_.origin() << nl
430  << " - r-axis = " << coordSys_.R().e1() << nl
431  << " - psi-axis = " << coordSys_.R().e2() << nl
432  << " - z-axis = " << coordSys_.R().e3() << endl;
433 }
434 
435 
436 void Foam::fv::rotorDiskSource::constructGeometry()
437 {
438  const vectorField& C = mesh().C();
439 
440  const labelUList cells = set_.cells();
441 
442  forAll(cells, i)
443  {
444  if (area_[i] > rootVSmall)
445  {
446  const label celli = cells[i];
447 
448  // Position in (planar) rotor co-ordinate system
449  x_[i] = coordSys_.localPosition(C[celli]);
450 
451  // Cache max radius
452  rMax_ = max(rMax_, x_[i].x());
453 
454  // Swept angle relative to rDir axis [radians] in range 0 -> 2*pi
455  scalar psi = x_[i].y();
456 
457  // Blade flap angle [radians]
458  scalar beta =
459  flap_.beta0 - flap_.beta1c*cos(psi) - flap_.beta2s*sin(psi);
460 
461  // Determine rotation tensor to convert from planar system into the
462  // rotor cone system
463  scalar c = cos(beta);
464  scalar s = sin(beta);
465  R_[i] = tensor(c, 0, -s, 0, 1, 0, s, 0, c);
466  invR_[i] = R_[i].T();
467  }
468  }
469 }
470 
471 
472 Foam::tmp<Foam::vectorField> Foam::fv::rotorDiskSource::inflowVelocity
473 (
474  const volVectorField& U
475 ) const
476 {
477  switch (inletFlow_)
478  {
480  case inletFlowType::surfaceNormal:
481  {
482  return tmp<vectorField>
483  (
484  new vectorField(mesh().nCells(), inletVelocity_)
485  );
486 
487  break;
488  }
489  case inletFlowType::local:
490  {
491  return U.primitiveField();
492 
493  break;
494  }
495  default:
496  {
498  << "Unknown inlet flow specification" << abort(FatalError);
499  }
500  }
501 
502  return tmp<vectorField>(new vectorField(mesh().nCells(), Zero));
503 }
504 
505 
506 // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
507 
509 (
510  const word& name,
511  const word& modelType,
512  const fvMesh& mesh,
513  const dictionary& dict
514 )
515 :
516  fvModel(name, modelType, mesh, dict),
517  set_(mesh, coeffs()),
518  UName_(word::null),
519  omega_(0),
520  nBlades_(0),
521  inletFlow_(inletFlowType::local),
522  inletVelocity_(Zero),
523  tipEffect_(1),
524  flap_(),
525  x_(set_.nCells(), Zero),
526  R_(set_.nCells(), I),
527  invR_(set_.nCells(), I),
528  area_(set_.nCells(), Zero),
529  coordSys_("rotorCoordSys", vector::zero, axesRotation(sphericalTensor::I)),
530  cylindrical_(),
531  rMax_(0),
532  trim_(trimModel::New(*this, coeffs())),
533  blade_(coeffs().subDict("blade")),
534  profiles_(coeffs().subDict("profiles")),
535  rhoRef_(1)
536 {
537  readCoeffs();
538 }
539 
540 
541 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
542 
544 {}
545 
546 
547 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
548 
550 {
551  return wordList(1, UName_);
552 }
553 
554 
556 (
557  fvMatrix<vector>& eqn,
558  const word& fieldName
559 ) const
560 {
561  volVectorField force
562  (
563  IOobject
564  (
565  name() + ":rotorForce",
566  mesh().time().name(),
567  mesh()
568  ),
569  mesh(),
571  (
572  "zero",
573  eqn.dimensions()/dimVolume,
574  Zero
575  )
576  );
577 
578  // Read the reference density for incompressible flow
579  coeffs().lookup("rhoRef") >> rhoRef_;
580 
581  const vectorField Uin(inflowVelocity(eqn.psi()));
582  trim_->correct(Uin, force);
583  calculate(geometricOneField(), Uin, trim_->thetag(), force);
584 
585  // Add source to rhs of eqn
586  eqn -= force;
587 
588  if (mesh().time().writeTime())
589  {
590  force.write();
591  }
592 }
593 
594 
596 (
597  const volScalarField& rho,
598  fvMatrix<vector>& eqn,
599  const word& fieldName
600 ) const
601 {
602  volVectorField force
603  (
604  IOobject
605  (
606  name() + ":rotorForce",
607  mesh().time().name(),
608  mesh()
609  ),
610  mesh(),
612  (
613  "zero",
614  eqn.dimensions()/dimVolume,
615  Zero
616  )
617  );
618 
619  const vectorField Uin(inflowVelocity(eqn.psi()));
620  trim_->correct(rho, Uin, force);
621  calculate(rho, Uin, trim_->thetag(), force);
622 
623  // Add source to rhs of eqn
624  eqn -= force;
625 
626  if (mesh().time().writeTime())
627  {
628  force.write();
629  }
630 }
631 
632 
634 {
635  set_.movePoints();
636  return true;
637 }
638 
639 
641 {
642  set_.topoChange(map);
643 }
644 
645 
647 {
648  set_.mapMesh(map);
649 }
650 
651 
653 {
654  set_.distribute(map);
655 }
656 
657 
659 {
660  if (fvModel::read(dict))
661  {
662  set_.read(coeffs());
663  readCoeffs();
664  return true;
665  }
666  else
667  {
668  return false;
669  }
670 }
671 
672 
673 // ************************************************************************* //
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
label n
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Generic GeometricField class.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:54
virtual Ostream & write(const char)=0
Write character.
Templated 3D SphericalTensor derived from VectorSpace adding construction from 1 component,...
A coordinate rotation specified using global axis.
Definition: axesRotation.H:67
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
VolField< Type > & psi()
Definition: fvMatrix.H:289
const dimensionSet & dimensions() const
Definition: fvMatrix.H:302
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
Finite volume model abstract base class.
Definition: fvModel.H:59
virtual bool read(const dictionary &dict)
Read source dictionary.
Definition: fvModel.C:187
Cell based momentum source which approximates the mean effects of rotor forces on a cylindrical regio...
virtual bool movePoints()
Update for mesh motion.
virtual wordList addSupFields() const
Return the list of fields for which the fvModel adds source term.
static const NamedEnum< inletFlowType, 3 > inletFlowTypeNames_
rotorDiskSource(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict)
Construct from components.
virtual void topoChange(const polyTopoChangeMap &)
Update topology using the given map.
virtual void distribute(const polyDistributionMap &)
Redistribute or update using the given distribution map.
virtual bool read(const dictionary &dict)
Read source dictionary.
virtual void addSup(fvMatrix< vector > &eqn, const word &fieldName) const
Source term to momentum equation.
virtual void mapMesh(const polyMeshMap &)
Update from another mesh using the given map.
virtual ~rotorDiskSource()
Destructor.
static const NamedEnum< geometryModeType, 2 > geometryModeTypeNames_
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
static const NamedEnum< selectionTypes, 4 > selectionTypeNames
Word list of selection type names.
Definition: polyCellSet.H:97
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
Class containing mesh-to-mesh mapping information.
Definition: polyMeshMap.H:51
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
virtual bool write(const bool write=true) const
Write using setting from DB.
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:436
A class for managing temporary objects.
Definition: tmp.H:55
Trim model base class.
Definition: trimModel.H:51
A class for handling words, derived from string.
Definition: word.H:62
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:50
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
A special matrix type and solver, designed for finite volume solutions of scalar equations.
label patchi
const cellShapeList & cells
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const volScalarField & psi
U
Definition: pEqn.H:72
void correct(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiCorr, const SpType &Sp, const SuType &Su)
const scalar twoPi(2 *pi)
const dimensionedScalar c
Speed of light in a vacuum.
Collection of constants.
label calculate(const fvMesh &mesh, const labelHashSet &patchIDs, const scalar minFaceFraction, GeometricField< scalar, PatchField, GeoMesh > &distance)
Calculate distance data from patches.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
static const zero Zero
Definition: zero.H:97
List< word > wordList
A List of words.
Definition: fileName.H:54
Type gSum(const FieldField< Field, Type > &f)
VolField< vector > volVectorField
Definition: volFieldsFwd.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensionedScalar sin(const dimensionedScalar &ds)
IOstream & fixed(IOstream &io)
Definition: IOstream.H:573
messageStream Info
static const Identity< scalar > I
Definition: Identity.H:93
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
const dimensionSet dimVolume
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:61
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Field< vector > vectorField
Specialisation of Field<T> for vector.
error FatalError
const dimensionSet dimArea
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
UList< label > labelUList
Definition: UList.H:65
static const char nl
Definition: Ostream.H:260
dimensionedScalar cos(const dimensionedScalar &ds)
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
labelList fv(nPoints)
dictionary dict