DAFvSourceActuatorPoint.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAFvSourceActuatorPoint.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAFvSourceActuatorPoint, 0);
16 addToRunTimeSelectionTable(DAFvSource, DAFvSourceActuatorPoint, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAFvSourceActuatorPoint::DAFvSourceActuatorPoint(
20  const word modelType,
21  const fvMesh& mesh,
22  const DAOption& daOption,
23  const DAModel& daModel,
24  const DAIndex& daIndex)
25  : DAFvSource(modelType, mesh, daOption, daModel, daIndex)
26 {
27  printIntervalUnsteady_ = daOption.getOption<label>("printIntervalUnsteady");
28  this->initFvSourcePars();
29 }
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 void DAFvSourceActuatorPoint::calcFvSource(volVectorField& fvSource)
33 {
34  /*
35  Description:
36  Compute the actuator point source term. We support two types of smooth functions:
37 
38  1. Hyperbolic:
39  thrust = xT*yT*zT where xT is
40  xT = tanh( eps*(meshCX-actCX+0.5*actSX) ) - tanh( eps*(meshCX-actCX-0.5*actSX) )
41  if the mesh cell center is within the actuator box, xT~2 other wise xT~0
42  eps is a parameter to smooth the distribution
43 
44  2. Gaussian:
45  We use a 2D Gaussian function for thrust source distribution:
46  thrust = 1/(2*pi*eps^2)*e^{-d^2/(2*eps^2)} where d is the
47  distance between the cell center and actuator source center
48 
49  NOTE: this option has some issues for U field... it has many spikes
50 
51  Example:
52  An example of the fvSource in pyOptions in pyDAFoam can be
53  defOptions =
54  {
55  "fvSource"
56  {
57  "line1"
58  {
59  "type": "actuatorPoint",
60  "smoothFunction": "hyperbolic",
61  "center": [0.0, 0.0, 0.0], # center and size define a rectangular
62  "size": [0.5, 0.3, 0.5],
63  "amplitude": [0.0, 1.0, 0.0],
64  "phase": 0.0,
65  "thrustDirIdx": 0,
66  "periodicity": 1.0,
67  "eps": 1.0,
68  "scale": 1.0 # scale the source such the integral equals desired thrust
69  },
70  "line2"
71  {
72  "type": "actuatorDisk",
73  "smoothFunction": "gaussian",
74  "center": [0.0, 0.0, 0.0],
75  "amplitude": [0.0, 1.0, 0.0],
76  "phase": 1.0,
77  "thrustDirIdx": 0,
78  "periodicity": 1.0,
79  "eps": 1.0,
80  "scale": 1.0 # scale the source such the integral equals desired thrust
81  }
82  }
83  }
84  */
85 
86  forAll(fvSource, idxI)
87  {
88  fvSource[idxI] = vector::zero;
89  }
90 
91  dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
92 
93  DAGlobalVar& globalVar =
94  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
95  HashTable<List<scalar>>& actuatorPointPars = globalVar.actuatorPointPars;
96 
97  // loop over all the cell indices for all actuator points
98  forAll(fvSourceSubDict.toc(), idxI)
99  {
100 
101  // name of this point model
102  word pointName = fvSourceSubDict.toc()[idxI];
103 
104  // sub dictionary with all parameters for this disk
105  dictionary pointSubDict = fvSourceSubDict.subDict(pointName);
106 
107  word smoothFunction = pointSubDict.get<word>("smoothFunction");
108 
109  if (smoothFunction == "hyperbolic")
110  {
111  // now read in all parameters for this actuator point
112  vector center = {
113  actuatorPointPars[pointName][0],
114  actuatorPointPars[pointName][1],
115  actuatorPointPars[pointName][2]};
116  vector size = {
117  actuatorPointPars[pointName][3],
118  actuatorPointPars[pointName][4],
119  actuatorPointPars[pointName][5]};
120  vector amp = {
121  actuatorPointPars[pointName][6],
122  actuatorPointPars[pointName][7],
123  actuatorPointPars[pointName][8]};
124  scalar period = actuatorPointPars[pointName][10];
125  scalar eps = actuatorPointPars[pointName][12];
126  scalar scale = actuatorPointPars[pointName][9];
127  label thrustDirIdx = pointSubDict.get<label>("thrustDirIdx");
128  scalar phase = actuatorPointPars[pointName][11];
129 
130  scalar t = mesh_.time().timeOutputValue();
131  center += amp * sin(constant::mathematical::twoPi * t / period + phase);
132 
133  scalar xTerm, yTerm, zTerm, s;
134  scalar thrustTotal = 0.0;
135  forAll(mesh_.C(), cellI)
136  {
137  const vector& meshC = mesh_.C()[cellI];
138  xTerm = (tanh(eps * (meshC[0] + 0.5 * size[0] - center[0])) - tanh(eps * (meshC[0] - 0.5 * size[0] - center[0])));
139  yTerm = (tanh(eps * (meshC[1] + 0.5 * size[1] - center[1])) - tanh(eps * (meshC[1] - 0.5 * size[1] - center[1])));
140  zTerm = (tanh(eps * (meshC[2] + 0.5 * size[2] - center[2])) - tanh(eps * (meshC[2] - 0.5 * size[2] - center[2])));
141 
142  s = xTerm * yTerm * zTerm;
143  // here we need to use += for multiple actuator points
144  fvSource[cellI][thrustDirIdx] += s * scale;
145 
146  thrustTotal += s * scale * mesh_.V()[cellI];
147  }
148  reduce(thrustTotal, sumOp<scalar>());
149 
150 #ifdef CODI_NO_AD
151  if (mesh_.time().timeIndex() % printIntervalUnsteady_ == 0 || mesh_.time().timeIndex() == 0)
152  {
153  Info << "Actuator point source: " << pointName << endl;
154  Info << "Total thrust source: " << thrustTotal << endl;
155  }
156 #endif
157  }
158  else if (smoothFunction == "gaussian")
159  {
160  // now read in all parameters for this actuator point
161  vector center = {
162  actuatorPointPars[pointName][0],
163  actuatorPointPars[pointName][1],
164  actuatorPointPars[pointName][2]};
165  // NOTE: the size variable is not used in the gaussian method, but we set it
166  // anyway. So the gaussian case will have empty DV keys for indices from 3 to 5
167  vector size = {
168  actuatorPointPars[pointName][3],
169  actuatorPointPars[pointName][4],
170  actuatorPointPars[pointName][5]};
171  vector amp = {
172  actuatorPointPars[pointName][6],
173  actuatorPointPars[pointName][7],
174  actuatorPointPars[pointName][8]};
175  scalar period = actuatorPointPars[pointName][10];
176  scalar eps = actuatorPointPars[pointName][12];
177  scalar scale = actuatorPointPars[pointName][9];
178  label thrustDirIdx = pointSubDict.get<label>("thrustDirIdx");
179  scalar phase = actuatorPointPars[pointName][11];
180 
181  scalar t = mesh_.time().timeOutputValue();
182  center += amp * sin(constant::mathematical::twoPi * t / period + phase);
183 
184  scalar thrustTotal = 0.0;
185  scalar coeff = 1.0 / constant::mathematical::twoPi / eps / eps;
186  forAll(mesh_.C(), cellI)
187  {
188  const vector& meshC = mesh_.C()[cellI];
189  scalar d = mag(meshC - center);
190  scalar s = coeff * exp(-d * d / 2.0 / eps / eps);
191  // here we need to use += for multiple actuator points
192  fvSource[cellI][thrustDirIdx] += s * scale;
193  thrustTotal += s * scale * mesh_.V()[cellI];
194  }
195  reduce(thrustTotal, sumOp<scalar>());
196 
197 #ifdef CODI_NO_AD
198  if (mesh_.time().timeIndex() % printIntervalUnsteady_ == 0 || mesh_.time().timeIndex() == 1)
199  {
200  Info << "Actuator point source: " << pointName << endl;
201  Info << "Total thrust source: " << thrustTotal << endl;
202  }
203 #endif
204  }
205  else
206  {
207  FatalErrorIn("") << "smoothFunction should be either hyperbolic or gaussian" << abort(FatalError);
208  }
209  }
210 
211  fvSource.correctBoundaryConditions();
212 }
213 
214 void DAFvSourceActuatorPoint::initFvSourcePars()
215 {
216  /*
217  Description:
218  Initialize the values for all types of fvSource in DAGlobalVar, including
219  actuatorDiskPars, heatSourcePars, etc
220  */
221 
222  // now we need to initialize actuatorDiskPars
223  dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
224 
225  forAll(fvSourceSubDict.toc(), idxI)
226  {
227  word pointName = fvSourceSubDict.toc()[idxI];
228  // sub dictionary with all parameters for this disk
229  dictionary pointSubDict = fvSourceSubDict.subDict(pointName);
230  word type = pointSubDict.getWord("type");
231 
232  if (type == "actuatorPoint")
233  {
234 
235  // now read in all parameters for this actuator disk
236  scalarList centerList;
237  pointSubDict.readEntry<scalarList>("center", centerList);
238 
239  scalarList sizeList;
240  word smoothFunction = pointSubDict.getWord("smoothFunction");
241  // read the size list for hyperbolic
242  if (smoothFunction == "hyperbolic")
243  {
244  pointSubDict.readEntry<scalarList>("size", sizeList);
245  }
246  else if (smoothFunction == "gaussian")
247  {
248  sizeList = {0.0, 0.0, 0.0};
249  }
250 
251  scalarList ampList;
252  pointSubDict.readEntry<scalarList>("amplitude", ampList);
253 
254  // we have 13 design variables for each line
255  scalarList dvList(13);
256  dvList[0] = centerList[0];
257  dvList[1] = centerList[1];
258  dvList[2] = centerList[2];
259  dvList[3] = sizeList[0];
260  dvList[4] = sizeList[1];
261  dvList[5] = sizeList[2];
262  dvList[6] = ampList[0];
263  dvList[7] = ampList[1];
264  dvList[8] = ampList[2];
265  dvList[9] = pointSubDict.getScalar("scale");
266  dvList[10] = pointSubDict.getScalar("periodicity");
267  dvList[11] = pointSubDict.getScalar("phase");
268  dvList[12] = pointSubDict.getScalar("eps");
269 
270  // set actuatorDiskPars
271  DAGlobalVar& globalVar =
272  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
273  globalVar.actuatorPointPars.set(pointName, dvList);
274  }
275  }
276 }
277 
278 } // End namespace Foam
279 
280 // ************************************************************************* //
Foam::DAFvSource
Definition: DAFvSource.H:34
Foam::DAFvSourceActuatorPoint::printIntervalUnsteady_
label printIntervalUnsteady_
print interval
Definition: DAFvSourceActuatorPoint.H:33
Foam::DAFvSource::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAFvSource.H:53
Foam::DAFvSourceActuatorPoint::DAFvSourceActuatorPoint
DAFvSourceActuatorPoint(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAFvSourceActuatorPoint.C:19
Foam::DAOption
Definition: DAOption.H:29
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAFvSourceActuatorPoint::calcFvSource
virtual void calcFvSource(volVectorField &fvSource)
compute the FvSource term
Definition: DAFvSourceActuatorPoint.C:32
Foam::DAOption::getOption
classType getOption(const word key) const
get an option from subDict and key
Definition: DAOption.H:92
fvSource
volScalarField & fvSource
Definition: createRefsHeatTransfer.H:7
Foam::DAOption::getAllOptions
const dictionary & getAllOptions() const
return a reference of allOptions_ dictionary
Definition: DAOption.H:56
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAIndex
Definition: DAIndex.H:32
DAFvSourceActuatorPoint.H
Foam::DAModel
Definition: DAModel.H:57
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam
Definition: checkGeometry.C:32
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAFvSource::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAFvSource.H:50
Foam::DAFvSourceActuatorPoint::initFvSourcePars
virtual void initFvSourcePars()
Initialize the values for all types of fvSource in DAGlobalVar, including actuatorDiskPars,...
Definition: DAFvSourceActuatorPoint.C:214
Foam::DAGlobalVar::actuatorPointPars
HashTable< List< scalar > > actuatorPointPars
the list of parameters for all the actuator points
Definition: DAGlobalVar.H:71
Foam::DAGlobalVar
Definition: DAGlobalVar.H:26

Generated by   doxygen 1.8.17