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