rotorDiskSource.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-2016 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"
28 #include "trimModel.H"
29 #include "fvMatrices.H"
30 #include "geometricOneField.H"
31 #include "syncTools.H"
32 
33 using namespace Foam::constant;
34 
35 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39  namespace fv
40  {
41  defineTypeNameAndDebug(rotorDiskSource, 0);
42  addToRunTimeSelectionTable(option, rotorDiskSource, dictionary);
43  }
44 
45  template<> const char* NamedEnum<fv::rotorDiskSource::geometryModeType, 2>::
46  names[] =
47  {
48  "auto",
49  "specified"
50  };
51 
52  const NamedEnum<fv::rotorDiskSource::geometryModeType, 2>
54 
55  template<> const char* NamedEnum<fv::rotorDiskSource::inletFlowType, 3>::
56  names[] =
57  {
58  "fixed",
59  "surfaceNormal",
60  "local"
61  };
62 
63  const NamedEnum<fv::rotorDiskSource::inletFlowType, 3>
65 }
66 
67 
68 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
69 
71 {
72  // Set inflow type
73  switch (selectionMode())
74  {
75  case smCellSet:
76  case smCellZone:
77  case smAll:
78  {
79  // Set the profile ID for each blade section
80  profiles_.connectBlades(blade_.profileName(), blade_.profileID());
81  switch (inletFlow_)
82  {
83  case ifFixed:
84  {
85  coeffs_.lookup("inletVelocity") >> inletVelocity_;
86  break;
87  }
88  case ifSurfaceNormal:
89  {
90  scalar UIn
91  (
92  readScalar(coeffs_.lookup("inletNormalVelocity"))
93  );
94  inletVelocity_ = -coordSys_.R().e3()*UIn;
95  break;
96  }
97  case ifLocal:
98  {
99  break;
100  }
101  default:
102  {
104  << "Unknown inlet velocity type" << abort(FatalError);
105  }
106  }
107 
108 
109  break;
110  }
111  default:
112  {
114  << "Source cannot be used with '"
115  << selectionModeTypeNames_[selectionMode()]
116  << "' mode. Please use one of: " << nl
117  << selectionModeTypeNames_[smCellSet] << nl
118  << selectionModeTypeNames_[smCellZone] << nl
119  << selectionModeTypeNames_[smAll]
120  << exit(FatalError);
121  }
122  }
123 }
124 
125 
127 {
128  area_ = 0.0;
129 
130  static const scalar tol = 0.8;
131 
132  const label nInternalFaces = mesh_.nInternalFaces();
133  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
134  const vectorField& Sf = mesh_.Sf();
135  const scalarField& magSf = mesh_.magSf();
136 
137  vector n = Zero;
138 
139  // Calculate cell addressing for selected cells
140  labelList cellAddr(mesh_.nCells(), -1);
141  UIndirectList<label>(cellAddr, cells_) = identity(cells_.size());
142  labelList nbrFaceCellAddr(mesh_.nFaces() - nInternalFaces, -1);
143  forAll(pbm, patchi)
144  {
145  const polyPatch& pp = pbm[patchi];
146 
147  if (pp.coupled())
148  {
149  forAll(pp, i)
150  {
151  label facei = pp.start() + i;
152  label nbrFacei = facei - nInternalFaces;
153  label own = mesh_.faceOwner()[facei];
154  nbrFaceCellAddr[nbrFacei] = cellAddr[own];
155  }
156  }
157  }
158 
159  // Correct for parallel running
160  syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
161 
162  // Add internal field contributions
163  for (label facei = 0; facei < nInternalFaces; facei++)
164  {
165  const label own = cellAddr[mesh_.faceOwner()[facei]];
166  const label nbr = cellAddr[mesh_.faceNeighbour()[facei]];
167 
168  if ((own != -1) && (nbr == -1))
169  {
170  vector nf = Sf[facei]/magSf[facei];
171 
172  if ((nf & axis) > tol)
173  {
174  area_[own] += magSf[facei];
175  n += Sf[facei];
176  }
177  }
178  else if ((own == -1) && (nbr != -1))
179  {
180  vector nf = Sf[facei]/magSf[facei];
181 
182  if ((-nf & axis) > tol)
183  {
184  area_[nbr] += magSf[facei];
185  n -= Sf[facei];
186  }
187  }
188  }
189 
190 
191  // Add boundary contributions
192  forAll(pbm, patchi)
193  {
194  const polyPatch& pp = pbm[patchi];
195  const vectorField& Sfp = mesh_.Sf().boundaryField()[patchi];
196  const scalarField& magSfp = mesh_.magSf().boundaryField()[patchi];
197 
198  if (pp.coupled())
199  {
200  forAll(pp, j)
201  {
202  const label facei = pp.start() + j;
203  const label own = cellAddr[mesh_.faceOwner()[facei]];
204  const label nbr = nbrFaceCellAddr[facei - nInternalFaces];
205  const vector nf = Sfp[j]/magSfp[j];
206 
207  if ((own != -1) && (nbr == -1) && ((nf & axis) > tol))
208  {
209  area_[own] += magSfp[j];
210  n += Sfp[j];
211  }
212  }
213  }
214  else
215  {
216  forAll(pp, j)
217  {
218  const label facei = pp.start() + j;
219  const label own = cellAddr[mesh_.faceOwner()[facei]];
220  const vector nf = Sfp[j]/magSfp[j];
221 
222  if ((own != -1) && ((nf & axis) > tol))
223  {
224  area_[own] += magSfp[j];
225  n += Sfp[j];
226  }
227  }
228  }
229  }
230 
231  if (correct)
232  {
233  reduce(n, sumOp<vector>());
234  axis = n/mag(n);
235  }
236 
237  if (debug)
238  {
239  volScalarField area
240  (
241  IOobject
242  (
243  name_ + ":area",
244  mesh_.time().timeName(),
245  mesh_,
248  ),
249  mesh_,
250  dimensionedScalar("0", dimArea, 0)
251  );
252  UIndirectList<scalar>(area.primitiveField(), cells_) = area_;
253 
254  Info<< type() << ": " << name_ << " writing field " << area.name()
255  << endl;
256 
257  area.write();
258  }
259 }
260 
261 
263 {
264  // Construct the local rotor co-prdinate system
265  vector origin(Zero);
266  vector axis(Zero);
267  vector refDir(Zero);
268 
269  geometryModeType gm =
270  geometryModeTypeNames_.read(coeffs_.lookup("geometryMode"));
271 
272  switch (gm)
273  {
274  case gmAuto:
275  {
276  // Determine rotation origin (cell volume weighted)
277  scalar sumV = 0.0;
278  const scalarField& V = mesh_.V();
279  const vectorField& C = mesh_.C();
280  forAll(cells_, i)
281  {
282  const label celli = cells_[i];
283  sumV += V[celli];
284  origin += V[celli]*C[celli];
285  }
286  reduce(origin, sumOp<vector>());
287  reduce(sumV, sumOp<scalar>());
288  origin /= sumV;
289 
290  // Determine first radial vector
291  vector dx1(Zero);
292  scalar magR = -GREAT;
293  forAll(cells_, i)
294  {
295  const label celli = cells_[i];
296  vector test = C[celli] - origin;
297  if (mag(test) > magR)
298  {
299  dx1 = test;
300  magR = mag(test);
301  }
302  }
303  reduce(dx1, maxMagSqrOp<vector>());
304  magR = mag(dx1);
305 
306  // Determine second radial vector and cross to determine axis
307  forAll(cells_, i)
308  {
309  const label celli = cells_[i];
310  vector dx2 = C[celli] - origin;
311  if (mag(dx2) > 0.5*magR)
312  {
313  axis = dx1 ^ dx2;
314  if (mag(axis) > SMALL)
315  {
316  break;
317  }
318  }
319  }
320  reduce(axis, maxMagSqrOp<vector>());
321  axis /= mag(axis);
322 
323  // Correct the axis direction using a point above the rotor
324  {
325  vector pointAbove(coeffs_.lookup("pointAbove"));
326  vector dir = pointAbove - origin;
327  dir /= mag(dir);
328  if ((dir & axis) < 0)
329  {
330  axis *= -1.0;
331  }
332  }
333 
334  coeffs_.lookup("refDirection") >> refDir;
335 
336  cylindrical_.reset
337  (
338  new cylindrical
339  (
340  mesh_,
341  axis,
342  origin,
343  cells_
344  )
345  );
346 
347  // Set the face areas and apply correction to calculated axis
348  // e.g. if cellZone is more than a single layer in thickness
349  setFaceArea(axis, true);
350 
351  break;
352  }
353  case gmSpecified:
354  {
355  coeffs_.lookup("origin") >> origin;
356  coeffs_.lookup("axis") >> axis;
357  coeffs_.lookup("refDirection") >> refDir;
358 
359  cylindrical_.reset
360  (
361  new cylindrical
362  (
363  mesh_,
364  axis,
365  origin,
366  cells_
367  )
368  );
369 
370  setFaceArea(axis, false);
371 
372  break;
373  }
374  default:
375  {
377  << "Unknown geometryMode " << geometryModeTypeNames_[gm]
378  << ". Available geometry modes include "
379  << geometryModeTypeNames_ << exit(FatalError);
380  }
381  }
382 
383  coordSys_ = cylindricalCS("rotorCoordSys", origin, axis, refDir, false);
384 
385  const scalar sumArea = gSum(area_);
386  const scalar diameter = Foam::sqrt(4.0*sumArea/mathematical::pi);
387  Info<< " Rotor gometry:" << nl
388  << " - disk diameter = " << diameter << nl
389  << " - disk area = " << sumArea << nl
390  << " - origin = " << coordSys_.origin() << nl
391  << " - r-axis = " << coordSys_.R().e1() << nl
392  << " - psi-axis = " << coordSys_.R().e2() << nl
393  << " - z-axis = " << coordSys_.R().e3() << endl;
394 }
395 
396 
398 {
399  const vectorField& C = mesh_.C();
400 
401  forAll(cells_, i)
402  {
403  if (area_[i] > ROOTVSMALL)
404  {
405  const label celli = cells_[i];
406 
407  // Position in (planar) rotor co-ordinate system
408  x_[i] = coordSys_.localPosition(C[celli]);
409 
410  // Cache max radius
411  rMax_ = max(rMax_, x_[i].x());
412 
413  // Swept angle relative to rDir axis [radians] in range 0 -> 2*pi
414  scalar psi = x_[i].y();
415 
416  // Blade flap angle [radians]
417  scalar beta =
418  flap_.beta0 - flap_.beta1c*cos(psi) - flap_.beta2s*sin(psi);
419 
420  // Determine rotation tensor to convert from planar system into the
421  // rotor cone system
422  scalar c = cos(beta);
423  scalar s = sin(beta);
424  R_[i] = tensor(c, 0, -s, 0, 1, 0, s, 0, c);
425  invR_[i] = R_[i].T();
426  }
427  }
428 }
429 
430 
432 (
433  const volVectorField& U
434 ) const
435 {
436  switch (inletFlow_)
437  {
438  case ifFixed:
439  case ifSurfaceNormal:
440  {
441  return tmp<vectorField>
442  (
443  new vectorField(mesh_.nCells(), inletVelocity_)
444  );
445 
446  break;
447  }
448  case ifLocal:
449  {
450  return U.primitiveField();
451 
452  break;
453  }
454  default:
455  {
457  << "Unknown inlet flow specification" << abort(FatalError);
458  }
459  }
460 
461  return tmp<vectorField>(new vectorField(mesh_.nCells(), Zero));
462 }
463 
464 
465 // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
466 
468 (
469  const word& name,
470  const word& modelType,
471  const dictionary& dict,
472  const fvMesh& mesh
473 
474 )
475 :
476  cellSetOption(name, modelType, dict, mesh),
477  rhoRef_(1.0),
478  omega_(0.0),
479  nBlades_(0),
480  inletFlow_(ifLocal),
481  inletVelocity_(Zero),
482  tipEffect_(1.0),
483  flap_(),
484  x_(cells_.size(), Zero),
485  R_(cells_.size(), I),
486  invR_(cells_.size(), I),
487  area_(cells_.size(), 0.0),
488  coordSys_(false),
489  cylindrical_(),
490  rMax_(0.0),
491  trim_(trimModel::New(*this, coeffs_)),
492  blade_(coeffs_.subDict("blade")),
493  profiles_(coeffs_.subDict("profiles"))
494 {
495  read(dict);
496 }
497 
498 
499 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
500 
502 {}
503 
504 
505 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
506 
508 (
509  fvMatrix<vector>& eqn,
510  const label fieldi
511 )
512 {
513  volVectorField force
514  (
515  IOobject
516  (
517  name_ + ":rotorForce",
518  mesh_.time().timeName(),
519  mesh_
520  ),
521  mesh_,
523  (
524  "zero",
525  eqn.dimensions()/dimVolume,
526  Zero
527  )
528  );
529 
530  // Read the reference density for incompressible flow
531  coeffs_.lookup("rhoRef") >> rhoRef_;
532 
533  const vectorField Uin(inflowVelocity(eqn.psi()));
534  trim_->correct(Uin, force);
535  calculate(geometricOneField(), Uin, trim_->thetag(), force);
536 
537  // Add source to rhs of eqn
538  eqn -= force;
539 
540  if (mesh_.time().writeTime())
541  {
542  force.write();
543  }
544 }
545 
546 
548 (
549  const volScalarField& rho,
550  fvMatrix<vector>& eqn,
551  const label fieldi
552 )
553 {
554  volVectorField force
555  (
556  IOobject
557  (
558  name_ + ":rotorForce",
559  mesh_.time().timeName(),
560  mesh_
561  ),
562  mesh_,
564  (
565  "zero",
566  eqn.dimensions()/dimVolume,
567  Zero
568  )
569  );
570 
571  const vectorField Uin(inflowVelocity(eqn.psi()));
572  trim_->correct(rho, Uin, force);
573  calculate(rho, Uin, trim_->thetag(), force);
574 
575  // Add source to rhs of eqn
576  eqn -= force;
577 
578  if (mesh_.time().writeTime())
579  {
580  force.write();
581  }
582 }
583 
584 
586 {
587  if (cellSetOption::read(dict))
588  {
589  coeffs_.lookup("fields") >> fieldNames_;
590  applied_.setSize(fieldNames_.size(), false);
591 
592  // Read co-ordinate system/geometry invariant properties
593  scalar rpm(readScalar(coeffs_.lookup("rpm")));
594  omega_ = rpm/60.0*mathematical::twoPi;
595 
596  coeffs_.lookup("nBlades") >> nBlades_;
597 
598  inletFlow_ = inletFlowTypeNames_.read(coeffs_.lookup("inletFlowType"));
599 
600  coeffs_.lookup("tipEffect") >> tipEffect_;
601 
602  const dictionary& flapCoeffs(coeffs_.subDict("flapCoeffs"));
603  flapCoeffs.lookup("beta0") >> flap_.beta0;
604  flapCoeffs.lookup("beta1c") >> flap_.beta1c;
605  flapCoeffs.lookup("beta2s") >> flap_.beta2s;
606  flap_.beta0 = degToRad(flap_.beta0);
607  flap_.beta1c = degToRad(flap_.beta1c);
608  flap_.beta2s = degToRad(flap_.beta2s);
609 
610 
611  // Create co-ordinate system
612  createCoordinateSystem();
613 
614  // Read co-odinate system dependent properties
615  checkData();
616 
617  constructGeometry();
618 
619  trim_->read(coeffs_);
620 
621  if (debug)
622  {
623  writeField("thetag", trim_->thetag()(), true);
624  writeField("faceArea", area_, true);
625  }
626 
627  return true;
628  }
629  else
630  {
631  return false;
632  }
633 }
634 
635 
636 // ************************************************************************* //
Collection of constants.
static autoPtr< trimModel > New(const fv::rotorDiskSource &rotor, const dictionary &dict)
Return a reference to the selected trim model.
Definition: trimModelNew.C:31
Graphite solid properties.
Definition: C.H:48
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
static const NamedEnum< inletFlowType, 3 > inletFlowTypeNames_
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:284
virtual bool read(const dictionary &dict)
Read source dictionary.
virtual ~rotorDiskSource()
Destructor.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
void read(Istream &, label &, const dictionary &)
In-place read with dictionary lookup.
Macros for easy insertion into run-time selection tables.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
rotorDiskSource(const word &name, const word &modelType, const dictionary &dict, const fvMesh &mesh)
Construct from components.
void constructGeometry()
Construct geometry.
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
Cylindrical coordinate system.
Definition: cylindricalCS.H:48
bool read(const char *, int32_t &)
Definition: int32IO.C:85
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.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
C()
Construct null.
Definition: C.C:40
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
Type gSum(const FieldField< Field, Type > &f)
dimensionedScalar cos(const dimensionedScalar &ds)
static const Identity< scalar > I
Definition: Identity.H:93
A class for handling words, derived from string.
Definition: word.H:59
labelList fv(nPoints)
A local coordinate rotation.
Definition: cylindrical.H:65
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
static const zero Zero
Definition: zero.H:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const scalar twoPi(2 *pi)
Foam::polyBoundaryMesh.
static const NamedEnum< geometryModeType, 2 > geometryModeTypeNames_
dimensionedScalar sin(const dimensionedScalar &ds)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
virtual void addSup(fvMatrix< vector > &eqn, const label fieldi)
Source term to momentum equation.
static const char nl
Definition: Ostream.H:262
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(psi *p+alphal *rhol0+((alphav *psiv+alphal *psil) - psi) *pSat, rhoMin);# 1 "/home/ubuntu/OpenFOAM-5.0/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:72
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
fileName::Type type(const fileName &, const bool followLink=true)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:485
label patchi
void checkData()
Check data.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A List with indirect addressing.
Definition: fvMatrix.H:106
const dimensionedScalar c
Speed of light in a vacuum.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
A special matrix type and solver, designed for finite volume solutions of scalar equations.
messageStream Info
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
Cell-set options abtract base class. Provides a base set of controls, e.g.:
Definition: cellSetOption.H:69
const volScalarField & psi
virtual Ostream & write(const token &)=0
Write next token to stream.
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
void setFaceArea(vector &axis, const bool correct)
Set the face areas per cell, and optionally correct the rotor axis.
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual bool write(const bool valid=true) const
Write using setting from DB.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
void createCoordinateSystem()
Create the co-ordinate system.
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
tmp< vectorField > inflowVelocity(const volVectorField &U) const
Return the inlet flow field.
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576
const dimensionSet & dimensions() const
Definition: fvMatrix.H:289