radialActuationDiskSourceTemplates.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 
27 #include "volFields.H"
28 #include "fvMatrix.H"
29 #include "fvm.H"
30 #include "mathematicalConstants.H"
31 
32 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
33 
34 template<class RhoFieldType>
35 void Foam::fv::radialActuationDiskSource::
36 addRadialActuationDiskAxialInertialResistance
37 (
38  vectorField& Usource,
39  const labelList& cells,
40  const scalarField& Vcells,
41  const RhoFieldType& rho,
42  const vectorField& U
43 ) const
44 {
45  scalar a = 1.0 - Cp_/Ct_;
46  scalarField Tr(cells.size());
47  const vector uniDiskDir = diskDir_/mag(diskDir_);
48 
49  tensor E(Zero);
50  E.xx() = uniDiskDir.x();
51  E.yy() = uniDiskDir.y();
52  E.zz() = uniDiskDir.z();
53 
54  const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
55  const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);
56 
57  const vector avgCentre = gSum(zoneCellVolumes*zoneCellCentres)/V();
58  const scalar maxR = gMax(mag(zoneCellCentres - avgCentre));
59 
60  scalar intCoeffs =
61  radialCoeffs_[0]
62  + radialCoeffs_[1]*sqr(maxR)/2.0
63  + radialCoeffs_[2]*pow4(maxR)/3.0;
64 
65  vector upU = vector(VGREAT, VGREAT, VGREAT);
66  scalar upRho = VGREAT;
67  if (upstreamCellId_ != -1)
68  {
69  upU = U[upstreamCellId_];
70  upRho = rho[upstreamCellId_];
71  }
72  reduce(upU, minOp<vector>());
73  reduce(upRho, minOp<scalar>());
74 
75  scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1.0 - a);
76  forAll(cells, i)
77  {
78  scalar r2 = magSqr(mesh().cellCentres()[cells[i]] - avgCentre);
79 
80  Tr[i] =
81  T
82  *(radialCoeffs_[0] + radialCoeffs_[1]*r2 + radialCoeffs_[2]*sqr(r2))
83  /intCoeffs;
84 
85  Usource[cells[i]] += ((Vcells[cells[i]]/V_)*Tr[i]*E) & upU;
86  }
87 
88  if (debug)
89  {
90  Info<< "Source name: " << name() << nl
91  << "Average centre: " << avgCentre << nl
92  << "Maximum radius: " << maxR << endl;
93  }
94 }
95 
96 
97 // ************************************************************************* //
vector diskDir_
Disk area normal.
label upstreamCellId_
Upstream cell ID.
scalar V_
Sum of cell volumes.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
volVectorField vectorField(fieldObject, mesh)
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const fvMesh & mesh() const
Return const access to the mesh database.
Definition: fvOptionI.H:34
scalar V() const
Return const access to the total cell volume.
Type gSum(const FieldField< Field, Type > &f)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
const word & name() const
Return const access to the source name.
Definition: fvOptionI.H:28
static const zero Zero
Definition: zero.H:91
dimensioned< scalar > magSqr(const dimensioned< Type > &)
volScalarField scalarField(fieldObject, mesh)
static const char nl
Definition: Ostream.H:262
scalar Cp_
Power coefficient.
Type gMax(const FieldField< Field, Type > &f)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
scalar Ct_
Thrust coefficient.
dimensionedScalar pow4(const dimensionedScalar &ds)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51