FreeFOAM The Cross-Platform CFD Toolkit
forces.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) 1991-2010 OpenCFD Ltd.
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 "forces.H"
27 #include <finiteVolume/volFields.H>
28 #include <OpenFOAM/dictionary.H>
29 #include <OpenFOAM/Time.H>
30 
34 
38 
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 namespace Foam
43 {
44  defineTypeNameAndDebug(forces, 0);
45 }
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
50 {
51  if (obr_.foundObject<compressible::RASModel>("RASProperties"))
52  {
53  const compressible::RASModel& ras
54  = obr_.lookupObject<compressible::RASModel>("RASProperties");
55 
56  return ras.devRhoReff();
57  }
58  else if (obr_.foundObject<incompressible::RASModel>("RASProperties"))
59  {
60  const incompressible::RASModel& ras
61  = obr_.lookupObject<incompressible::RASModel>("RASProperties");
62 
63  return rho()*ras.devReff();
64  }
65  else if (obr_.foundObject<compressible::LESModel>("LESProperties"))
66  {
67  const compressible::LESModel& les =
68  obr_.lookupObject<compressible::LESModel>("LESProperties");
69 
70  return les.devRhoBeff();
71  }
72  else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))
73  {
74  const incompressible::LESModel& les
75  = obr_.lookupObject<incompressible::LESModel>("LESProperties");
76 
77  return rho()*les.devBeff();
78  }
79  else if (obr_.foundObject<basicThermo>("thermophysicalProperties"))
80  {
81  const basicThermo& thermo =
82  obr_.lookupObject<basicThermo>("thermophysicalProperties");
83 
85 
86  return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
87  }
88  else if
89  (
90  obr_.foundObject<singlePhaseTransportModel>("transportProperties")
91  )
92  {
93  const singlePhaseTransportModel& laminarT =
95  ("transportProperties");
96 
98 
99  return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
100  }
101  else if (obr_.foundObject<dictionary>("transportProperties"))
102  {
104  obr_.lookupObject<dictionary>("transportProperties");
105 
106  dimensionedScalar nu(transportProperties.lookup("nu"));
107 
109 
110  return -rho()*nu*dev(twoSymm(fvc::grad(U)));
111  }
112  else
113  {
114  FatalErrorIn("forces::devRhoReff()")
115  << "No valid model for viscous stress calculation."
116  << exit(FatalError);
117 
118  return volSymmTensorField::null();
119  }
120 }
121 
122 
124 {
125  if (rhoName_ == "rhoInf")
126  {
127  const fvMesh& mesh = refCast<const fvMesh>(obr_);
128 
129  return tmp<volScalarField>
130  (
131  new volScalarField
132  (
133  IOobject
134  (
135  "rho",
136  mesh.time().timeName(),
137  mesh
138  ),
139  mesh,
140  dimensionedScalar("rho", dimDensity, rhoRef_)
141  )
142  );
143  }
144  else
145  {
146  return(obr_.lookupObject<volScalarField>(rhoName_));
147  }
148 }
149 
150 
151 Foam::scalar Foam::forces::rho(const volScalarField& p) const
152 {
153  if (p.dimensions() == dimPressure)
154  {
155  return 1.0;
156  }
157  else
158  {
159  if (rhoName_ != "rhoInf")
160  {
161  FatalErrorIn("forces::rho(const volScalarField& p)")
162  << "Dynamic pressure is expected but kinematic is provided."
163  << exit(FatalError);
164  }
165 
166  return rhoRef_;
167  }
168 }
169 
170 
171 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
172 
174 (
175  const word& name,
176  const objectRegistry& obr,
177  const dictionary& dict,
178  const bool loadFromFiles
179 )
180 :
181  name_(name),
182  obr_(obr),
183  active_(true),
184  log_(false),
185  patchSet_(),
186  pName_(word::null),
187  UName_(word::null),
188  rhoName_(word::null),
189  directForceDensity_(false),
190  fDName_(""),
191  rhoRef_(VGREAT),
192  pRef_(0),
193  CofR_(vector::zero),
194  forcesFilePtr_(NULL)
195 {
196  // Check if the available mesh is an fvMesh otherise deactivate
197  if (!isA<fvMesh>(obr_))
198  {
199  active_ = false;
200  WarningIn
201  (
202  "Foam::forces::forces"
203  "("
204  "const word&, "
205  "const objectRegistry&, "
206  "const dictionary&, "
207  "const bool"
208  ")"
209  ) << "No fvMesh available, deactivating."
210  << endl;
211  }
212 
213  read(dict);
214 }
215 
216 
217 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
218 
220 {}
221 
222 
223 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
224 
226 {
227  if (active_)
228  {
229  log_ = dict.lookupOrDefault<Switch>("log", false);
230 
231  const fvMesh& mesh = refCast<const fvMesh>(obr_);
232 
233  patchSet_ =
234  mesh.boundaryMesh().patchSet(wordList(dict.lookup("patches")));
235 
236  dict.readIfPresent("directForceDensity", directForceDensity_);
237 
238  if (directForceDensity_)
239  {
240  // Optional entry for fDName
241  fDName_ = dict.lookupOrDefault<word>("fDName", "fD");
242 
243  // Check whether fDName exists, if not deactivate forces
244  if
245  (
246  !obr_.foundObject<volVectorField>(fDName_)
247  )
248  {
249  active_ = false;
250  WarningIn("void forces::read(const dictionary& dict)")
251  << "Could not find " << fDName_ << " in database." << nl
252  << " De-activating forces."
253  << endl;
254  }
255  }
256  else
257  {
258  // Optional entries U and p
259  pName_ = dict.lookupOrDefault<word>("pName", "p");
260  UName_ = dict.lookupOrDefault<word>("UName", "U");
261  rhoName_ = dict.lookupOrDefault<word>("rhoName", "rho");
262 
263  // Check whether UName, pName and rhoName exists,
264  // if not deactivate forces
265  if
266  (
267  !obr_.foundObject<volVectorField>(UName_)
268  || !obr_.foundObject<volScalarField>(pName_)
269  || (
270  rhoName_ != "rhoInf"
271  && !obr_.foundObject<volScalarField>(rhoName_)
272  )
273  )
274  {
275  active_ = false;
276 
277  WarningIn("void forces::read(const dictionary& dict)")
278  << "Could not find " << UName_ << ", " << pName_;
279 
280  if (rhoName_ != "rhoInf")
281  {
282  Info<< " or " << rhoName_;
283  }
284 
285  Info<< " in database." << nl << " De-activating forces."
286  << endl;
287  }
288 
289  // Reference density needed for incompressible calculations
290  rhoRef_ = readScalar(dict.lookup("rhoInf"));
291 
292  // Reference pressure, 0 by default
293  pRef_ = dict.lookupOrDefault<scalar>("pRef", 0.0);
294  }
295 
296  // Centre of rotation for moment calculations
297  CofR_ = dict.lookup("CofR");
298  }
299 }
300 
301 
303 {
304  // Create the forces file if not already created
305  if (forcesFilePtr_.empty())
306  {
307  if (debug)
308  {
309  Info<< "Creating forces file." << endl;
310  }
311 
312  // File update
313  if (Pstream::master())
314  {
315  fileName forcesDir;
316  word startTimeName =
317  obr_.time().timeName(obr_.time().startTime().value());
318 
319  if (Pstream::parRun())
320  {
321  // Put in undecomposed case (Note: gives problems for
322  // distributed data running)
323  forcesDir = obr_.time().path()/".."/name_/startTimeName;
324  }
325  else
326  {
327  forcesDir = obr_.time().path()/name_/startTimeName;
328  }
329 
330  // Create directory if does not exist.
331  mkDir(forcesDir);
332 
333  // Open new file at start up
334  forcesFilePtr_.reset(new OFstream(forcesDir/(type() + ".dat")));
335 
336  // Add headers to output data
337  writeFileHeader();
338  }
339  }
340 }
341 
342 
344 {
345  if (forcesFilePtr_.valid())
346  {
347  forcesFilePtr_()
348  << "# Time" << tab
349  << "forces(pressure, viscous) moment(pressure, viscous)"
350  << endl;
351  }
352 }
353 
354 
356 {
357  // Do nothing - only valid on write
358 }
359 
360 
362 {
363  // Do nothing - only valid on write
364 }
365 
366 
368 {
369  if (active_)
370  {
371  // Create the forces file if not already created
372  makeFile();
373 
374  forcesMoments fm = calcForcesMoment();
375 
376  if (Pstream::master())
377  {
378  forcesFilePtr_() << obr_.time().value() << tab << fm << endl;
379 
380  if (log_)
381  {
382  Info<< "forces output:" << nl
383  << " forces(pressure, viscous)" << fm.first() << nl
384  << " moment(pressure, viscous)" << fm.second() << nl
385  << endl;
386  }
387  }
388  }
389 }
390 
391 
393 {
394  forcesMoments fm
395  (
398  );
399 
400  if (directForceDensity_)
401  {
402  const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
403 
404  const fvMesh& mesh = fD.mesh();
405 
407  mesh.Sf().boundaryField();
408 
409  forAllConstIter(labelHashSet, patchSet_, iter)
410  {
411  label patchi = iter.key();
412 
413  vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
414 
415  scalarField sA = mag(Sfb[patchi]);
416 
417  // Normal force = surfaceUnitNormal * (surfaceNormal & forceDensity)
418  vectorField fN =
419  Sfb[patchi]/sA
420  *(
421  Sfb[patchi] & fD.boundaryField()[patchi]
422  );
423 
424  fm.first().first() += sum(fN);
425  fm.second().first() += sum(Md ^ fN);
426 
427  // Tangential force (total force minus normal fN)
428  vectorField fT = sA*fD.boundaryField()[patchi] - fN;
429 
430  fm.first().second() += sum(fT);
431  fm.second().second() += sum(Md ^ fT);
432  }
433  }
434  else
435  {
436  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
437  const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
438 
439  const fvMesh& mesh = U.mesh();
440 
442  mesh.Sf().boundaryField();
443 
444  tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
445  const volSymmTensorField::GeometricBoundaryField& devRhoReffb
446  = tdevRhoReff().boundaryField();
447 
448  // Scale pRef by density for incompressible simulations
449  scalar pRef = pRef_/rho(p);
450 
451  forAllConstIter(labelHashSet, patchSet_, iter)
452  {
453  label patchi = iter.key();
454 
455  vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;
456 
457  vectorField pf = Sfb[patchi]*(p.boundaryField()[patchi] - pRef);
458 
459  fm.first().first() += rho(p)*sum(pf);
460  fm.second().first() += rho(p)*sum(Md ^ pf);
461 
462  vectorField vf = Sfb[patchi] & devRhoReffb[patchi];
463 
464  fm.first().second() += sum(vf);
465  fm.second().second() += sum(Md ^ vf);
466  }
467  }
468 
469  reduce(fm, sumOp());
470 
471  return fm;
472 }
473 
474 
475 // ************************ vim: set sw=4 sts=4 et: ************************ //