FreeFOAM The Cross-Platform CFD Toolkit
porousZone.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 "porousZone.H"
27 #include <finiteVolume/fvMesh.H>
30 
31 // * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
32 
33 // adjust negative resistance values to be multiplier of max value
34 void Foam::porousZone::adjustNegativeResistance(dimensionedVector& resist)
35 {
36  scalar maxCmpt = max(0, cmptMax(resist.value()));
37 
38  if (maxCmpt < 0)
39  {
41  (
42  "Foam::porousZone::porousZone::adjustNegativeResistance"
43  "(dimensionedVector&)"
44  ) << "negative resistances! " << resist
45  << exit(FatalError);
46  }
47  else
48  {
49  vector& val = resist.value();
50  for (label cmpt=0; cmpt < vector::nComponents; ++cmpt)
51  {
52  if (val[cmpt] < 0)
53  {
54  val[cmpt] *= -maxCmpt;
55  }
56  }
57  }
58 }
59 
60 
61 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
62 
63 Foam::porousZone::porousZone
64 (
65  const word& name,
66  const fvMesh& mesh,
67  const dictionary& dict
68 )
69 :
70  name_(name),
71  mesh_(mesh),
72  dict_(dict),
73  cellZoneID_(mesh_.cellZones().findZoneID(name)),
74  coordSys_(dict, mesh),
75  porosity_(1),
76  C0_(0),
77  C1_(0),
78  D_("D", dimensionSet(0, -2, 0, 0, 0), tensor::zero),
79  F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero)
80 {
81  Info<< "Creating porous zone: " << name_ << endl;
82 
83  bool foundZone = (cellZoneID_ != -1);
84  reduce(foundZone, orOp<bool>());
85 
86  if (!foundZone && Pstream::master())
87  {
89  (
90  "Foam::porousZone::porousZone"
91  "(const fvMesh&, const word&, const dictionary&)"
92  ) << "cannot find porous cellZone " << name_
93  << exit(FatalError);
94  }
95 
96 
97  // porosity
98  if (dict_.readIfPresent("porosity", porosity_))
99  {
100  if (porosity_ <= 0.0 || porosity_ > 1.0)
101  {
103  (
104  "Foam::porousZone::porousZone"
105  "(const fvMesh&, const word&, const dictionary&)",
106  dict_
107  )
108  << "out-of-range porosity value " << porosity_
109  << exit(FatalIOError);
110  }
111  }
112 
113  // powerLaw coefficients
114  if (const dictionary* dictPtr = dict_.subDictPtr("powerLaw"))
115  {
116  dictPtr->readIfPresent("C0", C0_);
117  dictPtr->readIfPresent("C1", C1_);
118  }
119 
120  // Darcy-Forchheimer coefficients
121  if (const dictionary* dictPtr = dict_.subDictPtr("Darcy"))
122  {
123  // local-to-global transformation tensor
124  const tensor& E = coordSys_.R();
125 
127  if (dictPtr->readIfPresent("d", d))
128  {
129  if (D_.dimensions() != d.dimensions())
130  {
132  (
133  "Foam::porousZone::porousZone"
134  "(const fvMesh&, const word&, const dictionary&)",
135  dict_
136  ) << "incorrect dimensions for d: " << d.dimensions()
137  << " should be " << D_.dimensions()
138  << exit(FatalIOError);
139  }
140 
141  adjustNegativeResistance(d);
142 
143  D_.value().xx() = d.value().x();
144  D_.value().yy() = d.value().y();
145  D_.value().zz() = d.value().z();
146  D_.value() = (E & D_ & E.T()).value();
147  }
148 
150  if (dictPtr->readIfPresent("f", f))
151  {
152  if (F_.dimensions() != f.dimensions())
153  {
155  (
156  "Foam::porousZone::porousZone"
157  "(const fvMesh&, const word&, const dictionary&)",
158  dict_
159  ) << "incorrect dimensions for f: " << f.dimensions()
160  << " should be " << F_.dimensions()
161  << exit(FatalIOError);
162  }
163 
164  adjustNegativeResistance(f);
165 
166  // leading 0.5 is from 1/2 * rho
167  F_.value().xx() = 0.5*f.value().x();
168  F_.value().yy() = 0.5*f.value().y();
169  F_.value().zz() = 0.5*f.value().z();
170  F_.value() = (E & F_ & E.T()).value();
171  }
172  }
173 
174  // provide some feedback for the user
175  // writeDict(Info, false);
176 
177  // it is an error not to define anything
178  if
179  (
180  C0_ <= VSMALL
181  && magSqr(D_.value()) <= VSMALL
182  && magSqr(F_.value()) <= VSMALL
183  )
184  {
186  (
187  "Foam::porousZone::porousZone"
188  "(const fvMesh&, const word&, const dictionary&)",
189  dict_
190  ) << "neither powerLaw (powerLaw C0/C1) "
191  "nor Darcy-Forchheimer law (Darcy d/f) specified"
192  << exit(FatalIOError);
193  }
194 }
195 
196 
197 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
198 
200 {
201  if (cellZoneID_ == -1)
202  {
203  return;
204  }
205 
206  bool compressible = false;
207  if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
208  {
209  compressible = true;
210  }
211 
212  const labelList& cells = mesh_.cellZones()[cellZoneID_];
213  const scalarField& V = mesh_.V();
214  scalarField& Udiag = UEqn.diag();
215  vectorField& Usource = UEqn.source();
216  const vectorField& U = UEqn.psi();
217 
218  if (C0_ > VSMALL)
219  {
220  if (compressible)
221  {
222  addPowerLawResistance
223  (
224  Udiag,
225  cells,
226  V,
227  mesh_.lookupObject<volScalarField>("rho"),
228  U
229  );
230  }
231  else
232  {
233  addPowerLawResistance
234  (
235  Udiag,
236  cells,
237  V,
239  U
240  );
241  }
242  }
243 
244  const tensor& D = D_.value();
245  const tensor& F = F_.value();
246 
247  if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
248  {
249  if (compressible)
250  {
251  addViscousInertialResistance
252  (
253  Udiag,
254  Usource,
255  cells,
256  V,
257  mesh_.lookupObject<volScalarField>("rho"),
258  mesh_.lookupObject<volScalarField>("mu"),
259  U
260  );
261  }
262  else
263  {
264  addViscousInertialResistance
265  (
266  Udiag,
267  Usource,
268  cells,
269  V,
271  mesh_.lookupObject<volScalarField>("nu"),
272  U
273  );
274  }
275  }
276 }
277 
278 
280 (
281  const fvVectorMatrix& UEqn,
283  bool correctAUprocBC
284 ) const
285 {
286  if (cellZoneID_ == -1)
287  {
288  return;
289  }
290 
291  bool compressible = false;
292  if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
293  {
294  compressible = true;
295  }
296 
297  const labelList& cells = mesh_.cellZones()[cellZoneID_];
298  const vectorField& U = UEqn.psi();
299 
300  if (C0_ > VSMALL)
301  {
302  if (compressible)
303  {
304  addPowerLawResistance
305  (
306  AU,
307  cells,
308  mesh_.lookupObject<volScalarField>("rho"),
309  U
310  );
311  }
312  else
313  {
314  addPowerLawResistance
315  (
316  AU,
317  cells,
319  U
320  );
321  }
322  }
323 
324  const tensor& D = D_.value();
325  const tensor& F = F_.value();
326 
327  if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
328  {
329  if (compressible)
330  {
331  addViscousInertialResistance
332  (
333  AU,
334  cells,
335  mesh_.lookupObject<volScalarField>("rho"),
336  mesh_.lookupObject<volScalarField>("mu"),
337  U
338  );
339  }
340  else
341  {
342  addViscousInertialResistance
343  (
344  AU,
345  cells,
347  mesh_.lookupObject<volScalarField>("nu"),
348  U
349  );
350  }
351  }
352 
353  if (correctAUprocBC)
354  {
355  // Correct the boundary conditions of the tensorial diagonal to ensure
356  // processor boundaries are correctly handled when AU^-1 is interpolated
357  // for the pressure equation.
359  }
360 }
361 
362 
363 void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
364 {
365  if (subDict)
366  {
367  os << indent << token::BEGIN_BLOCK << incrIndent << nl;
368  os.writeKeyword("name") << zoneName() << token::END_STATEMENT << nl;
369  }
370  else
371  {
372  os << indent << zoneName() << nl
374  }
375 
376  if (dict_.found("note"))
377  {
378  os.writeKeyword("note") << string(dict_.lookup("note"))
379  << token::END_STATEMENT << nl;
380  }
381 
382  coordSys_.writeDict(os, true);
383 
384  if (dict_.found("porosity"))
385  {
386  os.writeKeyword("porosity") << porosity() << token::END_STATEMENT << nl;
387  }
388 
389  // powerLaw coefficients
390  if (const dictionary* dictPtr = dict_.subDictPtr("powerLaw"))
391  {
392  os << indent << "powerLaw";
393  dictPtr->write(os);
394  }
395 
396  // Darcy-Forchheimer coefficients
397  if (const dictionary* dictPtr = dict_.subDictPtr("Darcy"))
398  {
399  os << indent << "Darcy";
400  dictPtr->write(os);
401  }
402 
403  os << decrIndent << indent << token::END_BLOCK << endl;
404 }
405 
406 
407 // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
408 
410 {
411  pZone.writeDict(os);
412  return os;
413 }
414 
415 // ************************ vim: set sw=4 sts=4 et: ************************ //