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 0 : Real64 pvstar(Real64 const T)
75 : {
76 : // Thermodynamic parameters
77 0 : constexpr Real64 Ttrip = 273.16; // K
78 0 : constexpr Real64 ptrip = 611.65; // Pa
79 0 : constexpr Real64 E0v = 2.3740e6; // J/kg
80 0 : constexpr Real64 E0s = 0.3337e6; // J/kg
81 0 : constexpr Real64 rgasa = 287.04; // J/kg/K
82 0 : constexpr Real64 rgasv = 461.; // J/kg/K
83 0 : constexpr Real64 cva = 719.; // J/kg/K
84 0 : constexpr Real64 cvv = 1418.; // J/kg/K
85 0 : constexpr Real64 cvl = 4119.; // J/kg/K
86 0 : constexpr Real64 cvs = 1861.; // J/kg/K
87 0 : constexpr Real64 cpa = cva + rgasa;
88 0 : constexpr Real64 cpv = cvv + rgasv;
89 :
90 0 : if (T == 0.0) {
91 0 : return 0.0;
92 0 : } else if (T < Ttrip) {
93 0 : return ptrip * pow((T / Ttrip), ((cpv - cvs) / rgasv)) * exp((E0v + E0s - (cvv - cvs) * Ttrip) / rgasv * (1. / Ttrip - 1. / T));
94 : } else {
95 0 : 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 0 : Real64 Le(Real64 const T)
102 : {
103 : // Thermodynamic parameters
104 0 : constexpr Real64 Ttrip = 273.16; // K
105 0 : constexpr Real64 E0v = 2.3740e6; // J/kg
106 0 : constexpr Real64 rgasv = 461.; // J/kg/K
107 0 : constexpr Real64 cvv = 1418.; // J/kg/K
108 0 : constexpr Real64 cvl = 4119.; // J/kg/K
109 0 : return (E0v + (cvv - cvl) * (T - Ttrip) + rgasv * T);
110 : }
111 :
112 : // Function to calculate respiratory heat loss, W/m^2
113 0 : Real64 Qv(Real64 const Ta, Real64 const Pa)
114 : {
115 0 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
116 0 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
117 0 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
118 0 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
119 0 : constexpr Real64 r = 124.; // Pa/K , Zf/Rf, steadman1979
120 0 : constexpr Real64 eta = 1.43e-6; // kg/J , "inhaled mass" / "metabolic rate", steadman1979
121 0 : Real64 L = Le(310.); // , latent heat of vaporization at 310 K
122 0 : constexpr Real64 p = 1.013e5; // Pa , atmospheric pressure
123 0 : constexpr Real64 rgasa = 287.04; // J/kg/K
124 0 : constexpr Real64 rgasv = 461.; // J/kg/K
125 0 : constexpr Real64 cva = 719.; // J/kg/K
126 0 : constexpr Real64 cvv = 1418.; // J/kg/K
127 0 : constexpr Real64 cpa = cva + rgasa;
128 0 : constexpr Real64 cpv = cvv + rgasv;
129 :
130 0 : 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 0 : Real64 Zs(Real64 const Rs)
135 : {
136 0 : 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 0 : Real64 Ra(Real64 const Ts, Real64 const Ta)
141 : {
142 0 : constexpr Real64 sigma = 5.67e-8; // W/m^2/K^4 , Stefan-Boltzmann constant
143 0 : constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979
144 0 : constexpr Real64 hc = 17.4;
145 0 : constexpr Real64 phi_rad = 0.85;
146 0 : Real64 const hr = epsilon * phi_rad * sigma * (std::pow(Ts, 2) + std::pow(Ta, 2)) * (Ts + Ta);
147 0 : 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 0 : Real64 Ra_bar(Real64 const Tf, Real64 const Ta)
152 : {
153 0 : constexpr Real64 sigma = 5.67e-8; // W/m^2/K^4 , Stefan-Boltzmann constant
154 0 : constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979
155 0 : constexpr Real64 hc = 11.6;
156 0 : constexpr Real64 phi_rad = 0.79;
157 0 : Real64 const hr = epsilon * phi_rad * sigma * (std::pow(Tf, 2) + std::pow(Ta, 2)) * (Tf + Ta);
158 0 : return 1.0 / (hc + hr);
159 : }
160 :
161 : // Function to calculate heat transfer resistance through air, when being naked, K m^2/W
162 0 : Real64 Ra_un(Real64 const Ts, Real64 const Ta)
163 : {
164 0 : constexpr Real64 sigma = 5.67e-8; // W/m^2/K^4 , Stefan-Boltzmann constant
165 0 : constexpr Real64 epsilon = 0.97; // , emissivity of surface, steadman1979
166 0 : constexpr Real64 hc = 12.3;
167 0 : constexpr Real64 phi_rad = 0.80;
168 0 : Real64 const hr = epsilon * phi_rad * sigma * (std::pow(Ts, 2) + std::pow(Ta, 2)) * (Ts + Ta);
169 0 : return 1.0 / (hc + hr);
170 : }
171 :
172 : constexpr Real64 tol = 1e-5;
173 : constexpr Real64 maxIter = 100;
174 0 : Real64 find_eqvar_phi(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
175 : {
176 0 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
177 0 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
178 0 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
179 0 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
180 0 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
181 :
182 0 : Real64 phi = 0.84;
183 0 : Real64 const Pa = RH * pvstar(Ta);
184 0 : constexpr Real64 Rs = 0.0387;
185 0 : Real64 const ZsRs = Zs(Rs);
186 0 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
187 : Real64 Ts;
188 : int SolFla;
189 0 : General::SolveRoot(
190 : state,
191 : tol,
192 : maxIter,
193 : SolFla,
194 : Ts,
195 0 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
196 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
197 0 : std::max(Tc, Ta) + Rs * std::abs(m));
198 0 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
199 0 : if (flux1 <= 0.0) {
200 0 : phi = 1.0 - (Q - Qv(Ta, Pa)) * Rs / (Tc - Ts);
201 : }
202 0 : return phi;
203 : }
204 :
205 0 : Real64 find_eqvar_Rf(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
206 : {
207 0 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
208 0 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
209 0 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
210 0 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
211 0 : constexpr Real64 r = 124.; // Pa/K , Zf/Rf, steadman1979
212 0 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
213 0 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
214 :
215 0 : Real64 Pa = RH * pvstar(Ta);
216 0 : constexpr Real64 Rs = 0.0387;
217 0 : constexpr Real64 phi = 0.84;
218 0 : Real64 const ZsRs = Zs(Rs);
219 0 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
220 0 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
221 : Real64 Ts;
222 : int SolFla;
223 0 : General::SolveRoot(
224 : state,
225 : tol,
226 : maxIter,
227 : SolFla,
228 : Ts,
229 0 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
230 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
231 0 : std::max(Tc, Ta) + Rs * std::abs(m));
232 : Real64 Tf;
233 0 : General::SolveRoot(
234 : state,
235 : tol,
236 : maxIter,
237 : SolFla,
238 : Tf,
239 0 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
240 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
241 0 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
242 0 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
243 0 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
244 : Real64 Rf;
245 0 : if (flux1 <= 0.0) {
246 0 : Rf = std::numeric_limits<Real64>::infinity();
247 0 : } else if (flux2 <= 0.0) {
248 0 : Real64 const Ts_bar = Tc - (Q - Qv(Ta, Pa)) * Rs / phi + (1.0 / phi - 1.0) * (Tc - Ts);
249 0 : General::SolveRoot(
250 : state,
251 : tol,
252 : maxIter,
253 : SolFla,
254 : Tf,
255 0 : [&](Real64 Tf) {
256 0 : 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 0 : (Tc - Ts_bar) / Rs;
258 : },
259 : Ta,
260 : Ts_bar);
261 0 : Rf = Ra_bar(Tf, Ta) * (Ts_bar - Tf) / (Tf - Ta);
262 : } else {
263 0 : Rf = 0.0;
264 : }
265 0 : return Rf;
266 : }
267 :
268 0 : Real64 find_eqvar_rs(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
269 : {
270 0 : constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018
271 0 : constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018
272 0 : Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014
273 0 : constexpr Real64 cpc = 3492.; // J/kg/K , specific heat capacity of core, gagge1972
274 0 : Real64 C = M * cpc / A; // , heat capacity of core
275 0 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
276 0 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
277 0 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
278 0 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
279 0 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
280 0 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
281 0 : constexpr Real64 Za_un = 60.6 / 12.3; // Pa m^2/W, mass transfer resistance through air, when being naked
282 :
283 0 : Real64 Pa = RH * pvstar(Ta);
284 0 : constexpr Real64 phi = 0.84;
285 0 : Real64 Rs = 0.0387;
286 0 : Real64 ZsRs = Zs(Rs);
287 0 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
288 0 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
289 : Real64 Ts;
290 : int SolFla;
291 0 : General::SolveRoot(
292 : state,
293 : tol,
294 : maxIter,
295 : SolFla,
296 : Ts,
297 0 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
298 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
299 0 : std::max(Tc, Ta) + Rs * std::abs(m));
300 : Real64 Tf;
301 0 : General::SolveRoot(
302 : state,
303 : tol,
304 : maxIter,
305 : SolFla,
306 : Tf,
307 0 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
308 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
309 0 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
310 0 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
311 0 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
312 0 : Real64 const flux3 = Q - Qv(Ta, Pa) - (Tc - Ta) / Ra_un(Tc, Ta) - (phi_salt * pvstar(Tc) - Pa) / Za_un;
313 0 : if (flux1 > 0 && flux2 > 0) {
314 0 : if (flux3 < 0.0) {
315 0 : General::SolveRoot(
316 : state,
317 : tol,
318 : maxIter,
319 : SolFla,
320 : Ts,
321 0 : [&](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 0 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
325 0 : ZsRs = Zs(Rs);
326 0 : Real64 Ps = Pc - (Pc - Pa) * ZsRs / (ZsRs + Za_un);
327 0 : if (Ps > phi_salt * pvstar(Ts)) {
328 0 : General::SolveRoot(
329 : state,
330 : tol,
331 : maxIter,
332 : SolFla,
333 : Ts,
334 0 : [&](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 0 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
338 : }
339 : } else {
340 0 : Rs = 0.0;
341 : }
342 : }
343 0 : return Rs;
344 : }
345 :
346 0 : Real64 find_eqvar_dTcdt(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
347 : {
348 0 : constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018
349 0 : constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018
350 0 : Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014
351 0 : constexpr Real64 cpc = 3492.; // J/kg/K , specific heat capacity of core, gagge1972
352 0 : Real64 C = M * cpc / A; // , heat capacity of core
353 0 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
354 0 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
355 0 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
356 0 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
357 0 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
358 0 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
359 0 : constexpr Real64 Za_un = 60.6 / 12.3; // Pa m^2/W, mass transfer resistance through air, when being naked
360 :
361 0 : Real64 dTcdt = 0.0;
362 0 : Real64 const Pa = RH * pvstar(Ta);
363 0 : constexpr Real64 Rs = 0.0387;
364 0 : Real64 const ZsRs = Zs(Rs);
365 0 : constexpr Real64 phi = 0.84;
366 0 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
367 0 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
368 : Real64 Ts;
369 : int SolFla;
370 0 : General::SolveRoot(
371 : state,
372 : tol,
373 : maxIter,
374 : SolFla,
375 : Ts,
376 0 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
377 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
378 0 : std::max(Tc, Ta) + Rs * std::abs(m));
379 : Real64 Tf;
380 0 : General::SolveRoot(
381 : state,
382 : tol,
383 : maxIter,
384 : SolFla,
385 : Tf,
386 0 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
387 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
388 0 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
389 0 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
390 0 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
391 0 : Real64 const flux3 = Q - Qv(Ta, Pa) - (Tc - Ta) / Ra_un(Tc, Ta) - (phi_salt * pvstar(Tc) - Pa) / Za_un;
392 0 : if (flux1 > 0.0 && flux2 > 0.0 && flux3 >= 0.0) {
393 0 : dTcdt = (1.0 / C) * flux3;
394 : }
395 0 : return dTcdt;
396 : }
397 :
398 : // given T and RH, returns a key and value pair
399 0 : Real64 find_eqvar_name_and_value(EnergyPlusData &state, Real64 const Ta, Real64 const RH, EqvarName &varname)
400 : {
401 0 : constexpr Real64 M = 83.6; // kg , mass of average US adults, fryar2018
402 0 : constexpr Real64 H = 1.69; // m , height of average US adults, fryar2018
403 0 : Real64 A = 0.202 * pow(M, 0.425) * pow(H, 0.725); // m^2 , DuBois formula, parson2014
404 0 : constexpr Real64 cpc = 3492.; // J/kg/K , specific heat capacity of core, gagge1972
405 0 : Real64 C = M * cpc / A; // , heat capacity of core
406 0 : constexpr Real64 Q = 180.; // W/m^2 , metabolic rate per skin area, steadman1979
407 0 : constexpr Real64 phi_salt = 0.9; // , vapor saturation pressure level of saline solution, steadman1979
408 0 : constexpr Real64 Tc = 310.; // K , core temperature, steadman1979
409 0 : Real64 Pc = phi_salt * pvstar(Tc); // , core vapor pressure
410 0 : constexpr Real64 r = 124.; // Pa/K , Zf/Rf, steadman1979
411 0 : constexpr Real64 Za = 60.6 / 17.4; // Pa m^2/W, mass transfer resistance through air, exposed part of skin
412 0 : constexpr Real64 Za_bar = 60.6 / 11.6; // Pa m^2/W, mass transfer resistance through air, clothed part of skin
413 0 : constexpr Real64 Za_un = 60.6 / 12.3; // Pa m^2/W, mass transfer resistance through air, when being naked
414 :
415 0 : Real64 Pa = RH * pvstar(Ta);
416 0 : Real64 Rs = 0.0387;
417 0 : Real64 phi = 0.84;
418 0 : Real64 dTcdt = 0.0;
419 0 : Real64 ZsRs = Zs(Rs);
420 0 : Real64 const m = (Pc - Pa) / (ZsRs + Za);
421 0 : Real64 const m_bar = (Pc - Pa) / (ZsRs + Za_bar);
422 :
423 : int SolFla;
424 : Real64 Ts;
425 0 : General::SolveRoot(
426 : state,
427 : tol,
428 : maxIter,
429 : SolFla,
430 : Ts,
431 0 : [&](Real64 Ts) { return (Ts - Ta) / Ra(Ts, Ta) + (Pc - Pa) / (ZsRs + Za) - (Tc - Ts) / Rs; },
432 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m)),
433 0 : std::max(Tc, Ta) + Rs * std::abs(m));
434 :
435 : Real64 Tf;
436 0 : General::SolveRoot(
437 : state,
438 : tol,
439 : maxIter,
440 : SolFla,
441 : Tf,
442 0 : [&](Real64 Tf) { return (Tf - Ta) / Ra_bar(Tf, Ta) + (Pc - Pa) / (ZsRs + Za_bar) - (Tc - Tf) / Rs; },
443 0 : std::max(0.0, std::min(Tc, Ta) - Rs * std::abs(m_bar)),
444 0 : std::max(Tc, Ta) + Rs * std::abs(m_bar));
445 0 : Real64 const flux1 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs;
446 0 : Real64 const flux2 = Q - Qv(Ta, Pa) - (1.0 - phi) * (Tc - Ts) / Rs - phi * (Tc - Tf) / Rs;
447 : Real64 Rf;
448 :
449 0 : if (flux1 <= 0.0) {
450 0 : varname = EqvarName::Phi;
451 0 : phi = 1.0 - (Q - Qv(Ta, Pa)) * Rs / (Tc - Ts);
452 0 : return phi;
453 0 : } else if (flux2 <= 0.0) {
454 0 : varname = EqvarName::Rf;
455 0 : Real64 const Ts_bar = Tc - (Q - Qv(Ta, Pa)) * Rs / phi + (1.0 / phi - 1.0) * (Tc - Ts);
456 0 : General::SolveRoot(
457 : state,
458 : tol,
459 : maxIter,
460 : SolFla,
461 : Tf,
462 0 : [&](Real64 Tf) {
463 0 : 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 0 : (Tc - Ts_bar) / Rs;
465 : },
466 : Ta,
467 : Ts_bar);
468 0 : Rf = Ra_bar(Tf, Ta) * (Ts_bar - Tf) / (Tf - Ta);
469 0 : return Rf;
470 : } else {
471 0 : Real64 const flux3 = Q - Qv(Ta, Pa) - (Tc - Ta) / Ra_un(Tc, Ta) - (phi_salt * pvstar(Tc) - Pa) / Za_un;
472 0 : if (flux3 < 0.0) {
473 0 : varname = EqvarName::Rs;
474 0 : General::SolveRoot(
475 : state,
476 : tol,
477 : maxIter,
478 : SolFla,
479 : Ts,
480 0 : [&](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 0 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
484 0 : ZsRs = Zs(Rs);
485 0 : Real64 const Ps = Pc - (Pc - Pa) * ZsRs / (ZsRs + Za_un);
486 0 : if (Ps > phi_salt * pvstar(Ts)) {
487 0 : General::SolveRoot(
488 : state,
489 : tol,
490 : maxIter,
491 : SolFla,
492 : Ts,
493 0 : [&](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 0 : Rs = (Tc - Ts) / (Q - Qv(Ta, Pa));
497 : }
498 0 : return Rs;
499 : } else {
500 0 : varname = EqvarName::DTcdt;
501 0 : Rs = 0.0;
502 0 : dTcdt = (1.0 / C) * flux3;
503 0 : return dTcdt;
504 : }
505 : }
506 : }
507 :
508 : // Convert the find_T function
509 0 : Real64 find_T(EnergyPlusData &state, EqvarName const eqvar_name, Real64 const eqvar)
510 : {
511 : Real64 T;
512 : int SolFla;
513 0 : constexpr Real64 Pa0 = 1.6e3; // Pa , reference air vapor pressure in regions III, IV, V, VI, steadman1979
514 :
515 0 : if (eqvar_name == EqvarName::Phi) {
516 0 : General::SolveRoot(state, tol, maxIter, SolFla, T, [&](Real64 T) { return find_eqvar_phi(state, T, 1.0) - eqvar; }, 0.0, 240.0);
517 0 : } else if (eqvar_name == EqvarName::Rf) {
518 0 : General::SolveRoot(
519 : state,
520 : tol,
521 : maxIter,
522 : SolFla,
523 : T,
524 0 : [&](Real64 T) { return (find_eqvar_Rf(state, T, std::min(1.0, Pa0 / pvstar(T)))) - eqvar; },
525 : 230.0,
526 : 300.0);
527 0 : } else if (eqvar_name == EqvarName::Rs) {
528 0 : General::SolveRoot(
529 0 : state, tol, maxIter, SolFla, T, [&](Real64 T) { return find_eqvar_rs(state, T, Pa0 / pvstar(T)) - eqvar; }, 295.0, 350.0);
530 : } else {
531 0 : General::SolveRoot(
532 0 : state, tol, maxIter, SolFla, T, [&](Real64 T) { return find_eqvar_dTcdt(state, T, Pa0 / pvstar(T)) - eqvar; }, 340.0, 1000.0);
533 : }
534 :
535 0 : return T;
536 : }
537 :
538 0 : Real64 heatindex(EnergyPlusData &state, Real64 const Ta, Real64 const RH)
539 : {
540 :
541 : // Ta: temperature in Kelvin
542 : // RH: relative humidity in range of 0.0 to 1.0
543 : // The function computes the extended heat index, in Kelvinn
544 :
545 0 : auto const HVACSystemRootSolverMethodBackup = state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod;
546 0 : state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverAlgorithm::ShortBisectionThenRegulaFalsi;
547 0 : EqvarName eqvar_name = EqvarName::Invalid;
548 0 : Real64 const eqvar_value = find_eqvar_name_and_value(state, Ta, RH, eqvar_name);
549 :
550 0 : Real64 T = find_T(state, eqvar_name, eqvar_value);
551 :
552 0 : if (Ta == 0.0) {
553 0 : T = 0.0;
554 : }
555 :
556 0 : state.dataRootFinder->HVACSystemRootFinding.HVACSystemRootSolverMethod = HVACSystemRootSolverMethodBackup;
557 0 : return T;
558 : }
559 :
560 : } // namespace ExtendedHI
561 : } // namespace EnergyPlus
|
