PDRobstacleTypes.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2019-2021 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "PDRobstacleTypes.H"
29#include "PDRobstacleTypes.H"
30#include "Enum.H"
31#include "unitConversion.H"
33
34using namespace Foam::constant;
35
36// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37
38#define addObstacleReader(obsType, obsName) \
39 namespace Foam \
40 { \
41 namespace PDRobstacles \
42 { \
43 addNamedToMemberFunctionSelectionTable \
44 ( \
45 PDRobstacle, \
46 obsType, \
47 read, \
48 dictionary, \
49 obsName \
50 ); \
51 } \
52 }
53
54
55// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
56
57namespace Foam
58{
59
60// Read porosity, change to blockage. Clamp values [0-1] silently
61static const scalarMinMax limits01(scalarMinMax::zero_one());
62
63// Volume porosity -> blockage
64inline scalar getPorosity(const dictionary& dict)
65{
66 return 1 - limits01.clip(dict.getOrDefault<scalar>("porosity", 0));
67}
68
69// Direction porosities -> blockage
70inline vector getPorosities(const dictionary& dict)
71{
72 vector blockage(vector::one);
73
74 if (dict.readIfPresent("porosities", blockage))
75 {
76 for (scalar& val : blockage)
77 {
78 val = 1 - limits01.clip(val);
79 }
80 }
81
82 return blockage;
83}
84
85
86// Check for "porosity", or "porosities"
87// inline static bool hasPorosity(const dictionary& dict)
88// {
89// return dict.found("porosity") || dict.found("porosities");
90// }
91
92
94vectorComponentsNames
95({
96 { vector::components::X, "x" },
97 { vector::components::Y, "y" },
98 { vector::components::Z, "z" },
99});
100
101
102enum inletDirnType
103{
104 _X = -1, // -ve x
105 _Y = -2, // -ve y
106 _Z = -3, // -ve z
107 X = 1, // +ve x
108 Y = 2, // +ve y
109 Z = 3, // +ve z
110};
111
112static const Foam::Enum<inletDirnType>
113inletDirnNames
114({
115 { inletDirnType::_X, "-x" },
116 { inletDirnType::_Y, "-y" },
117 { inletDirnType::_Z, "-z" },
118 { inletDirnType::_X, "_x" },
119 { inletDirnType::_Y, "_y" },
120 { inletDirnType::_Z, "_z" },
121 { inletDirnType::X, "+x" },
122 { inletDirnType::Y, "+y" },
123 { inletDirnType::Z, "+z" },
124 { inletDirnType::X, "x" },
125 { inletDirnType::Y, "y" },
126 { inletDirnType::Z, "z" },
127});
128
129} // End namespace Foam
130
131
132// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
133
134addObstacleReader(cylinder, cyl);
135addObstacleReader(cylinder, cylinder);
136
138(
139 PDRobstacle& obs,
140 const dictionary& dict
141)
142{
143 obs.PDRobstacle::readProperties(dict);
144 obs.typeId = enumTypeId;
145
146 // Enforce complete blockage
147 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
148 // if (hasPorosity(dict)) ... warn?
149
150
151 dict.readEntry("point", obs.pt);
152 dict.readEntry("length", obs.len());
153 dict.readEntry("diameter", obs.dia());
154
155 obs.orient = vectorComponentsNames.get("direction", dict);
156
157 // The sortBias for later position sorting
158 switch (obs.orient)
159 {
160 case vector::X:
161 obs.sortBias = obs.len();
162 break;
163
164 default:
165 obs.sortBias = 0.5*obs.dia();
166 break;
167 }
168}
169
170
171// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
172
173addObstacleReader(diagbeam, diag);
174addObstacleReader(diagbeam, diagbeam);
175
177(
178 PDRobstacle& obs,
179 const dictionary& dict
180)
181{
182 obs.PDRobstacle::readProperties(dict);
183 obs.typeId = enumTypeId;
184
185 // Enforce complete blockage
186 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
187 // if (hasPorosity(dict)) ... warn?
188
189
190 dict.readEntry("point", obs.pt);
191 dict.readEntry("length", obs.len());
192 obs.dia() = Zero;
193 obs.theta() = Zero; // Fix later on
194
195 obs.orient = vectorComponentsNames.get("direction", dict);
196
197 // Angle (degrees) on input, limit to range [0, PI]
198 scalar angle = dict.get<scalar>("angle");
199
200 while (angle > 180) angle -= 180;
201 while (angle < 0) angle += 180;
202
203 labelPair dims;
204 dict.readEntry("width", dims);
205
206 // Swap axes when theta > PI/2
207 // For 89-90 degrees it becomes -ve, which is picked up in following section
208 if (angle > 89)
209 {
210 angle -= 90;
211 dims.flip(); // Swap dimensions
212 }
213
214 obs.theta() = degToRad(angle);
215
216 obs.wa = dims.first();
217 obs.wb = dims.second();
218
219 const scalar ctheta = cos(obs.theta());
220 const scalar stheta = sin(obs.theta());
221
222 // The sortBias for later position sorting
223 switch (obs.orient)
224 {
225 case vector::X:
226 obs.sortBias = obs.len();
227 break;
228
229 case vector::Y:
230 obs.sortBias = 0.5*(obs.wa * stheta + obs.wb * ctheta);
231 break;
232
233 case vector::Z:
234 obs.sortBias = 0.5*(obs.wa * ctheta + obs.wb * stheta);
235 break;
236 }
237
238 // If very nearly aligned with axis, turn it into normal block,
239 // to avoid 1/tan(theta) blowing up
240 if (angle < 1)
241 {
242 Info<< "... changed diag-beam to box" << nl;
243
244 switch (obs.orient)
245 {
246 case vector::X:
247 obs.span = vector(obs.len(), obs.wa, obs.wb);
248 break;
249
250 case vector::Y:
251 obs.span = vector(obs.wb, obs.len(), obs.wa);
252 break;
253
254 case vector::Z:
255 obs.span = vector(obs.wa, obs.wb, obs.len());
256 break;
257 }
258
259 // The pt was end centre (cylinder), now lower corner
260 vector adjustPt = -0.5*obs.span;
261 adjustPt[obs.orient] = 0;
262
263 obs.pt -= adjustPt;
264
266 obs.sortBias = 0;
267 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1.0;
268 obs.blowoff_type = 0;
269 }
270}
271
272
273// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
274
275addObstacleReader(cuboid, box);
276
278(
279 PDRobstacle& obs,
280 const dictionary& dict
281)
282{
283 obs.PDRobstacle::readProperties(dict);
284 obs.typeId = enumTypeId;
285
286 // Default - full blockage
287 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
288
289
290 dict.readEntry("point", obs.pt);
291 dict.readEntry("size", obs.span);
292
293 // Optional
294 obs.vbkge = getPorosity(dict);
295
296 // Optional
297 const vector blockages = getPorosities(dict);
298 obs.xbkge = blockages.x();
299 obs.ybkge = blockages.y();
300 obs.zbkge = blockages.z();
301}
302
303
304// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
305
306addObstacleReader(wallbeam, wallbeam);
307
309(
310 PDRobstacle& obs,
311 const dictionary& dict
312)
313{
315 obs.typeId = enumTypeId;
316
317 // Enforce complete blockage
318 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
319
320 // if (hasPorosity(dict)) ... warn?
321
322 // Additional adjustment for drag etc.
323}
324
325
326// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
327
328addObstacleReader(grating, grating);
329addObstacleReader(grating, grate);
330
332(
333 PDRobstacle& obs,
334 const dictionary& dict
335)
336{
337 obs.PDRobstacle::readProperties(dict);
338 obs.typeId = enumTypeId;
339
340 // Initialize to full blockage
341 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
342
343 dict.readEntry("point", obs.pt);
344 dict.readEntry("size", obs.span);
345
346 // TODO: better error handling
347 // if (!equal(cmptProduct(obs.span), 0))
348 // {
349 // Info<< "Type " << typeId << " has non-zero thickness.";
350 // ReportLineInfo(lineNo, inputFile);
351 // }
352
353 obs.vbkge = getPorosity(dict);
354
355 const vector blockages = getPorosities(dict);
356 obs.xbkge = blockages.x();
357 obs.ybkge = blockages.y();
358 obs.zbkge = blockages.z();
359
360 // TODO: Warning if porosity was not specified?
361
362
363 // TBD: Default slat width from PDR.params
364 obs.slat_width = dict.getOrDefault<scalar>("slats", Zero);
365
366 // Determine the orientation
367 if (equal(obs.span.x(), 0))
368 {
369 obs.orient = vector::X;
370 }
371 else if (equal(obs.span.y(), 0))
372 {
373 obs.orient = vector::Y;
374 }
375 else
376 {
377 obs.orient = vector::Z;
378 }
379}
380
381
382// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
383
384addObstacleReader(louver, louver);
385addObstacleReader(louver, louvre);
386
388(
389 PDRobstacle& obs,
390 const dictionary& dict
391)
392{
393 obs.PDRobstacle::readProperties(dict);
394 obs.typeId = enumTypeId;
395
396 // Initialize to full blockage
397 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
398
399 dict.readEntry("point", obs.pt);
400 dict.readEntry("size", obs.span);
401
402 // TODO: better error handling
403 // if (!equal(cmptProduct(obs.span), 0))
404 // {
405 // Info<< "Type " << typeId << " has non-zero thickness.";
406 // ReportLineInfo(lineNo, inputFile);
407 // }
408
409 obs.vbkge = getPorosity(dict);
410
411 const vector blockages = getPorosities(dict);
412 obs.xbkge = blockages.x();
413 obs.ybkge = blockages.y();
414 obs.zbkge = blockages.z();
415
416 // TODO: Warning if porosity was not specified?
417
418
419 // Blowoff pressure [bar]
420 const scalar blowoffPress = dict.get<scalar>("pressure");
421
422 obs.blowoff_press = barToPa(blowoffPress);
423 obs.blowoff_time = dict.getOrDefault<scalar>("time", 0);
424 obs.blowoff_type = dict.getOrDefault<label>("type", 2);
425
426 if (obs.blowoff_type == 1)
427 {
428 Info<< "Louver : blowoff-type 1 not yet implemented." << nl;
429 // ReportLineInfo(lineNo, inputFile);
430
431 if (obs.blowoff_time != 0)
432 {
433 Info<< "Louver : has blowoff time set,"
434 << " not set to blow off cell-by-cell" << nl;
435 // ReportLineInfo(lineNo, inputFile);
436 }
437 }
438 else
439 {
440 if
441 (
442 (obs.blowoff_type == 1 || obs.blowoff_type == 2)
443 && (blowoffPress > 0)
444 )
445 {
446 if (blowoffPress > maxBlowoffPressure)
447 {
448 Info<< "Blowoff pressure (" << blowoffPress
449 << ") too high for blowoff type "
450 << obs.blowoff_type << nl;
451 // ReportLineInfo(lineNo, inputFile);
452 }
453 }
454 else
455 {
456 Info<< "Problem with blowoff parameters" << nl;
457 // ReportLineInfo(lineNo, inputFile);
458 Info<< "Pressure[bar] " << blowoffPress
459 << " Blowoff type " << obs.blowoff_type
460 << ", blowoff pressure " << obs.blowoff_press << nl;
461 }
462 }
463}
464
465
466// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
467
468addObstacleReader(patch, patch);
469
471(
472 PDRobstacle& obs,
473 const dictionary& dict
474)
475{
476 obs.PDRobstacle::readProperties(dict);
477 obs.typeId = enumTypeId;
478
479 const auto nameLen = obs.identifier.length();
480
481 word patchName = word::validate(obs.identifier);
482
483 if (patchName.empty())
484 {
486 << "RECT_PATCH without a patch name"
487 << exit(FatalError);
488 }
489 else if (patchName.length() != nameLen)
490 {
492 << "RECT_PATCH stripped invalid characters from patch name: "
493 << obs.identifier
494 << exit(FatalError);
495
496 obs.identifier = std::move(patchName);
497 }
498
499 // Enforce complete blockage
500 obs.xbkge = obs.ybkge = obs.zbkge = obs.vbkge = 1;
501
502 dict.readEntry("point", obs.pt);
503 dict.readEntry("size", obs.span);
504 obs.inlet_dirn = inletDirnNames.get("direction", dict);
505}
506
507
508// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
509
510#undef addObstacleReader
511
512// ************************************************************************* //
Macros for easy insertion into member function selection tables.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
T & first() noexcept
The first element of the list, position [0].
Definition: FixedListI.H:207
static MinMax< scalar > zero_one()
A 0-1 range corresponding to the pTraits zero, one.
Definition: MinMaxI.H:45
Obstacle definitions for PDR.
Definition: PDRobstacle.H:75
scalar sortBias
Bias for position sorting.
Definition: PDRobstacle.H:119
point pt
The obstacle location.
Definition: PDRobstacle.H:123
direction orient
The x/y/z orientation (0,1,2)
Definition: PDRobstacle.H:116
scalar theta() const
Definition: PDRobstacle.H:131
scalar len() const
Definition: PDRobstacle.H:132
scalar dia() const
Definition: PDRobstacle.H:130
vector span
The obstacle dimensions (for boxes)
Definition: PDRobstacle.H:126
int typeId
The obstacle type-id.
Definition: PDRobstacle.H:113
void flip()
Flip the Pair in-place.
Definition: PairI.H:158
const T & second() const noexcept
Return second element, which is also the last element.
Definition: PairI.H:120
virtual bool read()
Re-read model coefficients if they have changed.
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
Specialization of rigidBody to construct a cuboid given the mass and lengths of the sides.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
virtual void validate()
Validate the turbulence fields after construction.
Definition: kkLOmega.C:597
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
PtrList< volScalarField > & Y
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
#define WarningInFunction
Report a warning using Foam::Warning.
Different types of constants.
Namespace for OpenFOAM.
constexpr scalar barToPa(const scalar bar) noexcept
Conversion from bar to Pa.
dimensionedScalar sin(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
bool equal(const T &s1, const T &s2)
Compare two values for equality.
Definition: doubleFloat.H:46
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
Vector< scalar > vector
Definition: vector.H:61
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
static constexpr int enumTypeId
static constexpr int enumTypeId
Unit conversion functions.