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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
30 #include "volFields.H"
31 #include "fvMatrix.H"
32 #include "fvm.H"
33 #include "mathematicalConstants.H"
34 
35 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
36 
37 template<class RhoFieldType>
38 void Foam::fv::radialActuationDiskSource::
39 addRadialActuationDiskAxialInertialResistance
40 (
41  vectorField& Usource,
42  const labelList& cells,
43  const scalarField& Vcells,
44  const RhoFieldType& rho,
45  const vectorField& U
46 )
47 {
48  scalarField Tr(cells.size());
49 
50  tensor E(Zero);
51  E.xx() = diskDir_.x();
52  E.yy() = diskDir_.y();
53  E.zz() = diskDir_.z();
54 
55  const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
56  const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);
57 
58  const vector avgCentre = gSum(zoneCellVolumes*zoneCellCentres)/V();
59  const scalar maxR = gMax(mag(zoneCellCentres - avgCentre));
60 
61  const scalar intCoeffs =
62  radialCoeffs_[0]
63  + radialCoeffs_[1]*sqr(maxR)/2.0
64  + radialCoeffs_[2]*pow4(maxR)/3.0;
65 
66  if (mag(intCoeffs) < VSMALL)
67  {
69  << "Radial distribution coefficients lead to zero polynomial." << nl
70  << "radialCoeffs = " << radialCoeffs_
71  << exit(FatalIOError);
72  }
73 
74  // Compute upstream U and rho, spatial-averaged over monitor-region
75  vector Uref(Zero);
76  scalar rhoRef = 0.0;
77  label szMonitorCells = monitorCells_.size();
78 
79  for (const auto& celli : monitorCells_)
80  {
81  Uref += U[celli];
82  rhoRef = rhoRef + rho[celli];
83  }
84  reduce(Uref, sumOp<vector>());
85  reduce(rhoRef, sumOp<scalar>());
86  reduce(szMonitorCells, sumOp<label>());
87 
88  if (szMonitorCells == 0)
89  {
91  << "No cell is available for incoming velocity monitoring."
92  << exit(FatalError);
93  }
94 
95  Uref /= szMonitorCells;
96  rhoRef /= szMonitorCells;
97 
98  const scalar Ct = sink_*UvsCtPtr_->value(mag(Uref));
99  const scalar Cp = sink_*UvsCpPtr_->value(mag(Uref));
100 
101  if (Cp <= VSMALL || Ct <= VSMALL)
102  {
104  << "Cp and Ct must be greater than zero." << nl
105  << "Cp = " << Cp << ", Ct = " << Ct
106  << exit(FatalError);
107  }
108 
109  const scalar a = 1.0 - Cp/Ct;
110  const scalar T = 2.0*rhoRef*diskArea_*mag(Uref)*a*(1.0 - a);
111 
112  forAll(cells, i)
113  {
114  const scalar r2 = magSqr(mesh().cellCentres()[cells[i]] - avgCentre);
115 
116  Tr[i] =
117  T
118  *(radialCoeffs_[0] + radialCoeffs_[1]*r2 + radialCoeffs_[2]*sqr(r2))
119  /intCoeffs;
120 
121  Usource[cells[i]] += ((Vcells[cells[i]]/V_)*Tr[i]*E) & Uref;
122  }
123 
124  if
125  (
128  )
129  {
130  Ostream& os = file();
131  writeCurrentTime(os);
132 
133  os << Uref << tab << Cp << tab << Ct << tab << a << tab << T << tab
134  << endl;
135  }
136 
137  if (debug)
138  {
139  Info<< "Source name: " << name() << nl
140  << "Average centre: " << avgCentre << nl
141  << "Maximum radius: " << maxR << endl;
142  }
143 }
144 
145 
146 // ************************************************************************* //
Foam::fv::actuationDiskSource::diskDir_
vector diskDir_
Surface-normal vector of the actuator disk pointing downstream.
Definition: actuationDiskSource.H:416
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
volFields.H
Foam::TimeState::timeOutputValue
scalar timeOutputValue() const
Return current time value.
Definition: TimeStateI.H:30
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::Vector::x
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
mathematicalConstants.H
Foam::functionObjects::writeFile::writeCurrentTime
virtual void writeCurrentTime(Ostream &os) const
Write the current time to stream.
Definition: writeFile.C:309
Foam::functionObjects::writeFile::file
virtual OFstream & file()
Return access to the file (if only 1)
Definition: writeFile.C:237
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::fv::option::name
const word & name() const
Return const access to the source name.
Definition: fvOptionI.H:30
Foam::fv::cellSetOption::V_
scalar V_
Sum of cell volumes.
Definition: cellSetOption.H:207
Foam::FatalIOError
IOerror FatalIOError
Foam::fv::cellSetOption::V
scalar V() const
Return const access to the total cell volume.
Definition: cellSetOptionI.H:69
fvMatrix.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:594
Foam::fv::actuationDiskSource::diskArea_
scalar diskArea_
Actuator disk planar surface area [m2].
Definition: actuationDiskSource.H:413
rho
rho
Definition: readInitialConditions.H:88
Foam::Vector::z
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
Foam::fv::option::mesh_
const fvMesh & mesh_
Reference to the mesh database.
Definition: fvOption.H:139
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::fv::cellSetOption::cells
const labelList & cells() const
Return const access to the cell set.
Definition: cellSetOptionI.H:75
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::pow4
dimensionedScalar pow4(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:100
Foam::fv::actuationDiskSource::writeFileStart_
scalar writeFileStart_
Start time for file output.
Definition: actuationDiskSource.H:407
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:58
fvm.H
Foam::FatalError
error FatalError
Foam::fv::actuationDiskSource::UvsCtPtr_
autoPtr< Function1< scalar > > UvsCtPtr_
Velocity vs thrust coefficients.
Definition: actuationDiskSource.H:422
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::fv::actuationDiskSource::UvsCpPtr_
autoPtr< Function1< scalar > > UvsCpPtr_
Velocity vs power coefficients.
Definition: actuationDiskSource.H:419
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::fv::actuationDiskSource::writeFileEnd_
scalar writeFileEnd_
End time for file output.
Definition: actuationDiskSource.H:410
U
U
Definition: pEqn.H:72
Foam::Vector::y
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:381
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::tab
constexpr char tab
Definition: Ostream.H:384
Foam::nl
constexpr char nl
Definition: Ostream.H:385
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Cp
const volScalarField & Cp
Definition: EEqn.H:7
Foam::fv::option::mesh
const fvMesh & mesh() const
Return const access to the mesh database.
Definition: fvOptionI.H:36
radialActuationDiskSource.H
Foam::fvMesh::time
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:275
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::fv::actuationDiskSource::sink_
label sink_
Flag for body forces to act as a source (false) or a sink (true)
Definition: actuationDiskSource.H:404
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
Foam::fv::actuationDiskSource::monitorCells_
labelList monitorCells_
Set of cells whereat the incoming velocity is monitored.
Definition: actuationDiskSource.H:425
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592