comfort.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) 2019-2023 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 "comfort.H"
27 #include "wallFvPatch.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 namespace functionObjects
35 {
38 }
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::magU() const
45 {
46  tmp<volScalarField> tmagU = mag(lookupObject<volVectorField>("U"));
47 
48  // Switch to use the averaged velocity field in the domain.
49  // Consistent with EN ISO 7730 but does not make physical sense
50  if (meanVelocity_)
51  {
52  tmagU.ref() = tmagU->weightedAverage(mesh_.V());
53  }
54 
55  return tmagU;
56 }
57 
58 
59 Foam::dimensionedScalar Foam::functionObjects::comfort::Trad() const
60 {
61  dimensionedScalar Trad(Trad_);
62 
63  // The mean radiation is calculated by the mean wall temperatures
64  // which are summed and divided by the area | only walls are taken into
65  // account. This approach might be correct for a squared room but will
66  // defintely be inconsistent for complex room geometries. The norm does
67  // not provide any information about the calculation of this quantity.
68  if (!TradSet_)
69  {
70  const volScalarField::Boundary& TBf =
71  lookupObject<volScalarField>("T").boundaryField();
72 
73  scalar areaIntegral = 0;
74  scalar TareaIntegral = 0;
75 
76  forAll(TBf, patchi)
77  {
78  const fvPatchScalarField& pT = TBf[patchi];
79  const fvPatch& pTBf = TBf[patchi].patch();
80  const scalarField& pSf = pTBf.magSf();
81 
82  if (isType<wallFvPatch>(pTBf))
83  {
84  areaIntegral += gSum(pSf);
85  TareaIntegral += gSum(pSf*pT);
86  }
87  }
88 
89  Trad.value() = TareaIntegral/areaIntegral;
90  }
91 
92  // Bounds based on EN ISO 7730
93  if ((Trad.value() < 283.15) || (Trad.value() > 313.15))
94  {
96  << "The calculated mean wall radiation temperature is out of the\n"
97  << "bounds specified in EN ISO 7730:2006\n"
98  << "Valid range is 10 degC < T < 40 degC\n"
99  << "The actual value is: " << Trad.value() - 273.15 << nl << endl;
100  }
101 
102  return Trad;
103 }
104 
105 
106 Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::pSat() const
107 {
108  static const dimensionedScalar kPaToPa(dimPressure, 1000);
109  static const dimensionedScalar A(dimless, 16.6563);
110  static const dimensionedScalar B(dimTemperature, 4030.183);
111  static const dimensionedScalar C(dimTemperature, -38.15);
112 
113  tmp<volScalarField> tpSat(volScalarField::New("pSat", mesh_, pSat_));
114 
115  // Calculate the saturation pressure if no user input is given
116  if (pSat_.value() == 0)
117  {
118  const volScalarField& T = lookupObject<volScalarField>("T");
119 
120  // Equation based on ISO 7730:2006
121  tpSat = kPaToPa*exp(A - B/(T + C));
122  }
123 
124  return tpSat;
125 }
126 
127 
128 Foam::tmp<Foam::volScalarField> Foam::functionObjects::comfort::Tcloth
129 (
130  volScalarField& hc,
131  const dimensionedScalar& metabolicRateSI,
132  const dimensionedScalar& extWorkSI,
133  const volScalarField& T,
134  const dimensionedScalar& Trad
135 )
136 {
137  const dimensionedScalar factor1(dimTemperature, 308.85);
138 
139  const dimensionedScalar factor2
140  (
141  dimTemperature/metabolicRateSI.dimensions(),
142  0.028
143  );
144 
145  const dimensionedScalar factor3
146  (
147  dimensionSet(1, 0, -3, -4, 0, 0, 0),
148  3.96e-8
149  );
150 
151  // Heat transfer coefficient based on forced convection [W/m^2/K]
152  const volScalarField hcForced
153  (
154  dimensionedScalar(hc.dimensions()/sqrt(dimVelocity), 12.1)
155  *sqrt(magU())
156  );
157 
158  // Tcl [K] (surface cloth temperature)
159  tmp<volScalarField> tTcl
160  (
162  (
163  "Tcl",
164  T.mesh(),
166  )
167  );
168 
169  volScalarField& Tcl = tTcl.ref();
170 
171  // Initial guess
172  Tcl = T;
173 
174  label i = 0;
175 
176  Tcl.storePrevIter();
177 
178  // Iterative solving of equation (2)
179  do
180  {
181  Tcl = (Tcl + Tcl.prevIter())/2;
182  Tcl.storePrevIter();
183 
184  // Heat transfer coefficient based on natural convection
185  volScalarField hcNatural
186  (
187  dimensionedScalar(hc.dimensions()/pow025(dimTemperature), 2.38)
188  *pow025(mag(Tcl - T))
189  );
190 
191  // Set heat transfer coefficient based on equation (3)
192  hc =
193  pos(hcForced - hcNatural)*hcForced
194  + neg0(hcForced - hcNatural)*hcNatural;
195 
196  // Calculate surface temperature based on equation (2)
197  Tcl =
198  factor1
199  - factor2*(metabolicRateSI - extWorkSI)
200  - Icl_*factor3*fcl_*(pow4(Tcl) - pow4(Trad))
201  - Icl_*fcl_*hc*(Tcl - T);
202 
203  } while (!converged(Tcl) && i++ < maxClothIter_);
204 
205  if (i == maxClothIter_)
206  {
208  << "The surface cloth temperature did not converge within " << i
209  << " iterations\n";
210  }
211 
212  return tTcl;
213 }
214 
215 
216 bool Foam::functionObjects::comfort::converged
217 (
218  const volScalarField& phi
219 ) const
220 {
221  return
222  max(mag(phi.primitiveField() - phi.prevIter().primitiveField()))
223  < tolerance_;
224 }
225 
226 
227 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
228 
230 (
231  const word& name,
232  const Time& runTime,
233  const dictionary& dict
234 )
235 :
236  fvMeshFunctionObject(name, runTime, dict),
237  clothing_("clothing", dimless, 0),
238  metabolicRate_("metabolicRate", dimMass/pow3(dimTime), 0.8),
239  extWork_("extWork", dimMass/pow3(dimTime), 0),
240  TradSet_(false),
241  Trad_("Trad", dimTemperature, 0),
242  relHumidity_("relHumidity", dimless, 0.5),
243  pSat_("pSat", dimPressure, 0),
244  Icl_("Icl", dimensionSet(-1, 0, 3, 1, 0, 0, 0), 0),
245  fcl_("fcl", dimless, 0),
246  tolerance_(1e-4),
247  maxClothIter_(100),
248  meanVelocity_(false)
249 {
250  read(dict);
251 }
252 
253 
254 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
255 
257 {}
258 
259 
260 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
261 
263 {
264  clothing_.readIfPresent(dict);
265  metabolicRate_.readIfPresent(dict);
266  extWork_.readIfPresent(dict);
267  pSat_.readIfPresent(dict);
268  tolerance_ = dict.lookupOrDefault("tolerance", 1e-4);
269  maxClothIter_ = dict.lookupOrDefault("maxClothIter", 100);
270  meanVelocity_ = dict.lookupOrDefault<Switch>("meanVelocity", false);
271 
272  // Read relative humidity if provided and convert from % to fraction
273  if (dict.found(relHumidity_.name()))
274  {
275  relHumidity_.read(dict);
276  relHumidity_ /= 100;
277  }
278 
279  // Read radiation temperature if provided
280  if (dict.found(Trad_.name()))
281  {
282  TradSet_ = true;
283  Trad_.read(dict);
284  }
285  else
286  {
287  TradSet_ = false;
288  }
289 
290  Icl_ = dimensionedScalar(Icl_.dimensions(), 0.155)*clothing_;
291 
292  fcl_.value() =
293  Icl_.value() <= 0.078
294  ? 1.0 + 1.290*Icl_.value()
295  : 1.05 + 0.645*Icl_.value();
296 
297  return true;
298 }
299 
300 
302 {
303  return wordList{"U", "T"};
304 }
305 
306 
308 {
309  const dimensionedScalar Trad(this->Trad());
310  const volScalarField pSat(this->pSat());
311 
312  const dimensionedScalar metabolicRateSI(58.15*metabolicRate_);
313  const dimensionedScalar extWorkSI(58.15*extWork_);
314 
315  const volScalarField& T = lookupObject<volScalarField>("T");
316 
317  // Heat transfer coefficient [W/m^2/K]
318  // This field is updated in Tcloth()
319  volScalarField hc
320  (
321  IOobject
322  (
323  "hc",
324  time_.name(),
325  mesh_
326  ),
327  mesh_,
328  dimensionedScalar(dimensionSet(1, 0, -3, -1, 0, 0, 0), 0)
329  );
330 
331  // Calculate the surface temperature of the cloth by an iterative
332  // process using equation (2) from DIN EN ISO 7730 [degC]
333  const volScalarField Tcloth
334  (
335  this->Tcloth
336  (
337  hc,
338  metabolicRateSI,
339  extWorkSI,
340  T,
341  Trad
342  )
343  );
344 
345  // Calculate the PMV quantity
346  const dimensionedScalar factor1(dimensionSet(-1, 0, 3, 0, 0, 0, 0), 0.303);
347  const dimensionedScalar factor2
348  (
349  dimless/metabolicRateSI.dimensions(),
350  -0.036
351  );
352  const dimensionedScalar factor3(factor1.dimensions(), 0.028);
353  const dimensionedScalar factor4(dimLength/dimTime, 3.05e-3);
354  const dimensionedScalar factor5(dimPressure, 5733);
355  const dimensionedScalar factor6(dimTime/dimLength, 6.99);
356  const dimensionedScalar factor8(metabolicRateSI.dimensions(), 58.15);
357  const dimensionedScalar factor9(dimless/dimPressure, 1.7e-5);
358  const dimensionedScalar factor10(dimPressure, 5867);
359  const dimensionedScalar factor11(dimless/dimTemperature, 0.0014);
360  const dimensionedScalar factor12(dimTemperature, 307.15);
361  const dimensionedScalar factor13
362  (
363  dimensionSet(1, 0, -3, -4, 0, 0, 0),
364  3.96e-8
365  );
366 
367  const scalar factor7
368  (
369  // Special treatment of Term4
370  // if metaRate - extWork < factor8, set to zero
371  (metabolicRateSI - extWorkSI).value() < factor8.value() ? 0 : 0.42
372  );
373 
374  Info<< "Calculating the predicted mean vote (PMV)\n";
375 
376  // Equation (1)
378  (
380  (
381  "PMV",
382 
383  // Term1: Thermal sensation transfer coefficient
384  (factor1*exp(factor2*metabolicRateSI) + factor3)
385  *(
386  (metabolicRateSI - extWorkSI)
387 
388  // Term2: Heat loss difference through skin
389  - factor4
390  *(
391  factor5
392  - factor6*(metabolicRateSI - extWorkSI)
393  - pSat*relHumidity_
394  )
395 
396  // Term3: Heat loss through sweating
397  - factor7*(metabolicRateSI - extWorkSI - factor8)
398 
399  // Term4: Heat loss through latent respiration
400  - factor9*metabolicRateSI*(factor10 - pSat*relHumidity_)
401 
402  // Term5: Heat loss through dry respiration
403  - factor11*metabolicRateSI*(factor12 - T)
404 
405  // Term6: Heat loss through radiation
406  - factor13*fcl_*(pow4(Tcloth) - pow4(Trad))
407 
408  // Term7: Heat loss through convection
409  - fcl_*hc*(Tcloth - T)
410  )
411  )
412  );
413 
414  Info<< "Calculating the predicted percentage of dissatisfaction (PPD)\n";
415 
416  // Equation (5)
418  (
420  (
421  "PPD",
422  100 - 95*exp(-0.03353*pow4(PMV()) - 0.21790*sqr(PMV()))
423  )
424  );
425 
426  Info<< "Calculating the draught rating (DR)\n";
427 
428  const dimensionedScalar Umin("Umin", dimVelocity, 0.05);
429  const dimensionedScalar Umax("Umax", dimVelocity, 0.5);
430  const dimensionedScalar pre("preU", dimless, 0.37);
431  const dimensionedScalar C1("C1", dimVelocity, 3.14);
432 
433  // Limit the velocity field to the values given in EN ISO 7733
434  volScalarField Umag(mag(lookupObject<volVectorField>("U")));
435  Umag.maxMin(Umin, Umax);
436 
437  // Calculate the turbulent intensity if turbulent kinetic energy field k
438  // exists
439  volScalarField TI
440  (
441  IOobject
442  (
443  "TI",
444  time_.name(),
445  mesh_
446  ),
447  mesh_,
448  dimensionedScalar(dimensionSet(0, 0, 0, 0, 0, 0, 0), 0)
449  );
450 
451  if (foundObject<volScalarField>("k"))
452  {
453  const volScalarField& k = lookupObject<volScalarField>("k");
454  TI = sqrt(2/3*k)/Umag;
455  }
456 
457  // For unit correctness
458  const dimensionedScalar correctUnit
459  (
460  "correctUnit",
461  dimensionSet(0,- 1.62, 1.62, -1, 0, 0, 0),
462  1
463  );
464 
465  // Equation (6)
467  (
469  (
470  "DR",
471  correctUnit*(factor12 - T)*pow(Umag - Umin, 0.62)*(pre*Umag*TI + C1)
472  )
473  );
474 
475  return
476  store(PMV)
477  && store(PPD)
478  && store(DR);
479 }
480 
481 
483 {
484  return
485  writeObject("PMV")
486  && writeObject("PPD")
487  && writeObject("DR");
488 }
489 
490 
491 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
label k
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Generic GeometricField class.
void maxMin(const dimensioned< Type > &minDt, const dimensioned< Type > &maxDt)
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
static tmp< GeometricField< Type, PatchField, GeoMesh > > New(const word &name, const Internal &, const PtrList< PatchField< Type >> &)
Return a temporary field constructed from name,.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:61
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:76
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
Dimension set for the base types.
Definition: dimensionSet.H:122
const dimensionSet & dimensions() const
Return const reference to dimensions.
const Type & value() const
Return const reference to value.
Abstract base-class for Time/database functionObjects.
Calculates the thermal comfort quantities predicted mean vote (PMV), predicted percentage of dissatis...
Definition: comfort.H:174
virtual wordList fields() const
Return the list of fields required.
Definition: comfort.C:301
comfort(const word &name, const Time &runTime, const dictionary &dict)
Construct from Time and dictionary.
Definition: comfort.C:230
virtual bool execute()
Calculate the predicted mean vote (PMV)
Definition: comfort.C:307
virtual bool write()
Write the PPD and PMV fields.
Definition: comfort.C:482
virtual bool read(const dictionary &)
Read the data needed for the comfort calculation.
Definition: comfort.C:262
virtual ~comfort()
Destructor.
Definition: comfort.C:256
Specialisation of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
const fvMesh & mesh_
Reference to the fvMesh.
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
A class for managing temporary objects.
Definition: tmp.H:55
A class for handling words, derived from string.
Definition: word.H:62
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
defineTypeNameAndDebug(adjustTimeStepToCombustion, 0)
addToRunTimeSelectionTable(functionObject, adjustTimeStepToCombustion, dictionary)
Namespace for OpenFOAM.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
Type gSum(const FieldField< Field, Type > &f)
const doubleScalar e
Definition: doubleScalar.H:105
dimensionedScalar exp(const dimensionedScalar &ds)
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
const dimensionSet dimPressure
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
dimensionedScalar pow3(const dimensionedScalar &ds)
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
const dimensionSet dimless
messageStream Info
const dimensionSet dimLength
const dimensionSet dimTemperature
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionSet dimTime
dimensionedScalar sqrt(const dimensionedScalar &ds)
VolField< scalar > volScalarField
Definition: volFieldsFwd.H:61
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensionedScalar pow4(const dimensionedScalar &ds)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
const dimensionSet dimMass
const dimensionSet dimVelocity
dimensionedScalar neg0(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
fvPatchField< scalar > fvPatchScalarField
dimensionedScalar pow025(const dimensionedScalar &ds)
dictionary dict