DAFvSourceActuatorLine.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6 \*---------------------------------------------------------------------------*/
7 
8 #include "DAFvSourceActuatorLine.H"
9 
10 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
11 
12 namespace Foam
13 {
14 
15 defineTypeNameAndDebug(DAFvSourceActuatorLine, 0);
16 addToRunTimeSelectionTable(DAFvSource, DAFvSourceActuatorLine, dictionary);
17 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
18 
19 DAFvSourceActuatorLine::DAFvSourceActuatorLine(
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 DAFvSourceActuatorLine::calcFvSource(volVectorField& fvSource)
33 {
34  /*
35  Description:
36  Compute the actuator line source term.
37  Reference: Stokkermans et al. Validation and comparison of RANS propeller
38  modeling methods for tip-mounted applications
39  NOTE: rotDir = right means propeller rotates clockwise viewed from
40  the tail of the aircraft looking forward
41  */
42 
43  const scalar& pi = Foam::constant::mathematical::pi;
44 
45  scalar t = mesh_.time().timeOutputValue();
46 
47  forAll(fvSource, idxI)
48  {
49  fvSource[idxI] = vector::zero;
50  }
51 
52  dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
53 
54  DAGlobalVar& globalVar =
55  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
56  HashTable<List<scalar>>& actuatorLinePars = globalVar.actuatorLinePars;
57 
58  // loop over all the cell indices for all actuator lines
59  forAll(fvSourceSubDict.toc(), idxI)
60  {
61  // name of this line model
62  word lineName = fvSourceSubDict.toc()[idxI];
63 
64  // sub dictionary with all parameters for this disk
65  dictionary lineSubDict = fvSourceSubDict.subDict(lineName);
66  // center of the actuator line
67  vector center = {
68  actuatorLinePars[lineName][0],
69  actuatorLinePars[lineName][1],
70  actuatorLinePars[lineName][2]};
71  // thrust direction
72  vector direction = {
73  actuatorLinePars[lineName][3],
74  actuatorLinePars[lineName][4],
75  actuatorLinePars[lineName][5]};
76  direction = direction / mag(direction);
77  // initial vector for the actuator line
78  vector initial = {
79  actuatorLinePars[lineName][6],
80  actuatorLinePars[lineName][7],
81  actuatorLinePars[lineName][8]};
82  initial = initial / mag(initial);
83  if (fabs(direction & initial) > 1.0e-10)
84  {
85  FatalErrorIn(" ") << "direction and initial need to be orthogonal!" << abort(FatalError);
86  }
87  // rotation direction, can be either right or left
88  word rotDir = lineSubDict.getWord("rotDir");
89  // inner and outer radius of the lines
90  scalar innerRadius = actuatorLinePars[lineName][9];
91  scalar outerRadius = actuatorLinePars[lineName][10];
92  // rotation speed in rpm
93  scalar rpm = actuatorLinePars[lineName][15];
94  scalar radPerS = rpm / 60.0 * 2.0 * pi;
95  // smooth factor which should be of grid size
96  scalar eps = actuatorLinePars[lineName][17];
97  // scaling factor to ensure a desired integrated thrust
98  scalar scale = actuatorLinePars[lineName][11];
99  // phase (rad) of the rotation positive value means rotates ahead of time for phase rad
100  scalar phase = actuatorLinePars[lineName][16];
101  // number of blades for this line model
102  label nBlades = lineSubDict.get<label>("nBlades");
103  scalar POD = actuatorLinePars[lineName][12];
104  scalar expM = actuatorLinePars[lineName][13];
105  scalar expN = actuatorLinePars[lineName][14];
106  // Now we need to compute normalized eps in the radial direction, i.e. epsRStar this is because
107  // we need to smooth the radial distribution of the thrust, here the radial location is
108  // normalized as rStar = (r - rInner) / (rOuter - rInner), so to make epsRStar consistent with this
109  // we need to normalize eps with the demoninator of rStar, i.e. Outer - rInner
110  scalar epsRStar = eps / (outerRadius - innerRadius);
111  scalar rStarMin = epsRStar;
112  scalar rStarMax = 1.0 - epsRStar;
113  scalar fRMin = pow(rStarMin, expM) * pow(1.0 - rStarMin, expN);
114  scalar fRMax = pow(rStarMax, expM) * pow(1.0 - rStarMax, expN);
115 
116  scalar thrustTotal = 0.0;
117  scalar torqueTotal = 0.0;
118  forAll(mesh_.C(), cellI)
119  {
120  // the cell center coordinates of this cellI
121  vector cellC = mesh_.C()[cellI];
122  // cell center to disk center vector
123  vector cellC2AVec = cellC - center;
124  // tmp tensor for calculating the axial/radial components of cellC2AVec
125  tensor cellC2AVecE(tensor::zero);
126  cellC2AVecE.xx() = cellC2AVec.x();
127  cellC2AVecE.yy() = cellC2AVec.y();
128  cellC2AVecE.zz() = cellC2AVec.z();
129  // now we need to decompose cellC2AVec into axial and radial components
130  // the axial component of cellC2AVec vector
131  vector cellC2AVecA = cellC2AVecE & direction;
132  // the radial component of cellC2AVec vector
133  vector cellC2AVecR = cellC2AVec - cellC2AVecA;
134  // now we can use the cross product to compute the tangential
135  // (circ) direction of cellI
136  vector cellC2AVecC(vector::zero);
137  if (rotDir == "left")
138  {
139  // propeller rotates counter-clockwise viewed from the tail of the aircraft looking forward
140  cellC2AVecC = cellC2AVecR ^ direction; // circ
141  }
142  else if (rotDir == "right")
143  {
144  // propeller rotates clockwise viewed from the tail of the aircraft looking forward
145  cellC2AVecC = direction ^ cellC2AVecR; // circ
146  }
147  else
148  {
149  FatalErrorIn(" ") << "rotDir not valid" << abort(FatalError);
150  }
151  // the magnitude of radial component of cellC2AVecR
152  scalar cellC2AVecRLen = mag(cellC2AVecR);
153  // the magnitude of tangential component of cellC2AVecR
154  scalar cellC2AVecCLen = mag(cellC2AVecC);
155  // the magnitude of axial component of cellC2AVecA
156  scalar cellC2AVecALen = mag(cellC2AVecA);
157  // the normalized cellC2AVecC (tangential) vector
158  vector cellC2AVecCNorm = cellC2AVecC / (cellC2AVecCLen + SMALL);
159  // smooth coefficient in the axial direction
160  scalar etaAxial = exp(-sqr(cellC2AVecALen / eps));
161 
162  scalar etaTheta = 0.0;
163  for (label bb = 0; bb < nBlades; bb++)
164  {
165  scalar thetaBlade = bb * 2.0 * pi / nBlades + radPerS * t + phase;
166 #ifdef CODI_NO_AD
167  if (cellI == 0)
168  {
169  if (mesh_.time().timeIndex() % printIntervalUnsteady_ == 0
170  || mesh_.time().timeIndex() == 1)
171  {
172  scalar twoPi = 2.0 * pi;
173  Info << "blade " << bb << " theta: "
174  //<< fmod(thetaBlade.getValue(), twoPi.getValue()) * 180.0 / pi.getValue()
175  << fmod(thetaBlade, twoPi) * 180.0 / pi
176  << " deg" << endl;
177  }
178  }
179 #endif
180  // compute the rotated vector of initial by thetaBlade degree
181  // We use a simplified version of Rodrigues rotation formulation
182  vector rotatedVec = vector::zero;
183  if (rotDir == "right")
184  {
185  rotatedVec = initial * cos(thetaBlade)
186  + (direction ^ initial) * sin(thetaBlade);
187  }
188  else if (rotDir == "left")
189  {
190  rotatedVec = initial * cos(thetaBlade)
191  + (initial ^ direction) * sin(thetaBlade);
192  }
193  else
194  {
195  FatalErrorIn(" ") << "rotDir not valid" << abort(FatalError);
196  }
197  // scale the rotated vector to have the same length as cellC2AVecR
198  rotatedVec *= cellC2AVecRLen;
199  // now we can compute the distance between the cellC2AVecR and the rotatedVec
200  scalar dS_Theta = mag(cellC2AVecR - rotatedVec);
201  // smooth coefficient in the theta direction
202  etaTheta += exp(-sqr(dS_Theta / eps));
203  }
204 
205  // now we can use Hoekstra's formulation to compute radial thrust distribution
206  scalar rPrime = cellC2AVecRLen / outerRadius;
207  scalar rPrimeHub = innerRadius / outerRadius;
208  // rStar is normalized radial location
209  scalar rStar = (rPrime - rPrimeHub) / (1.0 - rPrimeHub);
210 
211  scalar fAxial = 0.0;
212  if (rStar < rStarMin)
213  {
214  scalar dR2 = (rStar - rStarMin) * (rStar - rStarMin);
215  fAxial = fRMin * exp(-dR2 / epsRStar / epsRStar);
216  }
217  else if (rStar >= rStarMin && rStar <= rStarMax)
218  {
219  fAxial = pow(rStar, expM) * pow(1.0 - rStar, expN);
220  }
221  else
222  {
223  scalar dR2 = (rStar - rStarMax) * (rStar - rStarMax);
224  fAxial = fRMax * exp(-dR2 / epsRStar / epsRStar);
225  }
226  // we use Hoekstra's method to calculate the fCirc based on fAxial
227  // here we add 0.01*eps/outerRadius to avoid diving a zero rPrime
228  // this might happen if a cell center is very close to actuator center
229  scalar fCirc = fAxial * POD / pi / (rPrime + 0.01 * eps / outerRadius);
230 
231  vector sourceVec = (fAxial * direction + fCirc * cellC2AVecCNorm);
232 
233  fvSource[cellI] += scale * etaAxial * etaTheta * sourceVec;
234 
235  thrustTotal += scale * etaAxial * etaTheta * fAxial * mesh_.V()[cellI];
236  torqueTotal += scale * etaAxial * etaTheta * fCirc * mesh_.V()[cellI];
237  }
238  reduce(thrustTotal, sumOp<scalar>());
239  reduce(torqueTotal, sumOp<scalar>());
240 
241 #ifdef CODI_NO_AD
242  if (mesh_.time().timeIndex() % printIntervalUnsteady_ == 0 || mesh_.time().timeIndex() == 1)
243  {
244  Info << "Actuator line source: " << lineName << endl;
245  Info << "Total thrust source: " << thrustTotal << endl;
246  Info << "Total torque source: " << torqueTotal << endl;
247  }
248 #endif
249  }
250 
251  fvSource.correctBoundaryConditions();
252 }
253 
254 void DAFvSourceActuatorLine::initFvSourcePars()
255 {
256  /*
257  Description:
258  Initialize the values for all types of fvSource in DAGlobalVar, including
259  actuatorDiskPars, heatSourcePars, etc
260  */
261 
262  // now we need to initialize actuatorDiskPars
263  dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
264 
265  forAll(fvSourceSubDict.toc(), idxI)
266  {
267  word diskName = fvSourceSubDict.toc()[idxI];
268  // sub dictionary with all parameters for this disk
269  dictionary lineSubDict = fvSourceSubDict.subDict(diskName);
270  word type = lineSubDict.getWord("type");
271 
272  if (type == "actuatorLine")
273  {
274 
275  // now read in all parameters for this actuator disk
276  scalarList centerList;
277  lineSubDict.readEntry<scalarList>("center", centerList);
278 
279  scalarList dirList;
280  lineSubDict.readEntry<scalarList>("direction", dirList);
281 
282  scalarList initList;
283  lineSubDict.readEntry<scalarList>("initial", initList);
284 
285  // we have 18 design variables for each line
286  scalarList dvList(18);
287  dvList[0] = centerList[0];
288  dvList[1] = centerList[1];
289  dvList[2] = centerList[2];
290  dvList[3] = dirList[0];
291  dvList[4] = dirList[1];
292  dvList[5] = dirList[2];
293  dvList[6] = initList[0];
294  dvList[7] = initList[1];
295  dvList[8] = initList[2];
296  dvList[9] = lineSubDict.getScalar("innerRadius");
297  dvList[10] = lineSubDict.getScalar("outerRadius");
298  dvList[11] = lineSubDict.getScalar("scale");
299  dvList[12] = lineSubDict.getScalar("POD");
300  dvList[13] = lineSubDict.getScalar("expM");
301  dvList[14] = lineSubDict.getScalar("expN");
302  dvList[15] = lineSubDict.getScalar("rpm");
303  dvList[16] = lineSubDict.getScalar("phase");
304  dvList[17] = lineSubDict.getScalar("eps");
305 
306  // set actuatorDiskPars
307  DAGlobalVar& globalVar =
308  const_cast<DAGlobalVar&>(mesh_.thisDb().lookupObject<DAGlobalVar>("DAGlobalVar"));
309  globalVar.actuatorLinePars.set(diskName, dvList);
310  }
311  }
312 }
313 
314 } // End namespace Foam
315 
316 // ************************************************************************* //
Foam::DAFvSource
Definition: DAFvSource.H:34
Foam::DAFvSourceActuatorLine::DAFvSourceActuatorLine
DAFvSourceActuatorLine(const word modelType, const fvMesh &mesh, const DAOption &daOption, const DAModel &daModel, const DAIndex &daIndex)
Definition: DAFvSourceActuatorLine.C:19
DAFvSourceActuatorLine.H
Foam::DAFvSource::daOption_
const DAOption & daOption_
DAOption object.
Definition: DAFvSource.H:53
Foam::DAOption
Definition: DAOption.H:29
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
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::DAFvSourceActuatorLine::calcFvSource
virtual void calcFvSource(volVectorField &fvSource)
compute the FvSource term
Definition: DAFvSourceActuatorLine.C:32
Foam::DAIndex
Definition: DAIndex.H:32
Foam::DAModel
Definition: DAModel.H:57
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam
Definition: checkGeometry.C:32
Foam::DAFvSourceActuatorLine::initFvSourcePars
virtual void initFvSourcePars()
Initialize the values for all types of fvSource in DAGlobalVar, including actuatorDiskPars,...
Definition: DAFvSourceActuatorLine.C:254
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAFvSource::mesh_
const fvMesh & mesh_
fvMesh
Definition: DAFvSource.H:50
Foam::DAFvSourceActuatorLine::printIntervalUnsteady_
label printIntervalUnsteady_
print interval
Definition: DAFvSourceActuatorLine.H:33
Foam::DAGlobalVar
Definition: DAGlobalVar.H:26
Foam::DAGlobalVar::actuatorLinePars
HashTable< List< scalar > > actuatorLinePars
the list of parameters for all the actuator lines
Definition: DAGlobalVar.H:68

Generated by   doxygen 1.8.17