comfortFoam.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 held by original author
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 the
13  Free Software Foundation; either version 2 of the License, or (at your
14  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, write to the Free Software Foundation,
23  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
24 
25 Application
26  comfortFoam
27 
28 Description
29  Das Comfort-Tool berechnet folgende Kennwerte:
30  - Operative Raumlufttemperatur (TOp)
31  - DR (Draugth Rating oder Draft Risk), Wertebereich: 0 bis 100% (DR)
32  - PMV (Predicted Percentage of Dissatisfied), Wertebereich: 0 bis 100% (PMV)
33  - PPD (Predicted Mean Vote), Wertebereich: -3 (kalt) bis 3 (warm) (PPD)
34  - Mittlere Verweildauer der Luft - Mean Age of Air (AoA)
35  - Lüftungseffektivität (AE)
36  - Berechnet die relative Raumluftfeuchtigkeit (Humidity), Wertebereich: 0 bis 100%
37 
38  Grundlage: "DIN EN ISO 7730"
39  Autor: "Dipl.-Ing. Thomas Tian"
40  Version: "1.5"
41 
42 Source files:
43  createFields.H
44 
45 \*---------------------------------------------------------------------------*/
46 
47 #include "fvCFD.H"
49 #include "wallFvPatch.H"
50 
51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
52 
54 {
55  const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
56  vector midU(0,0,0);
57  scalar cellsum(0);
58  scalar volume(0);
59 
60  forAll (mesh_.cells(), cellI)
61  {
62  midU += U[cellI] * mesh_.V()[cellI];
63  volume += mesh_.V()[cellI];
64  };
65 
66  Info<< "Average velocity = "<< mag(midU/volume) << " m/s" << nl;
67 
68  return midU/(volume);
69 }
70 
71 Foam::scalar Taverage(const Foam::fvMesh& mesh_)
72 {
73  const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
74  scalar midT(0);
75  scalar cellsum(0);
76  scalar volume(0);
77 
78  forAll (mesh_.cells(), cellI)
79  {
80  midT += T[cellI] * mesh_.V()[cellI];
81  volume += mesh_.V()[cellI];
82  };
83  Info<< "Average Temperature = "<< (midT/volume)-273.15 << " °C" << nl;
84 
85  return (midT/volume);
86 }
87 
88 Foam::scalar radiationTemperature(const Foam::fvMesh& mesh_, const Foam::fvPatchList& Patches_)
89 {
90  // Berechnet die Strahlungstemperatur der Umgebung
91  // radiationTemperatur = Sum(Mittlere Wandoberflächentemperaturen / Anzahl Flächen)
92  // Rückgabe: Temperatur in Kelvin
93 
94  const volScalarField& T = mesh_.lookupObject<volScalarField>("T");
95  scalar PatchSumTemp(0), area(0);
96  int counter = 0;
97 
98  forAll (Patches_, patchI)
99  {
100 
101  const label curPatch = Patches_[patchI].index();
102 
103  if (isType<wallFvPatch>( Patches_[patchI] ))
104  {
105 
106  area = gSum(mesh_.magSf().boundaryField()[curPatch]);
107 
108  if (area != 0)
109  {
110  PatchSumTemp +=
111  gSum
112  (
113  mesh_.magSf().boundaryField()[curPatch]
114  * T.boundaryField()[curPatch]
115  ) / area;
116 
117  counter++;
118  }
119  }
120  }
121  return PatchSumTemp / counter - 273.15;
122 }
123 
124 int main(int argc, char *argv[])
125 {
126 
127  argList::noParallel();
128 
129 # include "addTimeOptions.H"
130 # include "setRootCase.H"
131 # include "createTime.H"
132 
133  // Get times list
134  instantList Times = runTime.times();
135 
136  // set startTime and endTime depending on -time and -latestTime options
137 # include "checkTimeOptions.H"
138 
139  runTime.setTime(Times[startTime], startTime);
141 
142 # include "createMesh.H"
143 
144  forAll(timeDirs, timeI)
145  {
146  runTime.setTime(timeDirs[timeI], timeI);
147 
148  Info<< "Time = " << runTime.timeName() << endl;
149 
150  #include "createFields.H"
151 
152  // Operative Raumluftemperatur Top
154  const fvPatchList& Patches = T.mesh().boundary();
155 
156  scalar STemp(20);
157 
158  if ((w.headerOk()==1) )
159  {
160  RHswitch = true;
161  }
162 
163  if (QrHeader.headerOk()==1 )
164  {
166  const fvPatchList& Patches_ = Qr.mesh().boundary();
167  scalar PatchSumTemp(0);
168 
169  forAll (Patches_, patchI)
170  {
171  const label curPatch = Patches_[patchI].index();
172 
173  if (isType<wallFvPatch>( Patches_[patchI] ))
174  {
175  PatchSumTemp +=
176  gSum
177  (
178  mesh.magSf().boundaryField()[curPatch]
179  * Qr.boundaryField()[curPatch]
180  );
181  }
182  }
183  Info<< "Sum of all heat flows: "<< PatchSumTemp-273.15 << " °C" << endl;
184  }
185  else
186  {
187  STemp = radiationTemperature(mesh, Patches);
188  Info << "Average Radiation Temperature: " << STemp << " °C" << endl;
189  };
190 
191  if (AoAHeader.headerOk()==1 )
192  {
193 
195  scalar cellsum(0);
196  scalar midAoA(0);
197 
198  forAll (mesh.cells(), cellI)
199  {
200  midAoA += AoA[cellI];
201  cellsum++;
202  };
203 
204  scalar maxValue = max(AoA).value();
205 
206  Info << "Average age of air " << (midAoA/cellsum) << " s" << endl;
207  Info << "Maximum age of air " << maxValue << " s" << endl;
208 
209 
210  // Lüftungseffektivität AE = (Tn/2TM) x 100 %
211  Info << "Room volume " << gSum( mesh.V() ) << " m3" << endl;
212  Info << "Ventilation efficiency AE = " << ( gSum( mesh.V() ) / sumZuluft ) / (2 * ( maxValue / 3600) ) << " %" << endl;
213  };
214 
216  vector wl = Uaverage(mesh);
217  scalar Tr = Taverage(mesh);
219 
220  scalar cellsum(0);
221  scalar P1(0), P2(0), P3(0), P4(0), P5(0), XN(0), XF(0), HCN(0), HCF(0), HC(0), PA(0), FCL(0), EPS(0), ICL(0);
222  scalar Tu(0), HL1(0), HL2(0), HL3(0), HL4(0), HL5(0), HL6(0), TCL(0), TS(0), TCLA(0), midPMV(0), midPPD(0), midDR(0), midTO(0), midRH(0);
223  scalar p_w(0), p_ws(0.01);
224 
225  scalar volume(0);
226 
227  int N;
228  double a_korr; // Korrekturfaktor a nach DIN EN 7730;
229 
230  // FANGER - Gleichungen (PMV)
231 
232  forAll (mesh.cells(), cellI)
233  {
234 
235  // Zum Testen
236  /*
237  STemp= 21;
238  Tr = 23.0+273.15;
239  T[cellI] = Tr;
240  U[cellI].x() = 0.3;
241  U[cellI].y() = 0;
242  U[cellI].z() = 0;
243  */
244 
245  if (T[cellI] < 0)
246  T[cellI] = 273;
247 
248  if (T[cellI] > 350)
249  T[cellI] = 350;
250 
251 
252  if (RHswitch)
253  {
254  // Berechne die Raumluftfeuchte
255  p_w = ((101325 + p_rgh[cellI] ) * w[cellI])/(0.62198 + 0.37802 * w[cellI]);
256 
257  p_ws = Foam::exp(-0.58002206*Foam::pow(10.0,4)*Foam::pow(T[cellI],-1)
258  +0.13914993*10
259  -0.48640239*Foam::pow(10.0,-1)*T[cellI]
260  +0.41764768*Foam::pow(10.0,-4)*Foam::pow(T[cellI],2)
261  -0.14452093*Foam::pow(10.0,-7)*Foam::pow(T[cellI],3)
262  +0.65459673*10*Foam::log(T[cellI]) );
263 
264  RH[cellI] = p_w/p_ws*100;
265  PA = RH[cellI] * 10 * Foam::exp( 16.6563 - (4030.183 / (T[cellI] - 273.15 + 235) )); // water vapour pressure, Pa
266  midRH += RH[cellI];
267  } else {
268  PA = RH1 * 10 * Foam::exp( 16.6563 - (4030.183 / (T[cellI] - 273.15 + 235) )); // water vapour pressure, Pa
269  midRH += RH1;
270  }
271 
272  // Berechne den Turbulenzgrad %
273  // Tu = Wurzel ( 1/3 * (Ux'² + Uy'² + Uz'²) ) / Mittelwert U
274  if ( mag(wl) >0 )
275  {
276  Tu = (Foam::sqrt(scalar(1)/scalar(3) * (
277  Foam::pow( U[cellI].x() ,2) +
278  Foam::pow( U[cellI].y() ,2) +
279  Foam::pow( U[cellI].z() ,2)
280  )) / mag(wl)) * 100;
281  // ToDo!!!!
282  Tu = 40;
283 
284  }
285  else
286  {
287  Tu = 0;
288  };
289 
290  // Berechne DR-Wert
291  // bei wl<0,05 m/s ist wl=0,05 m/s einzusetzen!
292  if (mag(U[cellI]) >= 0.05)
293  {
294  DR[cellI] = ( 34 - (T[cellI] - 273.15) ) * ( Foam::pow( mag(U[cellI]) - 0.05 ,0.62 ) * ( (0.37 * mag(U[cellI]) * Tu) + 3.14) );
295  }
296  else
297  {
298  DR[cellI] = (34 - T[cellI] - 273.15) * Foam::pow( 0.05 ,0.6223) * ((0.37 * 0.05 * Tu) + 3.14);
299  };
300  if (DR[cellI] > 100)
301  {
302  DR[cellI] = 100;
303  };
304  if (DR[cellI] < 0)
305  {
306  DR[cellI] = 0;
307  };
308 
309 
310 
311  ICL = 0.155 * clo;
312 
313  if (ICL < 0.078) // ok
314  { FCL = 1 + 1.29 * ICL; } // ok
315  else
316  { FCL = 1.05 + 0.645 * ICL; }; // clothing area factor
317 
318  HCF = 12.1 * Foam::sqrt(mag( U[cellI] )); // heat transf. coeff. by forced convection
319 
320  TCLA = T[cellI] + (35.5 - (T[cellI] - 273.15)) / (3.5 * 6.45 * ( ICL + 0.1));
321  XN = TCLA / 100; // ok
322 
323  P1 = ICL * FCL; // ok
324  P2 = P1 * 3.96; // ok
325  P3 = P1 * 100; // ok
326  P4 = P1 * T[cellI]; // ok
327  P5 = 308.7 - 0.028 * (met*58.15 - wme*58.15) + P2 * Foam::pow( (STemp + 273.0)/100, 4);
328 
329  // geändert Tian 06.10.2012
330  // XF = XN;
331  XF = TCLA / 50;
332  EPS = 0.0015;
333 
334  N=0;
335 
336  do
337  {
338 
339  N++;
340  XF = (XF + XN)/2; // stop criteria by iteration
341  HCF = 12.1 * Foam::sqrt(mag( U[cellI] ));
342  HCN = 2.38 * Foam::pow( mag( 100 * XF - T[cellI] ), 0.25); // heat transf. coeff. by natural convection;
343 
344  if (HCF>HCN) { HC = HCF; } else { HC = HCN; }; // ok
345  XN = (P5 + P4 * HC - P2 * Foam::pow(XF,4.0) ) / (100 + P3 * HC);
346  if (N>150) { break; };
347  } while (mag(XN-XF)>EPS);
348 
349 
350  TCL = 100 * XN - 273; // surface temperature of clothing
351  HL1 = 3.05 * 0.001 * (5733 - (6.99 * (met*58.15 - wme*58.15)) - PA); // heat loss diff. through skin
352  if ( (met*58.15 - wme*58.15) > 58.15) { HL2 = 0.42 * ( (met*58.15 - wme*58.15) - 58.15); }
353  else { HL2 = 0; }; // heat loss by sweating (comfort)
354  HL3 = 1.7 * 0.00001 * met * 58.15 * (5867 - PA); // latent respiration heat loss
355  HL4 = 0.0014 * met * 58.15 * (34 - (T[cellI] - 273.15) ); // dry respiration heat loss
356  HL5 = 3.96 * FCL * (Foam::pow(XN,4) - Foam::pow( ((STemp + 273.0)/100),4)) ; // heat lose by radiation
357  HL6 = FCL * HC * (TCL - (T[cellI] - 273.15)); // heat lose by convection
358  TS = 0.303 * Foam::exp(-0.036 * met * 58.15) + 0.028; // thermal sensation trans coeff
359 
360  // Berechne PMV-Wert
361  PMV[cellI] = TS * ( (met*58.15 - wme*58.15) - HL1 - HL2 - HL3 - HL4 - HL5 - HL6 );
362 
363  // Berechne PPD-Wert
364  PPD[cellI] = 100 - 95 * Foam::exp( -0.03353*pow(PMV[cellI],4) - 0.2179 * pow(PMV[cellI],2) );
365 
366  if (mag(U[cellI])<0.2) { a_korr=0.5; };
367  if ( (mag(U[cellI])>=0.2) || (mag(U[cellI])<=0.6) ) { a_korr=0.6; };
368  if (mag(U[cellI])>0.6) { a_korr=0.7; };
369  // Berechne die operative Raumlufttemperatur
370  TOp[cellI] = a_korr * T[cellI] + (1 - a_korr) * (STemp+273.15);
371 
372  midPMV += PMV[cellI];
373  midPPD += PPD[cellI];
374  midDR += DR[cellI];
375  midTO += TOp[cellI];
376 
377  volume += mesh.V()[cellI];
378  cellsum++;
379 
380 
381  };
382 
383  DR.write();
384  PMV.write();
385  PPD.write();
386  TOp.write();
387  RH.write();
388 
389  Info<< "Mean Radiation temperature "<< STemp << " °C" << endl;
390  Info<< "Water vapour pressure "<< PA << " Pa" <<nl;
391  Info << "Average PMV-Value = " << (midPMV/cellsum) << nl;
392  Info << "Average PPD-Value = " << (midPPD/cellsum) << " %" << nl;
393  Info << "Average DR-Value = " << (midDR/cellsum) << " %" << nl;
394  Info<< "Average air velocity = "<< mag(wl) << " m/s" << nl;
395  Info<< "Average room temperature = "<< mag(Tr)-273.15 << " °C" << nl;
396  Info<< "Average relative room humidity = "<< mag(midRH/cellsum) << " %" << nl;
397  Info << "Average operative temperature = " << (midTO/cellsum)-273.15 << " °C" << endl;
398 
399  if ( ((midPMV/cellsum > -0.2) || (midPMV/cellsum < 0.2)) && (midDR/cellsum < 10) && (midPPD/cellsum < 6))
400  {
401  Info << "Analysis: Category A" << endl;
402  } else {
403  if ( ((midPMV/cellsum > -0.5) || (midPMV/cellsum < 0.5)) && (midDR/cellsum < 20) && (midPPD/cellsum < 10))
404  {
405  Info << "Analysis: Category B" << endl;
406  } else {
407  if ( ((midPMV/cellsum > -0.7) || (midPMV/cellsum < 0.7)) && (midDR/cellsum < 30) && (midPPD/cellsum < 15))
408  {
409  Info << "Analysis: Category C" << endl;
410  } else
411  {
412  Info << "Analysis: Category D" << endl;
413  }}};
414 
415  // Info<< "Average Ux = "<< midUx / cellsum << endl;
416  // Info<< "Average Uy = "<< midUy / cellsum << endl;
417  // Info<< "Average Uz = "<< midUz / cellsum << endl;
418  }
419 
420  Info<< "\nEnd\n" << endl;
421 
422  return(0);
423 }
424 
425 
426 // ************************************************************************* //
427 
w
volScalarField w(IOobject("w", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE), mesh, dimensionedScalar("w", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0))
PMV
volScalarField PMV(IOobject("PMV", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar("PMV", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0))
RHswitch
bool RHswitch
Definition: createFields.H:177
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
singlePhaseTransportModel.H
Uaverage
Foam::vector Uaverage(const Foam::fvMesh &mesh_)
Definition: comfortFoam.C:64
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:136
AoAHeader
IOobject AoAHeader("AoA", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)
wallFvPatch.H
main
int main(int argc, char *argv[])
Definition: comfortFoam.C:149
Foam::GeometricField::boundaryField
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
Definition: GeometricField.C:735
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Foam::dimensioned::value
const Type & value() const
Return const reference to value.
Definition: dimensionedType.C:261
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:564
Foam::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:252
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
addTimeOptions.H
wme
const scalar wme(readScalar(comfortFoamDict.lookup("wme")))
U
U
Definition: pEqn.H:46
Foam::fvMesh::magSf
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
Definition: fvMeshGeometry.C:358
Foam::IOobject::headerOk
bool headerOk()
Read and check header info.
Definition: IOobject.C:439
UHeader
IOobject UHeader("U", runTime.timeName(), mesh, IOobject::MUST_READ)
Foam::fvMesh::V
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
Definition: fvMeshGeometry.C:199
Foam::label
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
RH
volScalarField RH(IOobject("RH", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar("RH", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0))
Foam::nl
static const char nl
Definition: Ostream.H:260
Foam::Info
messageStream Info
RH1
const scalar RH1(readScalar(comfortFoamDict.lookup("RH")))
Taverage
Foam::scalar Taverage(const Foam::fvMesh &mesh_)
Definition: comfortFoam.C:87
THeader
IOobject THeader("T", runTime.timeName(), mesh, IOobject::MUST_READ)
radiationTemperature
Foam::scalar radiationTemperature(const Foam::fvMesh &mesh_, const Foam::fvPatchList &Patches_)
Definition: comfortFoam.C:108
Qr
volScalarField Qr(IOobject("Qr", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE), mesh, dimensionedScalar("Qr", dimMass/pow3(dimTime), 0.0))
Foam::PtrList
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:61
p_rghHeader
IOobject p_rghHeader("p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ)
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:73
Foam::log
dimensionedScalar log(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:253
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
met
const scalar met(readScalar(comfortFoamDict.lookup("met")))
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
checkTimeOptions.H
clo
const scalar clo(readScalar(comfortFoamDict.lookup("clo")))
Foam::max
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
setRootCase.H
timeDirs
static instantList timeDirs
Definition: globalFoam.H:44
createFields.H
T
const volScalarField & T
Definition: createFields.H:25
Foam::Vector
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:57
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:142
createMesh.H
createTime.H
x
x
Definition: LISASMDCalcMethod2.H:52
TOp
volScalarField TOp(IOobject("TOp", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar("TOp", dimensionSet(0, 0, 0, 1, 0, 0, 0), 0.0))
fvCFD.H
startTime
Foam::label startTime
Definition: checkTimeOptions.H:5
DR
volScalarField DR(IOobject("DR", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar("DR", dimensionSet(0, 1,-1, 0, 0, 0, 0), 0.0))
Foam::objectRegistry::lookupObject
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
Definition: objectRegistryTemplates.C:165
Foam::timeSelector::select0
static instantList select0(Time &runTime, const argList &args)
Return the set of times selected based on the argList options.
Definition: timeSelector.C:253
Foam::GeometricField
Generic GeometricField class.
Definition: surfaceFieldsFwd.H:52
args
Foam::argList args(argc, argv)
p_rgh
volScalarField & p_rgh
Definition: setRegionFluidFields.H:31
sumZuluft
const scalar sumZuluft(readScalar(comfortFoamDict.lookup("Zuluft")))
y
scalar y
Definition: LISASMDCalcMethod1.H:14
maxValue
scalar maxValue
Definition: LISASMDCalcMethod1.H:5
PPD
volScalarField PPD(IOobject("PPD", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar("PPD", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0))
QrHeader
IOobject QrHeader("Qr", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE)
N
label N
Definition: createFields.H:22