FreeFOAM The Cross-Platform CFD Toolkit
Pe.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 Application
25  Pe
26 
27 Description
28  Calculates and writes the Pe number as a surfaceScalarField obtained from
29  field phi.
30 
31  The -noWrite option just outputs the max/min values without writing
32  the field.
33 
34 Usage
35 
36  - Pe [OPTIONS]
37 
38  @param -noWrite \n
39  Suppress output to files.
40 
41  @param -dict <dictionary name>\n
42  Use named dictionary instead of system/controlDict.
43 
44  @param -noZero \n
45  Ignore timestep 0.
46 
47  @param -constant \n
48  Include the constant directory.
49 
50  @param -time <time>\n
51  Apply only to specific time.
52 
53  @param -latestTime \n
54  Only apply to latest time step.
55 
56  @param -case <dir>\n
57  Case directory.
58 
59  @param -parallel \n
60  Run in parallel.
61 
62  @param -help \n
63  Display help message.
64 
65  @param -doc \n
66  Display Doxygen API documentation page for this application.
67 
68  @param -srcDoc \n
69  Display Doxygen source documentation page for this application.
70 
71 \*---------------------------------------------------------------------------*/
72 
73 #include <postCalc/calc.H>
74 #include <finiteVolume/fvc.H>
75 
82 
83 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
84 
85 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
86 {
87  bool writeResults = !args.optionFound("noWrite");
88 
89  IOobject phiHeader
90  (
91  "phi",
92  runTime.timeName(),
93  mesh,
94  IOobject::MUST_READ
95  );
96 
97  if (phiHeader.headerOk())
98  {
99  autoPtr<surfaceScalarField> PePtr;
100 
101  Info<< " Reading phi" << endl;
102  surfaceScalarField phi(phiHeader, mesh);
103 
105  (
106  IOobject
107  (
108  "U",
109  runTime.timeName(),
110  mesh,
111  IOobject::MUST_READ
112  ),
113  mesh
114  );
115 
116  IOobject RASPropertiesHeader
117  (
118  "RASProperties",
119  runTime.constant(),
120  mesh,
121  IOobject::MUST_READ,
122  IOobject::NO_WRITE
123  );
124 
125  IOobject LESPropertiesHeader
126  (
127  "LESProperties",
128  runTime.constant(),
129  mesh,
130  IOobject::MUST_READ,
131  IOobject::NO_WRITE
132  );
133 
134  Info<< " Calculating Pe" << endl;
135 
136  if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
137  {
138  if (RASPropertiesHeader.headerOk())
139  {
140  IOdictionary RASProperties(RASPropertiesHeader);
141 
142  singlePhaseTransportModel laminarTransport(U, phi);
143 
144  autoPtr<incompressible::RASModel> RASModel
145  (
146  incompressible::RASModel::New
147  (
148  U,
149  phi,
151  )
152  );
153 
154  PePtr.set
155  (
157  (
158  IOobject
159  (
160  "Pe",
161  runTime.timeName(),
162  mesh,
163  IOobject::NO_READ
164  ),
165  mag(phi)
166  /(
167  mesh.magSf()
168  * mesh.surfaceInterpolation::deltaCoeffs()
169  * fvc::interpolate(RASModel->nuEff())
170  )
171  )
172  );
173  }
174  else if (LESPropertiesHeader.headerOk())
175  {
176  IOdictionary LESProperties(LESPropertiesHeader);
177 
178  singlePhaseTransportModel laminarTransport(U, phi);
179 
180  autoPtr<incompressible::LESModel> sgsModel
181  (
182  incompressible::LESModel::New(U, phi, laminarTransport)
183  );
184 
185  PePtr.set
186  (
188  (
189  IOobject
190  (
191  "Pe",
192  runTime.timeName(),
193  mesh,
194  IOobject::NO_READ
195  ),
196  mag(phi)
197  /(
198  mesh.magSf()
199  * mesh.surfaceInterpolation::deltaCoeffs()
200  * fvc::interpolate(sgsModel->nuEff())
201  )
202  )
203  );
204  }
205  else
206  {
207  IOdictionary transportProperties
208  (
209  IOobject
210  (
211  "transportProperties",
212  runTime.constant(),
213  mesh,
214  IOobject::MUST_READ,
215  IOobject::NO_WRITE
216  )
217  );
218 
220 
221  PePtr.set
222  (
224  (
225  IOobject
226  (
227  "Pe",
228  runTime.timeName(),
229  mesh,
230  IOobject::NO_READ
231  ),
232  mesh.surfaceInterpolation::deltaCoeffs()
233  * (mag(phi)/mesh.magSf())*(runTime.deltaT()/nu)
234  )
235  );
236  }
237  }
238  else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
239  {
240  if (RASPropertiesHeader.headerOk())
241  {
242  IOdictionary RASProperties(RASPropertiesHeader);
243 
244  autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
245 
247  (
248  IOobject
249  (
250  "rho",
251  runTime.timeName(),
252  mesh
253  ),
254  thermo->rho()
255  );
256 
257  autoPtr<compressible::RASModel> RASModel
258  (
259  compressible::RASModel::New
260  (
261  rho,
262  U,
263  phi,
264  thermo()
265  )
266  );
267 
268  PePtr.set
269  (
271  (
272  IOobject
273  (
274  "Pe",
275  runTime.timeName(),
276  mesh,
277  IOobject::NO_READ
278  ),
279  mag(phi)
280  /(
281  mesh.magSf()
282  * mesh.surfaceInterpolation::deltaCoeffs()
283  * fvc::interpolate(RASModel->muEff())
284  )
285  )
286  );
287  }
288  else if (LESPropertiesHeader.headerOk())
289  {
290  IOdictionary LESProperties(LESPropertiesHeader);
291 
292  autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
293 
295  (
296  IOobject
297  (
298  "rho",
299  runTime.timeName(),
300  mesh
301  ),
302  thermo->rho()
303  );
304 
305  autoPtr<compressible::LESModel> sgsModel
306  (
307  compressible::LESModel::New(rho, U, phi, thermo())
308  );
309 
310  PePtr.set
311  (
313  (
314  IOobject
315  (
316  "Pe",
317  runTime.timeName(),
318  mesh,
319  IOobject::NO_READ
320  ),
321  mag(phi)
322  /(
323  mesh.magSf()
324  * mesh.surfaceInterpolation::deltaCoeffs()
325  * fvc::interpolate(sgsModel->muEff())
326  )
327  )
328  );
329  }
330  else
331  {
332  IOdictionary transportProperties
333  (
334  IOobject
335  (
336  "transportProperties",
337  runTime.constant(),
338  mesh,
339  IOobject::MUST_READ,
340  IOobject::NO_WRITE
341  )
342  );
343 
345 
346  PePtr.set
347  (
349  (
350  IOobject
351  (
352  "Pe",
353  runTime.timeName(),
354  mesh,
355  IOobject::NO_READ
356  ),
357  mesh.surfaceInterpolation::deltaCoeffs()
358  * (mag(phi)/(mesh.magSf()))*(runTime.deltaT()/mu)
359  )
360  );
361  }
362  }
363  else
364  {
365  FatalErrorIn(args.executable())
366  << "Incorrect dimensions of phi: " << phi.dimensions()
367  << abort(FatalError);
368  }
369 
370 
371  // can also check how many cells exceed a particular Pe limit
372  /*
373  {
374  label count = 0;
375  label PeLimit = 200;
376  forAll(PePtr(), i)
377  {
378  if (PePtr()[i] > PeLimit)
379  {
380  count++;
381  }
382 
383  }
384 
385  Info<< "Fraction > " << PeLimit << " = "
386  << scalar(count)/Pe.size() << endl;
387  }
388  */
389 
390  Info << "Pe max : " << max(PePtr()).value() << endl;
391 
392  if (writeResults)
393  {
394  PePtr().write();
395  }
396  }
397  else
398  {
399  Info<< " No phi" << endl;
400  }
401 
402  Info<< "\nEnd\n" << endl;
403 }
404 
405 
406 // ************************ vim: set sw=4 sts=4 et: ************************ //