Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // Adapted from the code Version 1.1 released by Yi-Chuan Lu on February 23, 2023.
49 : //
50 : // @article{20heatindex,
51 : // Title = {Extending the Heat Index},
52 : // Author = {Yi-Chuan Lu and David M. Romps},
53 : // Journal = {Journal of Applied Meteorology and Climatology},
54 : // Year = {2022},
55 : // Volume = {61},
56 : // Number = {10},
57 : // Pages = {1367--1383},
58 : // Year = {2022},
59 : // }
60 : //
61 :
62 : #include <cmath>
63 :
64 : #include <EnergyPlus/Data/EnergyPlusData.hh>
65 : #include <EnergyPlus/ExtendedHeatIndex.hh>
66 : #include <EnergyPlus/General.hh>
67 : #include <EnergyPlus/HVACSystemRootFindingAlgorithm.hh>
68 :
69 : namespace EnergyPlus {
70 :
71 : namespace ExtendedHI {
72 :
73 : // The saturation vapor pressure
74 61878 : Real64 pvstar(Real64 const T)
75 : {
76 : // Thermodynamic parameters
77 61878 : constexpr Real64 Ttrip = 273.16; // K
78 61878 : constexpr Real64 ptrip = 611.65; // Pa
79 61878 : constexpr Real64 E0v = 2.3740e6; // J/kg
80 61878 : constexpr Real64 E0s = 0.3337e6; // J/kg
81 61878 : constexpr Real64 rgasa = 287.04; // J/kg/K
82 61878 : constexpr Real64 rgasv = 461.; // J/kg/K
83 61878 : constexpr Real64 cva = 719.; // J/kg/K
84 61878 : constexpr Real64 cvv = 1418.; // J/kg/K
85 61878 : constexpr Real64 cvl = 4119.; // J/kg/K
86 61878 : constexpr Real64 cvs = 1861.; // J/kg/K
87 61878 : constexpr Real64 cpa = cva + rgasa;
88 61878 : constexpr Real64 cpv = cvv + rgasv;
89 :
90 61878 : if (T == 0.0) {
91 278 : return 0.0;
92 61600 : } else if (T < Ttrip) {
93 2030 : return ptrip * pow((T / Ttrip), ((cpv - cvs) / rgasv)) * exp((E0v + E0s - (cvv - cvs) * Ttrip) / rgasv * (1. / Ttrip - 1. / T));
94 : } else {
95 59570 : return ptrip * pow((T / Ttrip), ((cpv - cvl) / rgasv)) * exp((E0v - (cvv - cvl) * Ttrip) / rgasv * (1. / Ttrip - 1. / T));
96 : }
97 : return 0.0;
98 : }
99 :
100 : // The latent heat of vaporization of water
101 47072 : Real64 Le(Real64 const T)
102 : {
103 : // Thermodynamic parameters
104 47072 : constexpr Real64 Ttrip = 273.16; // K
105 47072 : constexpr Real64 E0v = 2.3740e6; // J/kg
106 47072 : constexpr Real64 rgasv = 461.; // J/kg/K
107 47072 : constexpr Real64 cvv = 1418.; // J/kg/K
108 47072 : constexpr Real64 cvl = 4119.; // J/kg/K
109 47072 : return (E0v + (cvv - cvl) * (T - Ttrip) + rgasv * T);
110 : }
111 :
112 : // Function to calculate respiratory heat loss, W/m^2
113 47054 : Real64 Qv(Real64 const Ta, Real64 const Pa)
114 : {
115 47054 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
116 47054 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
117 47054 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
118 47054 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
119 47054 : constexpr Real64 r = 124.; // Pa/K , Zf/Rf, steadman1979
120 47054 : constexpr Real64 eta = 1.43e-6; // kg/J , "inhaled mass" / "metabolic rate", steadman1979
121 47054 : Real64 L = Le(310.); // , latent heat of vaporization at 310 K
122 47054 : constexpr Real64 p = 1.013e5; // Pa , atmospheric pressure
123 47054 : constexpr Real64 rgasa = 287.04; // J/kg/K
124 47054 : constexpr Real64 rgasv = 461.; // J/kg/K
125 47054 : constexpr Real64 cva = 719.; // J/kg/K
126 47054 : constexpr Real64 cvv = 1418.; // J/kg/K
127 47054 : constexpr Real64 cpa = cva + rgasa;
128 47054 : constexpr Real64 cpv = cvv + rgasv;
129 :
130 47054 : return eta * Q * (cpa * (Tc - Ta) + L * rgasa / (p * rgasv) * (Pc - Pa));
131 : }
132 :
133 : // Function to calculate mass transfer resistance through skin, Pa m^2/W
134 20943 : Real64 Zs(Real64 const Rs)
135 : {
136 20943 : return (Rs == 0.0387) ? 52.1 : 6.0e8 * std::pow(Rs, 5);
137 : }
138 :
139 : // Function to calculate heat transfer resistance through air, exposed part of skin, K m^2/W
140 31253 : Real64 Ra(Real64 const Ts, Real64 const Ta)
141 : {
142 31253 : constexpr Real64 sigma = 5.67e-8; // W/m^2/K^4 , Stefan-Boltzmann constant
143 31253 : constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979
144 31253 : constexpr Real64 hc = 17.4;
145 31253 : constexpr Real64 phi_rad = 0.85;
146 31253 : Real64 const hr = epsilon * phi_rad * sigma * (std::pow(Ts, 2) + std::pow(Ta, 2)) * (Ts + Ta);
147 31253 : return 1.0 / (hc + hr);
148 : }
149 :
150 : // Function to calculate heat transfer resistance through air, clothed part of skin, K m^2/W
151 44971 : Real64 Ra_bar(Real64 const Tf, Real64 const Ta)
152 : {
153 44971 : constexpr Real64 sigma = 5.67e-8; // W/m^2/K^4 , Stefan-Boltzmann constant
154 44971 : constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979
155 44971 : constexpr Real64 hc = 11.6;
156 44971 : constexpr Real64 phi_rad = 0.79;
157 44971 : Real64 const hr = epsilon * phi_rad * sigma * (std::pow(Tf, 2) + std::pow(Ta, 2)) * (Tf + Ta);
158 44971 : return 1.0 / (hc + hr);
159 : }
160 :
161 : // Function to calculate heat transfer resistance through air, when being naked, K m^2/W
162 22233 : Real64 Ra_un(Real64 const Ts, Real64 const Ta)
163 : {
164 22233 : constexpr Real64 sigma = 5.67e-8; // W/m^2/K^4 , Stefan-Boltzmann constant
165 22233 : constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979
166 22233 : constexpr Real64 hc = 12.3;
167 22233 : constexpr Real64 phi_rad = 0.80;
168 22233 : Real64 const hr = epsilon * phi_rad * sigma * (std::pow(Ts, 2) + std::pow(Ta, 2)) * (Ts + Ta);
169 22233 : return 1.0 / (hc + hr);
170 : }
171 :
172 : constexpr Real64 tol = 1e-5;
173 : constexpr Real64 maxIter = 100;
174 177 : Real64 find_eqvar_phi(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
175 : {
176 177 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
177 177 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
178 177 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
179 177 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
180 177 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
181 :
182 177 : Real64 phi = 0.84;
183 177 : Real64 const Pa = RH * pvstar(Ta);
184 177 : constexpr Real64 Rs = 0.0387;
185 177 : Real64 const ZsRs = Zs(Rs);
186 177 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
187 : Real64 Ts;
188 : int SolFla;
189 177 : General::SolveRoot(
190 : state,
191 : tol,
192 : maxIter,
193 : SolFla,
194 : Ts,
195 3445 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
196 177 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
197 177 : std::max(Tc, Ta) + Rs * std::abs(m));
198 177 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
199 177 : if (flux1 <= 0.0) {
200 159 : phi = 1.0 - (Q - Qv(Ta, Pa)) * Rs / (Tc - Ts);
201 : }
202 177 : return phi;
203 : }
204 :
205 604 : Real64 find_eqvar_Rf(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
206 : {
207 604 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
208 604 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
209 604 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
210 604 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
211 604 : constexpr Real64 r = 124.; // Pa/K , Zf/Rf, steadman1979
212 604 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
213 604 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
214 :
215 604 : Real64 Pa = RH * pvstar(Ta);
216 604 : constexpr Real64 Rs = 0.0387;
217 604 : constexpr Real64 phi = 0.84;
218 604 : Real64 const ZsRs = Zs(Rs);
219 604 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
220 604 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
221 : Real64 Ts;
222 : int SolFla;
223 604 : General::SolveRoot(
224 : state,
225 : tol,
226 : maxIter,
227 : SolFla,
228 : Ts,
229 8649 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
230 604 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
231 604 : std::max(Tc, Ta) + Rs * std::abs(m));
232 : Real64 Tf;
233 604 : General::SolveRoot(
234 : state,
235 : tol,
236 : maxIter,
237 : SolFla,
238 : Tf,
239 8604 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
240 604 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
241 604 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
242 604 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
243 604 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
244 : Real64 Rf;
245 604 : if (flux1 <= 0.0) {
246 24 : Rf = std::numeric_limits<Real64>::infinity();
247 580 : } else if (flux2 <= 0.0) {
248 555 : Real64 const Ts_bar = Tc - (Q - Qv(Ta, Pa)) * Rs / phi + (1.0 / phi - 1.0) * (Tc - Ts);
249 555 : General::SolveRoot(
250 : state,
251 : tol,
252 : maxIter,
253 : SolFla,
254 : Tf,
255 555 : [&](Real64 Tf) {
256 7982 : return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) * (Tf - Ta) / ((ZsRs + Za_bar) * (Tf - Ta) + r * Ra_bar(Tf, Ta) * (Ts_bar - Tf)) -
257 7982 : (Tc - Ts_bar) / Rs;
258 : },
259 : Ta,
260 : Ts_bar);
261 555 : Rf = Ra_bar(Tf, Ta) * (Ts_bar - Tf) / (Tf - Ta);
262 : } else {
263 25 : Rf = 0.0;
264 : }
265 604 : return Rf;
266 : }
267 :
268 1206 : Real64 find_eqvar_rs(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
269 : {
270 1206 : constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018
271 1206 : constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018
272 1206 : Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014
273 1206 : constexpr Real64 cpc = 3492.; // J/kg/K , specific heat capacity of core, gagge1972
274 1206 : Real64 C = M * cpc / A; // , heat capacity of core
275 1206 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
276 1206 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
277 1206 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
278 1206 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
279 1206 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
280 1206 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
281 1206 : constexpr Real64 Za_un = 60.6 / 12.3; // Pa m^2/W, mass transfer resistance through air, when being naked
282 :
283 1206 : Real64 Pa = RH * pvstar(Ta);
284 1206 : constexpr Real64 phi = 0.84;
285 1206 : Real64 Rs = 0.0387;
286 1206 : Real64 ZsRs = Zs(Rs);
287 1206 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
288 1206 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
289 : Real64 Ts;
290 : int SolFla;
291 1206 : General::SolveRoot(
292 : state,
293 : tol,
294 : maxIter,
295 : SolFla,
296 : Ts,
297 11252 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
298 1206 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
299 1206 : std::max(Tc, Ta) + Rs * std::abs(m));
300 : Real64 Tf;
301 1206 : General::SolveRoot(
302 : state,
303 : tol,
304 : maxIter,
305 : SolFla,
306 : Tf,
307 10802 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
308 1206 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
309 1206 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
310 1206 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
311 1206 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
312 1206 : Real64 const flux3 = Q - Qv(Ta, Pa) - (Tc - Ta) / Ra_un(Tc, Ta) - (phi_salt * pvstar(Tc) - Pa) / Za_un;
313 1206 : if (flux1 > 0 && flux2 > 0) {
314 1047 : if (flux3 < 0.0) {
315 870 : General::SolveRoot(
316 : state,
317 : tol,
318 : maxIter,
319 : SolFla,
320 : Ts,
321 15308 : [&](Real64 Ts) { return (Ts - Ta) / Ra_un(Ts, Ta) + (Pc - Pa) / (Zs((Tc - Ts) / (Q - Qv(Ta, Pa))) + Za_un) - (Q - Qv(Ta, Pa)); },
322 : 0.0,
323 : Tc);
324 870 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
325 870 : ZsRs = Zs(Rs);
326 870 : Real64 Ps = Pc - (Pc - Pa) * ZsRs / (ZsRs + Za_un);
327 870 : if (Ps > phi_salt * pvstar(Ts)) {
328 198 : General::SolveRoot(
329 : state,
330 : tol,
331 : maxIter,
332 : SolFla,
333 : Ts,
334 2608 : [&](Real64 Ts) { return (Ts - Ta) / Ra_un(Ts, Ta) + (phi_salt * pvstar(Ts) - Pa) / Za_un - (Q - Qv(Ta, Pa)); },
335 : 0.0,
336 : Tc);
337 198 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
338 : }
339 : } else {
340 177 : Rs = 0.0;
341 : }
342 : }
343 1206 : return Rs;
344 : }
345 :
346 703 : Real64 find_eqvar_dTcdt(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
347 : {
348 703 : constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018
349 703 : constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018
350 703 : Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014
351 703 : constexpr Real64 cpc = 3492.; // J/kg/K , specific heat capacity of core, gagge1972
352 703 : Real64 C = M * cpc / A; // , heat capacity of core
353 703 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
354 703 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
355 703 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
356 703 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
357 703 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
358 703 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
359 703 : constexpr Real64 Za_un = 60.6 / 12.3; // Pa m^2/W, mass transfer resistance through air, when being naked
360 :
361 703 : Real64 dTcdt = 0.0;
362 703 : Real64 const Pa = RH * pvstar(Ta);
363 703 : constexpr Real64 Rs = 0.0387;
364 703 : Real64 const ZsRs = Zs(Rs);
365 703 : constexpr Real64 phi = 0.84;
366 703 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
367 703 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
368 : Real64 Ts;
369 : int SolFla;
370 703 : General::SolveRoot(
371 : state,
372 : tol,
373 : maxIter,
374 : SolFla,
375 : Ts,
376 8376 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
377 703 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
378 703 : std::max(Tc, Ta) + Rs * std::abs(m));
379 : Real64 Tf;
380 703 : General::SolveRoot(
381 : state,
382 : tol,
383 : maxIter,
384 : SolFla,
385 : Tf,
386 8533 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
387 703 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
388 703 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
389 703 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
390 703 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
391 703 : Real64 const flux3 = Q - Qv(Ta, Pa) - (Tc - Ta) / Ra_un(Tc, Ta) - (phi_salt * pvstar(Tc) - Pa) / Za_un;
392 703 : if (flux1 > 0.0 && flux2 > 0.0 && flux3 >= 0.0) {
393 599 : dTcdt = (1.0 / C) * flux3;
394 : }
395 703 : return dTcdt;
396 : }
397 :
398 : // given T and RH, returns a key and value pair
399 304 : Real64 find_eqvar_name_and_value(EnergyPlusData &state, Real64 const Ta, Real64 const RH, EqvarName &varname)
400 : {
401 304 : constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018
402 304 : constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018
403 304 : Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014
404 304 : constexpr Real64 cpc = 3492.; // J/kg/K , specific heat capacity of core, gagge1972
405 304 : Real64 C = M * cpc / A; // , heat capacity of core
406 304 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
407 304 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
408 304 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
409 304 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
410 304 : constexpr Real64 r = 124.; // Pa/K , Zf/Rf, steadman1979
411 304 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
412 304 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
413 304 : constexpr Real64 Za_un = 60.6 / 12.3; // Pa m^2/W, mass transfer resistance through air, when being naked
414 :
415 304 : Real64 Pa = RH * pvstar(Ta);
416 304 : Real64 Rs = 0.0387;
417 304 : Real64 phi = 0.84;
418 304 : Real64 dTcdt = 0.0;
419 304 : Real64 ZsRs = Zs(Rs);
420 304 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
421 304 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
422 :
423 : int SolFla;
424 : Real64 Ts;
425 304 : General::SolveRoot(
426 : state,
427 : tol,
428 : maxIter,
429 : SolFla,
430 : Ts,
431 2476 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
432 304 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
433 304 : std::max(Tc, Ta) + Rs * std::abs(m));
434 :
435 : Real64 Tf;
436 304 : General::SolveRoot(
437 : state,
438 : tol,
439 : maxIter,
440 : SolFla,
441 : Tf,
442 2463 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
443 304 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
444 304 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
445 304 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
446 304 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
447 : Real64 Rf;
448 :
449 304 : if (flux1 <= 0.0) {
450 12 : varname = EqvarName::Phi;
451 12 : phi = 1.0 - (Q - Qv(Ta, Pa)) * Rs / (Tc - Ts);
452 12 : return phi;
453 292 : } else if (flux2 <= 0.0) {
454 38 : varname = EqvarName::Rf;
455 38 : Real64 const Ts_bar = Tc - (Q - Qv(Ta, Pa)) * Rs / phi + (1.0 / phi - 1.0) * (Tc - Ts);
456 38 : General::SolveRoot(
457 : state,
458 : tol,
459 : maxIter,
460 : SolFla,
461 : Tf,
462 38 : [&](Real64 Tf) {
463 390 : return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) * (Tf - Ta) / ((ZsRs + Za_bar) * (Tf - Ta) + r * Ra_bar(Tf, Ta) * (Ts_bar - Tf)) -
464 390 : (Tc - Ts_bar) / Rs;
465 : },
466 : Ta,
467 : Ts_bar);
468 38 : Rf = Ra_bar(Tf, Ta) * (Ts_bar - Tf) / (Tf - Ta);
469 38 : return Rf;
470 : } else {
471 254 : Real64 const flux3 = Q - Qv(Ta, Pa) - (Tc - Ta) / Ra_un(Tc, Ta) - (phi_salt * pvstar(Tc) - Pa) / Za_un;
472 254 : if (flux3 < 0.0) {
473 163 : varname = EqvarName::Rs;
474 163 : General::SolveRoot(
475 : state,
476 : tol,
477 : maxIter,
478 : SolFla,
479 : Ts,
480 2637 : [&](Real64 Ts) { return (Ts - Ta) / Ra_un(Ts, Ta) + (Pc - Pa) / (Zs((Tc - Ts) / (Q - Qv(Ta, Pa))) + Za_un) - (Q - Qv(Ta, Pa)); },
481 : 0.0,
482 : Tc);
483 163 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
484 163 : ZsRs = Zs(Rs);
485 163 : Real64 const Ps = Pc - (Pc - Pa) * ZsRs / (ZsRs + Za_un);
486 163 : if (Ps > phi_salt * pvstar(Ts)) {
487 62 : General::SolveRoot(
488 : state,
489 : tol,
490 : maxIter,
491 : SolFla,
492 : Ts,
493 761 : [&](Real64 Ts) { return (Ts - Ta) / Ra_un(Ts, Ta) + (phi_salt * pvstar(Ts) - Pa) / Za_un - (Q - Qv(Ta, Pa)); },
494 : 0.0,
495 : Tc);
496 62 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
497 : }
498 163 : return Rs;
499 : } else {
500 91 : varname = EqvarName::DTcdt;
501 91 : Rs = 0.0;
502 91 : dTcdt = (1.0 / C) * flux3;
503 91 : return dTcdt;
504 : }
505 : }
506 : }
507 :
508 : // Convert the find_T function
509 281 : Real64 find_T(EnergyPlusData &state, EqvarName const eqvar_name, Real64 const eqvar)
510 : {
511 : Real64 T;
512 : int SolFla;
513 281 : constexpr Real64 Pa0 = 1.6e3; // Pa , reference air vapor pressure in regions III, IV, V, VI, steadman1979
514 :
515 281 : if (eqvar_name == EqvarName::Phi) {
516 18 : General::SolveRoot(
517 213 : state, tol, maxIter, SolFla, T, [&](Real64 T) { return find_eqvar_phi(state, T, 1.0) - eqvar; }, 0.0, 240.0);
518 263 : } else if (eqvar_name == EqvarName::Rf) {
519 24 : General::SolveRoot(
520 : state,
521 : tol,
522 : maxIter,
523 : SolFla,
524 : T,
525 652 : [&](Real64 T) { return (find_eqvar_Rf(state, T, std::min(1.0, Pa0 / pvstar(T)))) - eqvar; },
526 : 230.0,
527 : 300.0);
528 239 : } else if (eqvar_name == EqvarName::Rs) {
529 157 : General::SolveRoot(
530 1520 : state, tol, maxIter, SolFla, T, [&](Real64 T) { return find_eqvar_rs(state, T, Pa0 / pvstar(T)) - eqvar; }, 295.0, 350.0);
531 : } else {
532 82 : General::SolveRoot(
533 867 : state, tol, maxIter, SolFla, T, [&](Real64 T) { return find_eqvar_dTcdt(state, T, Pa0 / pvstar(T)) - eqvar; }, 340.0, 1000.0);
534 : }
535 :
536 281 : return T;
537 : }
538 :
539 262 : Real64 heatindex(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
540 : {
541 :
542 : // Ta: temperature in Kelvin
543 : // RH: relative humidity in range of 0.0 to 1.0
544 : // The function computes the extended heat index, in Kelvinn
545 :
546 262 : auto const HVACSystemRootSolverMethodBackup = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod;
547 262 : state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverAlgorithm::ShortBisectionThenRegulaFalsi;
548 262 : EqvarName eqvar_name = EqvarName::Invalid;
549 262 : Real64 const eqvar_value = find_eqvar_name_and_value(state, Ta, RH, eqvar_name);
550 :
551 262 : Real64 T = find_T(state, eqvar_name, eqvar_value);
552 :
553 262 : if (Ta == 0.0) T = 0.0;
554 :
555 262 : state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverMethodBackup;
556 262 : return T;
557 : }
558 :
559 : } // namespace ExtendedHI
560 : } // namespace EnergyPlus
|
