targetCoeffTrim.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2021 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 "targetCoeffTrim.H"
27 #include "geometricOneField.H"
29 
30 using namespace Foam::constant;
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(targetCoeffTrim, 0);
37 
38  addToRunTimeSelectionTable(trimModel, targetCoeffTrim, dictionary);
39 }
40 
41 
42 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
43 
44 template<class RhoFieldType>
46 (
47  const RhoFieldType& rho,
48  const vectorField& U,
49  const scalarField& thetag,
50  vectorField& force
51 ) const
52 {
53  rotor_.calculate(rho, U, thetag, force, false, false);
54 
55  const labelList& cells = rotor_.set().cells();
56  const vectorField& C = rotor_.mesh().C();
57  const List<point>& x = rotor_.x();
58 
59  const vector& origin = rotor_.coordSys().origin();
60  const vector& rollAxis = rotor_.coordSys().R().e1();
61  const vector& pitchAxis = rotor_.coordSys().R().e2();
62  const vector& yawAxis = rotor_.coordSys().R().e3();
63 
64  scalar coeff1 = alpha_*sqr(rotor_.omega())*mathematical::pi;
65 
66  vector cf(Zero);
67  forAll(cells, i)
68  {
69  label celli = cells[i];
70 
71  vector fc = force[celli];
72  vector mc = fc^(C[celli] - origin);
73 
74  if (useCoeffs_)
75  {
76  scalar radius = x[i].x();
77  scalar coeff2 = rho[celli]*coeff1*pow4(radius);
78 
79  // add to coefficient vector
80  cf[0] += (fc & yawAxis)/(coeff2 + rootVSmall);
81  cf[1] += (mc & pitchAxis)/(coeff2*radius + rootVSmall);
82  cf[2] += (mc & rollAxis)/(coeff2*radius + rootVSmall);
83  }
84  else
85  {
86  cf[0] += fc & yawAxis;
87  cf[1] += mc & pitchAxis;
88  cf[2] += mc & rollAxis;
89  }
90  }
91 
92  reduce(cf, sumOp<vector>());
93 
94  return cf;
95 }
96 
97 
98 template<class RhoFieldType>
100 (
101  const RhoFieldType& rho,
102  const vectorField& U,
103  vectorField& force
104 ) const
105 {
106  if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
107  {
108  word calcType = "forces";
109  if (useCoeffs_)
110  {
111  calcType = "coefficients";
112  }
113 
114  Info<< type() << ":" << nl
115  << " solving for target trim " << calcType << nl;
116 
117  const scalar rhoRef = rotor_.rhoRef();
118 
119  // iterate to find new pitch angles to achieve target force
120  scalar err = great;
121  label iter = 0;
122  tensor J(Zero);
123 
124  vector old = Zero;
125  while ((err > tol_) && (iter < nIter_))
126  {
127  // cache initial theta vector
128  vector theta0(theta_);
129 
130  // set initial values
131  old = calcCoeffs(rho, U, thetag(), force);
132 
133  // construct Jacobian by perturbing the pitch angles
134  // by +/-(dTheta_/2)
135  for (label pitchI = 0; pitchI < 3; pitchI++)
136  {
137  theta_[pitchI] -= dTheta_/2.0;
138  vector cf0 = calcCoeffs(rho, U, thetag(), force);
139 
140  theta_[pitchI] += dTheta_;
141  vector cf1 = calcCoeffs(rho, U, thetag(), force);
142 
143  vector ddTheta = (cf1 - cf0)/dTheta_;
144 
145  J[pitchI + 0] = ddTheta[0];
146  J[pitchI + 3] = ddTheta[1];
147  J[pitchI + 6] = ddTheta[2];
148 
149  theta_ = theta0;
150  }
151 
152  // calculate the change in pitch angle vector
153  vector dt = inv(J) & (target_/rhoRef - old);
154 
155  // update pitch angles
156  vector thetaNew = theta_ + relax_*dt;
157 
158  // update error
159  err = mag(thetaNew - theta_);
160 
161  // update for next iteration
162  theta_ = thetaNew;
163  iter++;
164  }
165 
166  if (iter == nIter_)
167  {
168  Info<< " solution not converged in " << iter
169  << " iterations, final residual = " << err
170  << "(" << tol_ << ")" << endl;
171  }
172  else
173  {
174  Info<< " final residual = " << err << "(" << tol_
175  << "), iterations = " << iter << endl;
176  }
177 
178  Info<< " current and target " << calcType << nl
179  << " thrust = " << old[0]*rhoRef << ", " << target_[0] << nl
180  << " pitch = " << old[1]*rhoRef << ", " << target_[1] << nl
181  << " roll = " << old[2]*rhoRef << ", " << target_[2] << nl
182  << " new pitch angles [deg]:" << nl
183  << " theta0 = " << radToDeg(theta_[0]) << nl
184  << " theta1c = " << radToDeg(theta_[1]) << nl
185  << " theta1s = " << radToDeg(theta_[2]) << nl
186  << endl;
187  }
188 }
189 
190 
191 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
192 
194 (
195  const fv::rotorDiskSource& rotor,
196  const dictionary& dict
197 )
198 :
199  trimModel(rotor, dict, typeName),
200  calcFrequency_(-1),
201  useCoeffs_(true),
202  target_(Zero),
203  theta_(Zero),
204  nIter_(50),
205  tol_(1e-8),
206  relax_(1.0),
207  dTheta_(degToRad(0.1)),
208  alpha_(1.0)
209 {
210  read(dict);
211 }
212 
213 
214 // * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * //
215 
217 {}
218 
219 
220 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
221 
223 {
224  trimModel::read(dict);
225 
226  const dictionary& targetDict(coeffs_.subDict("target"));
227  useCoeffs_ = targetDict.lookupOrDefault<bool>("useCoeffs", true);
228  word ext = "";
229  if (useCoeffs_)
230  {
231  ext = "Coeff";
232  }
233 
234  target_[0] = targetDict.lookup<scalar>("thrust" + ext);
235  target_[1] = targetDict.lookup<scalar>("pitch" + ext);
236  target_[2] = targetDict.lookup<scalar>("roll" + ext);
237 
238  const dictionary& pitchAngleDict(coeffs_.subDict("pitchAngles"));
239  theta_[0] = degToRad(pitchAngleDict.lookup<scalar>("theta0Ini"));
240  theta_[1] = degToRad(pitchAngleDict.lookup<scalar>("theta1cIni"));
241  theta_[2] = degToRad(pitchAngleDict.lookup<scalar>("theta1sIni"));
242 
243  coeffs_.lookup("calcFrequency") >> calcFrequency_;
244 
245  coeffs_.readIfPresent("nIter", nIter_);
246  coeffs_.readIfPresent("tol", tol_);
247  coeffs_.readIfPresent("relax", relax_);
248 
249  if (coeffs_.readIfPresent("dTheta", dTheta_))
250  {
251  dTheta_ = degToRad(dTheta_);
252  }
253 
254  alpha_ = coeffs_.lookup<scalar>("alpha");
255 }
256 
257 
259 {
260  const List<vector>& x = rotor_.x();
261 
262  tmp<scalarField> ttheta(new scalarField(x.size()));
263  scalarField& t = ttheta.ref();
264 
265  forAll(t, i)
266  {
267  scalar psi = x[i].y();
268  t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi);
269  }
270 
271  return ttheta;
272 }
273 
274 
276 (
277  const vectorField& U,
278  vectorField& force
279 ) const
280 {
281  correctTrim(geometricOneField(), U, force);
282 }
283 
284 
286 (
287  const volScalarField rho,
288  const vectorField& U,
289  vectorField& force
290 ) const
291 {
292  correctTrim(rho, U, force);
293 }
294 
295 
296 // ************************************************************************* //
Collection of constants.
Graphite solid properties.
Definition: C.H:48
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
void correctTrim(const RhoFieldType &rho, const vectorField &U, vectorField &force) const
Correct the model.
virtual void read(const dictionary &dict)
Read.
Definition: trimModel.C:61
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Trim model base class.
Definition: trimModel.H:50
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Macros for easy insertion into run-time selection tables.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Cell based momentum source which approximates the mean effects of rotor forces on a cylindrical regio...
bool read(const char *, int32_t &)
Definition: int32IO.C:85
virtual ~targetCoeffTrim()
Destructor.
C()
Construct null.
Definition: C.C:40
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
#define calcType(Type, GeoField)
virtual void correct(const vectorField &U, vectorField &force) const
Correct the model.
static const zero Zero
Definition: zero.H:97
vector calcCoeffs(const RhoFieldType &rho, const vectorField &U, const scalarField &alphag, vectorField &force) const
Calculate the rotor force and moment coefficients vector.
virtual tmp< scalarField > thetag() const
Return the geometric angle of attack [rad].
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross product operators.
Definition: Vector.H:57
dimensionedScalar sin(const dimensionedScalar &ds)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
targetCoeffTrim(const fv::rotorDiskSource &rotor, const dictionary &dict)
Constructor.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
dimensionedScalar pow4(const dimensionedScalar &ds)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
void read(const dictionary &dict)
Read.
const volScalarField & psi
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A class for managing temporary objects.
Definition: PtrList.H:53
Namespace for OpenFOAM.