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 : #include "AirflowNetwork/Elements.hpp"
49 : #include "AirflowNetwork/Properties.hpp"
50 : #include "AirflowNetwork/Solver.hpp"
51 :
52 : #include "../../Data/EnergyPlusData.hh"
53 : #include "../../DataAirLoop.hh"
54 : #include "../../DataEnvironment.hh"
55 : #include "../../DataHVACGlobals.hh"
56 : #include "../../DataLoopNode.hh"
57 : #include "../../DataSurfaces.hh"
58 :
59 : namespace EnergyPlus {
60 :
61 : namespace AirflowNetwork {
62 :
63 : // MODULE INFORMATION:
64 : // AUTHOR Lixing Gu, Don Shirey, and Muthusamy V. Swami
65 : // DATE WRITTEN Aug. 2003
66 : // MODIFIED na
67 : // RE-ENGINEERED na
68 :
69 : // PURPOSE OF THIS MODULE:
70 : // This module should contain the information that is needed to simulate
71 : // performance of air distribution system, including pressure, temperature
72 : // and moisture levels at each node, and airflow and sensible and latent energy losses
73 : // at each element
74 :
75 490145 : static Real64 square(Real64 x)
76 : {
77 490145 : return x * x;
78 : }
79 :
80 1799403 : int Duct::calculate([[maybe_unused]] EnergyPlusData &state,
81 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
82 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
83 : [[maybe_unused]] int const i, // Linkage number
84 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
85 : [[maybe_unused]] const Real64 control, // Element control signal
86 : const AirState &propN, // Node 1 properties
87 : const AirState &propM, // Node 2 properties
88 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
89 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
90 : )
91 : {
92 :
93 : // SUBROUTINE INFORMATION:
94 : // AUTHOR George Walton
95 : // DATE WRITTEN Extracted from AIRNET
96 : // MODIFIED Lixing Gu, 2/1/04
97 : // Revised the subroutine to meet E+ needs
98 : // MODIFIED Lixing Gu, 6/8/05
99 : // RE-ENGINEERED na
100 :
101 : // PURPOSE OF THIS SUBROUTINE:
102 : // This subroutine solves airflow for a duct/pipe component using Colebrook equation for the
103 : // turbulent friction factor
104 :
105 : // SUBROUTINE PARAMETER DEFINITIONS:
106 1799403 : Real64 constexpr C(0.868589);
107 1799403 : Real64 constexpr EPS(0.001);
108 :
109 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
110 : Real64 A0;
111 : Real64 A1;
112 : Real64 A2;
113 : Real64 B;
114 : Real64 D;
115 : Real64 S2;
116 : Real64 CDM;
117 : Real64 FL; // friction factor for laminar flow.
118 : Real64 FT; // friction factor for turbulent flow.
119 : Real64 FTT;
120 : Real64 RE; // Reynolds number.
121 : Real64 ed;
122 : Real64 ld;
123 : Real64 g;
124 : Real64 AA1;
125 :
126 : // Formats
127 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
128 :
129 : // CompNum = state.afn->AirflowNetworkCompData(j).TypeNum;
130 1799403 : ed = roughness / hydraulicDiameter;
131 1799403 : ld = L / hydraulicDiameter;
132 1799403 : g = 1.14 - 0.868589 * std::log(ed);
133 1799403 : AA1 = g;
134 :
135 1799403 : if (LFLAG) {
136 : // Initialization by linear relation.
137 0 : if (PDROP >= 0.0) {
138 0 : DF[0] = (2.0 * propN.density * A * hydraulicDiameter) / (propN.viscosity * InitLamCoef * ld);
139 : } else {
140 0 : DF[0] = (2.0 * propM.density * A * hydraulicDiameter) / (propM.viscosity * InitLamCoef * ld);
141 : }
142 0 : F[0] = -DF[0] * PDROP;
143 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwi:" << i << InitLamCoef << F[0] << DF[0];
144 : } else {
145 : // Standard calculation.
146 1799403 : if (PDROP >= 0.0) {
147 : // Flow in positive direction.
148 : // Laminar flow coefficient !=0
149 1702618 : if (LamFriCoef >= 0.001) {
150 939704 : A2 = LamFriCoef / (2.0 * propN.density * A * A);
151 939704 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * A * hydraulicDiameter);
152 939704 : A0 = -PDROP;
153 939704 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
154 939704 : FL = (CDM - A1) / (2.0 * A2);
155 939704 : CDM = 1.0 / CDM;
156 : } else {
157 762914 : CDM = (2.0 * propN.density * A * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
158 762914 : FL = CDM * PDROP;
159 : }
160 1702618 : RE = FL * hydraulicDiameter / (propN.viscosity * A);
161 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
162 : // Turbulent flow; test when Re>10.
163 1702618 : if (RE >= 10.0) {
164 1689985 : S2 = std::sqrt(2.0 * propN.density * PDROP) * A;
165 1689985 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
166 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
167 : while (true) {
168 3788372 : FT = FTT;
169 3788372 : B = (9.3 * propN.viscosity * A) / (FT * roughness);
170 3788372 : D = 1.0 + g * B;
171 3788372 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
172 3788372 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
173 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
174 3788372 : if (std::abs(FTT - FT) / FTT < EPS) break;
175 : }
176 1689985 : FT = FTT;
177 : } else {
178 12633 : FT = FL;
179 : }
180 : } else {
181 : // Flow in negative direction.
182 : // Laminar flow coefficient !=0
183 96785 : if (LamFriCoef >= 0.001) {
184 43450 : A2 = LamFriCoef / (2.0 * propM.density * A * A);
185 43450 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * A * hydraulicDiameter);
186 43450 : A0 = PDROP;
187 43450 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
188 43450 : FL = -(CDM - A1) / (2.0 * A2);
189 43450 : CDM = 1.0 / CDM;
190 : } else {
191 53335 : CDM = (2.0 * propM.density * A * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
192 53335 : FL = CDM * PDROP;
193 : }
194 96785 : RE = -FL * hydraulicDiameter / (propM.viscosity * A);
195 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
196 : // Turbulent flow; test when Re>10.
197 96785 : if (RE >= 10.0) {
198 91387 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * A;
199 91387 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
200 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
201 : while (true) {
202 319448 : FT = FTT;
203 319448 : B = (9.3 * propM.viscosity * A) / (FT * roughness);
204 319448 : D = 1.0 + g * B;
205 319448 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
206 319448 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
207 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
208 319448 : if (std::abs(FTT - FT) / FTT < EPS) break;
209 : }
210 91387 : FT = -FTT;
211 : } else {
212 5398 : FT = FL;
213 : }
214 : }
215 : // Select laminar or turbulent flow.
216 1799403 : if (std::abs(FL) <= std::abs(FT)) {
217 64777 : F[0] = FL;
218 64777 : DF[0] = CDM;
219 : } else {
220 1734626 : F[0] = FT;
221 1734626 : DF[0] = 0.5 * FT / PDROP;
222 : }
223 : }
224 1799403 : return 1;
225 : }
226 :
227 0 : int Duct::calculate([[maybe_unused]] EnergyPlusData &state,
228 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
229 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
230 : [[maybe_unused]] const Real64 control, // Element control signal
231 : const AirState &propN, // Node 1 properties
232 : const AirState &propM, // Node 2 properties
233 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
234 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
235 : )
236 : {
237 :
238 : // SUBROUTINE INFORMATION:
239 : // AUTHOR George Walton
240 : // DATE WRITTEN Extracted from AIRNET
241 : // MODIFIED Lixing Gu, 2/1/04
242 : // Revised the subroutine to meet E+ needs
243 : // MODIFIED Lixing Gu, 6/8/05
244 : // RE-ENGINEERED na
245 :
246 : // PURPOSE OF THIS SUBROUTINE:
247 : // This subroutine solves airflow for a duct/pipe component using Colebrook equation for the
248 : // turbulent friction factor
249 :
250 : // SUBROUTINE PARAMETER DEFINITIONS:
251 0 : Real64 constexpr C(0.868589);
252 0 : Real64 constexpr EPS(0.001);
253 :
254 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
255 : Real64 A0;
256 : Real64 A1;
257 : Real64 A2;
258 : Real64 B;
259 : Real64 D;
260 : Real64 S2;
261 : Real64 CDM;
262 : Real64 FL; // friction factor for laminar flow.
263 : Real64 FT; // friction factor for turbulent flow.
264 : Real64 FTT;
265 : Real64 RE; // Reynolds number.
266 : Real64 ed;
267 : Real64 ld;
268 : Real64 g;
269 : Real64 AA1;
270 :
271 : // Formats
272 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
273 :
274 : // CompNum = state.afn->AirflowNetworkCompData(j).TypeNum;
275 0 : ed = roughness / hydraulicDiameter;
276 0 : ld = L / hydraulicDiameter;
277 0 : g = 1.14 - 0.868589 * std::log(ed);
278 0 : AA1 = g;
279 :
280 : // Standard calculation.
281 0 : if (PDROP >= 0.0) {
282 : // Flow in positive direction.
283 : // Laminar flow coefficient !=0
284 0 : if (LamFriCoef >= 0.001) {
285 0 : A2 = LamFriCoef / (2.0 * propN.density * A * A);
286 0 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * A * hydraulicDiameter);
287 0 : A0 = -PDROP;
288 0 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
289 0 : FL = (CDM - A1) / (2.0 * A2);
290 0 : CDM = 1.0 / CDM;
291 : } else {
292 0 : CDM = (2.0 * propN.density * A * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
293 0 : FL = CDM * PDROP;
294 : }
295 0 : RE = FL * hydraulicDiameter / (propN.viscosity * A);
296 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
297 : // Turbulent flow; test when Re>10.
298 0 : if (RE >= 10.0) {
299 0 : S2 = std::sqrt(2.0 * propN.density * PDROP) * A;
300 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
301 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
302 : while (true) {
303 0 : FT = FTT;
304 0 : B = (9.3 * propN.viscosity * A) / (FT * roughness);
305 0 : D = 1.0 + g * B;
306 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
307 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
308 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
309 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
310 : }
311 0 : FT = FTT;
312 : } else {
313 0 : FT = FL;
314 : }
315 : } else {
316 : // Flow in negative direction.
317 : // Laminar flow coefficient !=0
318 0 : if (LamFriCoef >= 0.001) {
319 0 : A2 = LamFriCoef / (2.0 * propM.density * A * A);
320 0 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * A * hydraulicDiameter);
321 0 : A0 = PDROP;
322 0 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
323 0 : FL = -(CDM - A1) / (2.0 * A2);
324 0 : CDM = 1.0 / CDM;
325 : } else {
326 0 : CDM = (2.0 * propM.density * A * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
327 0 : FL = CDM * PDROP;
328 : }
329 0 : RE = -FL * hydraulicDiameter / (propM.viscosity * A);
330 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
331 : // Turbulent flow; test when Re>10.
332 0 : if (RE >= 10.0) {
333 0 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * A;
334 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
335 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
336 : while (true) {
337 0 : FT = FTT;
338 0 : B = (9.3 * propM.viscosity * A) / (FT * roughness);
339 0 : D = 1.0 + g * B;
340 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
341 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
342 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
343 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
344 : }
345 0 : FT = -FTT;
346 : } else {
347 0 : FT = FL;
348 : }
349 : }
350 : // Select laminar or turbulent flow.
351 0 : if (std::abs(FL) <= std::abs(FT)) {
352 0 : F[0] = FL;
353 0 : DF[0] = CDM;
354 : } else {
355 0 : F[0] = FT;
356 0 : DF[0] = 0.5 * FT / PDROP;
357 : }
358 0 : return 1;
359 : }
360 :
361 2583 : int SurfaceCrack::calculate([[maybe_unused]] EnergyPlusData &state,
362 : bool const linear, // Initialization flag. If true, use linear relationship
363 : Real64 const pdrop, // Total pressure drop across a component (P1 - P2) [Pa]
364 : [[maybe_unused]] int const i, // Linkage number
365 : const Real64 multiplier, // Element multiplier
366 : const Real64 control, // Element control signal
367 : const AirState &propN, // Node 1 properties
368 : const AirState &propM, // Node 2 properties
369 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
370 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
371 : )
372 : {
373 : // SUBROUTINE INFORMATION:
374 : // AUTHOR George Walton
375 : // DATE WRITTEN Extracted from AIRNET
376 : // MODIFIED Lixing Gu, 2/1/04
377 : // Revised the subroutine to meet E+ needs
378 : // MODIFIED Lixing Gu, 6/8/05
379 : // RE-ENGINEERED Jason DeGraw
380 :
381 : // PURPOSE OF THIS SUBROUTINE:
382 : // This subroutine solves airflow for a surface crack component
383 :
384 : // Real64 rhoz_norm = AIRDENSITY(StandardP, StandardT, StandardW);
385 : // Real64 viscz_norm = 1.71432e-5 + 4.828e-8 * StandardT;
386 :
387 2583 : Real64 VisAve{0.5 * (propN.viscosity + propM.viscosity)};
388 2583 : Real64 Tave{0.5 * (propN.temperature + propM.temperature)};
389 :
390 2583 : Real64 sign{1.0};
391 2583 : Real64 upwind_temperature{propN.temperature};
392 2583 : Real64 upwind_density{propN.density};
393 2583 : Real64 upwind_viscosity{propN.viscosity};
394 2583 : Real64 upwind_sqrt_density{propN.sqrt_density};
395 2583 : Real64 abs_pdrop = pdrop;
396 :
397 2583 : if (pdrop < 0.0) {
398 614 : sign = -1.0;
399 614 : upwind_temperature = propM.temperature;
400 614 : upwind_density = propM.density;
401 614 : upwind_viscosity = propM.viscosity;
402 614 : upwind_sqrt_density = propM.sqrt_density;
403 614 : abs_pdrop = -pdrop;
404 : }
405 :
406 2583 : Real64 coef = coefficient * control * multiplier / upwind_sqrt_density;
407 :
408 : // Laminar calculation
409 2583 : Real64 RhoCor{TOKELVIN(upwind_temperature) / TOKELVIN(Tave)};
410 2583 : Real64 Ctl{std::pow(reference_density / upwind_density / RhoCor, exponent - 1.0) *
411 2583 : std::pow(reference_viscosity / VisAve, 2.0 * exponent - 1.0)};
412 2583 : Real64 CDM{coef * upwind_density / upwind_viscosity * Ctl};
413 2583 : Real64 FL{CDM * pdrop};
414 : Real64 abs_FT;
415 :
416 2583 : if (linear) {
417 2 : DF[0] = CDM;
418 2 : F[0] = FL;
419 : } else {
420 : // Turbulent flow.
421 2581 : if (exponent == 0.5) {
422 0 : abs_FT = coef * upwind_sqrt_density * std::sqrt(abs_pdrop) * Ctl;
423 : } else {
424 2581 : abs_FT = coef * upwind_sqrt_density * std::pow(abs_pdrop, exponent) * Ctl;
425 : }
426 : // Select linear or turbulent flow.
427 2581 : if (std::abs(FL) <= abs_FT) {
428 102 : F[0] = FL;
429 102 : DF[0] = CDM;
430 : } else {
431 2479 : F[0] = sign * abs_FT;
432 2479 : DF[0] = F[0] * exponent / pdrop;
433 : }
434 : }
435 2583 : return 1;
436 : }
437 :
438 0 : int SurfaceCrack::calculate(EnergyPlusData &state,
439 : Real64 const pdrop, // Total pressure drop across a component (P1 - P2) [Pa]
440 : const Real64 multiplier, // Element multiplier
441 : const Real64 control, // Element control signal
442 : const AirState &propN, // Node 1 properties
443 : const AirState &propM, // Node 2 properties
444 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
445 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
446 : )
447 : {
448 : // SUBROUTINE INFORMATION:
449 : // AUTHOR George Walton
450 : // DATE WRITTEN Extracted from AIRNET
451 : // MODIFIED Lixing Gu, 2/1/04
452 : // Revised the subroutine to meet E+ needs
453 : // MODIFIED Lixing Gu, 6/8/05
454 : // RE-ENGINEERED Jason DeGraw
455 :
456 : // PURPOSE OF THIS SUBROUTINE:
457 : // This subroutine solves airflow for a surface crack component
458 :
459 : // Real64 rhoz_norm = AIRDENSITY(StandardP, StandardT, StandardW);
460 : // Real64 viscz_norm = 1.71432e-5 + 4.828e-8 * StandardT;
461 :
462 0 : Real64 VisAve{0.5 * (propN.viscosity + propM.viscosity)};
463 0 : Real64 Tave{0.5 * (propN.temperature + propM.temperature)};
464 :
465 0 : Real64 sign{1.0};
466 0 : Real64 upwind_temperature{propN.temperature};
467 0 : Real64 upwind_density{propN.density};
468 0 : Real64 upwind_viscosity{propN.viscosity};
469 0 : Real64 upwind_sqrt_density{propN.sqrt_density};
470 0 : Real64 abs_pdrop = pdrop;
471 :
472 0 : if (pdrop < 0.0) {
473 0 : sign = -1.0;
474 0 : upwind_temperature = propM.temperature;
475 0 : upwind_density = propM.density;
476 0 : upwind_viscosity = propM.viscosity;
477 0 : upwind_sqrt_density = propM.sqrt_density;
478 0 : abs_pdrop = -pdrop;
479 : }
480 :
481 0 : Real64 coef = coefficient * control * multiplier / upwind_sqrt_density;
482 :
483 : // Laminar calculation
484 0 : Real64 RhoCor{TOKELVIN(upwind_temperature) / TOKELVIN(Tave)};
485 0 : Real64 Ctl{std::pow(reference_density / upwind_density / RhoCor, exponent - 1.0) *
486 0 : std::pow(reference_viscosity / VisAve, 2.0 * exponent - 1.0)};
487 0 : Real64 CDM{coef * upwind_density / upwind_viscosity * Ctl};
488 0 : Real64 FL{CDM * pdrop};
489 : Real64 abs_FT;
490 :
491 : // Turbulent flow.
492 0 : if (exponent == 0.5) {
493 0 : abs_FT = coef * upwind_sqrt_density * std::sqrt(abs_pdrop) * Ctl;
494 : } else {
495 0 : abs_FT = coef * upwind_sqrt_density * std::pow(abs_pdrop, exponent) * Ctl;
496 : }
497 : // Select linear or turbulent flow.
498 0 : if (std::abs(FL) <= abs_FT) {
499 0 : F[0] = FL;
500 0 : DF[0] = CDM;
501 : } else {
502 0 : F[0] = sign * abs_FT;
503 0 : DF[0] = F[0] * exponent / pdrop;
504 : }
505 :
506 0 : return 1;
507 : }
508 :
509 172 : int DuctLeak::calculate(EnergyPlusData &state,
510 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
511 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
512 : [[maybe_unused]] int const i, // Linkage number
513 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
514 : [[maybe_unused]] const Real64 control, // Element control signal
515 : const AirState &propN, // Node 1 properties
516 : const AirState &propM, // Node 2 properties
517 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
518 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
519 : )
520 : {
521 : // SUBROUTINE INFORMATION:
522 : // AUTHOR George Walton
523 : // DATE WRITTEN Extracted from AIRNET
524 : // MODIFIED Lixing Gu, 2/1/04
525 : // Revised the subroutine to meet E+ needs
526 : // MODIFIED Lixing Gu, 6/8/05
527 : // RE-ENGINEERED na
528 :
529 : // PURPOSE OF THIS SUBROUTINE:
530 : // This subroutine solves airflow for a power law resistance airflow component
531 :
532 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
533 : Real64 CDM;
534 : Real64 FL;
535 : Real64 FT;
536 : Real64 Ctl;
537 :
538 : // Formats
539 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
540 :
541 : // Crack standard condition: T=20C, p=101325 Pa and 0 g/kg
542 172 : Real64 RhozNorm = state.afn->properties.density(101325.0, 20.0, 0.0);
543 172 : Real64 VisczNorm = 1.71432e-5 + 4.828e-8 * 20.0;
544 172 : Real64 coef = FlowCoef;
545 :
546 172 : if (PDROP >= 0.0) {
547 158 : coef /= propN.sqrt_density;
548 : } else {
549 14 : coef /= propM.sqrt_density;
550 : }
551 :
552 172 : if (LFLAG) {
553 : // Initialization by linear relation.
554 0 : if (PDROP >= 0.0) {
555 0 : Ctl = std::pow(RhozNorm / propN.density, FlowExpo - 1.0) * std::pow(VisczNorm / propN.viscosity, 2.0 * FlowExpo - 1.0);
556 0 : DF[0] = coef * propN.density / propN.viscosity * Ctl;
557 : } else {
558 0 : Ctl = std::pow(RhozNorm / propM.density, FlowExpo - 1.0) * std::pow(VisczNorm / propM.viscosity, 2.0 * FlowExpo - 1.0);
559 0 : DF[0] = coef * propM.density / propM.viscosity * Ctl;
560 : }
561 0 : F[0] = -DF[0] * PDROP;
562 : } else {
563 : // Standard calculation.
564 172 : if (PDROP >= 0.0) {
565 : // Flow in positive direction for laminar flow.
566 158 : Ctl = std::pow(RhozNorm / propN.density, FlowExpo - 1.0) * std::pow(VisczNorm / propN.viscosity, 2.0 * FlowExpo - 1.0);
567 158 : CDM = coef * propN.density / propN.viscosity * Ctl;
568 158 : FL = CDM * PDROP;
569 : // Flow in positive direction for turbulent flow.
570 158 : if (FlowExpo == 0.5) {
571 0 : FT = coef * propN.sqrt_density * std::sqrt(PDROP);
572 : } else {
573 158 : FT = coef * propN.sqrt_density * std::pow(PDROP, FlowExpo);
574 : }
575 : } else {
576 : // Flow in negative direction for laminar flow
577 14 : CDM = coef * propM.density / propM.viscosity;
578 14 : FL = CDM * PDROP;
579 : // Flow in negative direction for turbulent flow
580 14 : if (FlowExpo == 0.5) {
581 0 : FT = -coef * propM.sqrt_density * std::sqrt(-PDROP);
582 : } else {
583 14 : FT = -coef * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
584 : }
585 : }
586 : // Select laminar or turbulent flow.
587 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " plr: " << i << PDROP << FL << FT;
588 172 : if (std::abs(FL) <= std::abs(FT)) {
589 12 : F[0] = FL;
590 12 : DF[0] = CDM;
591 : } else {
592 160 : F[0] = FT;
593 160 : DF[0] = FT * FlowExpo / PDROP;
594 : }
595 : }
596 172 : return 1;
597 : }
598 :
599 0 : int DuctLeak::calculate(EnergyPlusData &state,
600 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
601 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
602 : [[maybe_unused]] const Real64 control, // Element control signal
603 : const AirState &propN, // Node 1 properties
604 : const AirState &propM, // Node 2 properties
605 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
606 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
607 : )
608 : {
609 : // SUBROUTINE INFORMATION:
610 : // AUTHOR George Walton
611 : // DATE WRITTEN Extracted from AIRNET
612 : // MODIFIED Lixing Gu, 2/1/04
613 : // Revised the subroutine to meet E+ needs
614 : // MODIFIED Lixing Gu, 6/8/05
615 : // RE-ENGINEERED na
616 :
617 : // PURPOSE OF THIS SUBROUTINE:
618 : // This subroutine solves airflow for a power law resistance airflow component
619 :
620 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
621 : Real64 CDM;
622 : Real64 FL;
623 : Real64 FT;
624 : Real64 Ctl;
625 :
626 : // Formats
627 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
628 :
629 : // Crack standard condition: T=20C, p=101325 Pa and 0 g/kg
630 0 : Real64 RhozNorm = state.afn->properties.density(101325.0, 20.0, 0.0);
631 0 : Real64 VisczNorm = 1.71432e-5 + 4.828e-8 * 20.0;
632 0 : Real64 coef = FlowCoef;
633 :
634 0 : if (PDROP >= 0.0) {
635 0 : coef /= propN.sqrt_density;
636 : } else {
637 0 : coef /= propM.sqrt_density;
638 : }
639 :
640 : // Standard calculation.
641 0 : if (PDROP >= 0.0) {
642 : // Flow in positive direction for laminar flow.
643 0 : Ctl = std::pow(RhozNorm / propN.density, FlowExpo - 1.0) * std::pow(VisczNorm / propN.viscosity, 2.0 * FlowExpo - 1.0);
644 0 : CDM = coef * propN.density / propN.viscosity * Ctl;
645 0 : FL = CDM * PDROP;
646 : // Flow in positive direction for turbulent flow.
647 0 : if (FlowExpo == 0.5) {
648 0 : FT = coef * propN.sqrt_density * std::sqrt(PDROP);
649 : } else {
650 0 : FT = coef * propN.sqrt_density * std::pow(PDROP, FlowExpo);
651 : }
652 : } else {
653 : // Flow in negative direction for laminar flow
654 0 : CDM = coef * propM.density / propM.viscosity;
655 0 : FL = CDM * PDROP;
656 : // Flow in negative direction for turbulent flow
657 0 : if (FlowExpo == 0.5) {
658 0 : FT = -coef * propM.sqrt_density * std::sqrt(-PDROP);
659 : } else {
660 0 : FT = -coef * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
661 : }
662 : }
663 : // Select laminar or turbulent flow.
664 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " plr: " << i << PDROP << FL << FT;
665 0 : if (std::abs(FL) <= std::abs(FT)) {
666 0 : F[0] = FL;
667 0 : DF[0] = CDM;
668 : } else {
669 0 : F[0] = FT;
670 0 : DF[0] = FT * FlowExpo / PDROP;
671 : }
672 0 : return 1;
673 : }
674 :
675 163425 : int ConstantVolumeFan::calculate(EnergyPlusData &state,
676 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
677 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
678 : int const i, // Linkage number
679 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
680 : [[maybe_unused]] const Real64 control, // Element control signal
681 : const AirState &propN, // Node 1 properties
682 : const AirState &propM, // Node 2 properties
683 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
684 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
685 : )
686 : {
687 :
688 : // SUBROUTINE INFORMATION:
689 : // AUTHOR George Walton
690 : // DATE WRITTEN Extracted from AIRNET
691 : // MODIFIED Lixing Gu, 2/1/04
692 : // Revised the subroutine to meet E+ needs
693 : // MODIFIED Lixing Gu, 6/8/05
694 : // RE-ENGINEERED na
695 :
696 : // PURPOSE OF THIS SUBROUTINE:
697 : // This subroutine solves airflow for a constant flow rate airflow component -- using standard interface.
698 :
699 : // Using/Aliasing
700 163425 : auto &NumPrimaryAirSys = state.dataHVACGlobal->NumPrimaryAirSys;
701 :
702 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
703 : int k;
704 : int k1;
705 : Real64 SumTermFlow; // Sum of all Terminal flows [kg/s]
706 : Real64 SumFracSuppLeak; // Sum of all supply leaks as a fraction of supply fan flow rate
707 : int Node1;
708 : int Node2;
709 :
710 163425 : int NF(1);
711 :
712 163425 : int AirLoopNum = state.afn->AirflowNetworkLinkageData(i).AirLoopNum;
713 :
714 163425 : if (fanType == HVAC::FanType::OnOff) {
715 326586 : if (state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopFanOperationMode == HVAC::FanOp::Cycling &&
716 163293 : state.dataLoopNodes->Node(InletNode).MassFlowRate == 0.0) {
717 0 : NF = GenericDuct(0.1, 0.001, LFLAG, PDROP, propN, propM, F, DF);
718 326586 : } else if (state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopFanOperationMode == HVAC::FanOp::Cycling &&
719 163293 : state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopSystemOnMassFlowrate > 0.0) {
720 163293 : F[0] = state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopSystemOnMassFlowrate;
721 : } else {
722 0 : F[0] = state.dataLoopNodes->Node(InletNode).MassFlowRate * Ctrl;
723 0 : if (state.afn->MultiSpeedHPIndicator == 2) {
724 0 : F[0] = state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopSystemOnMassFlowrate *
725 0 : state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopCompCycRatio +
726 0 : state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopSystemOffMassFlowrate *
727 0 : (1.0 - state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopCompCycRatio);
728 : }
729 : }
730 132 : } else if (fanType == HVAC::FanType::Constant) {
731 132 : if (state.dataLoopNodes->Node(InletNode).MassFlowRate > 0.0) {
732 132 : F[0] = FlowRate * Ctrl;
733 0 : } else if (NumPrimaryAirSys > 1 && state.dataLoopNodes->Node(InletNode).MassFlowRate <= 0.0) {
734 0 : NF = GenericDuct(0.1, 0.001, LFLAG, PDROP, propN, propM, F, DF);
735 : }
736 :
737 132 : if (state.afn->MultiSpeedHPIndicator == 2) {
738 0 : F[0] = state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopSystemOnMassFlowrate;
739 : }
740 0 : } else if (fanType == HVAC::FanType::VAV) {
741 : // Check VAV termals with a damper
742 0 : SumTermFlow = 0.0;
743 0 : SumFracSuppLeak = 0.0;
744 0 : for (k = 1; k <= state.afn->ActualNumOfLinks; ++k) {
745 0 : if (state.afn->AirflowNetworkLinkageData(k).VAVTermDamper && state.afn->AirflowNetworkLinkageData(k).AirLoopNum == AirLoopNum) {
746 0 : k1 = state.afn->AirflowNetworkNodeData(state.afn->AirflowNetworkLinkageData(k).NodeNums[0]).EPlusNodeNum;
747 0 : if (state.dataLoopNodes->Node(k1).MassFlowRate > 0.0) {
748 0 : SumTermFlow += state.dataLoopNodes->Node(k1).MassFlowRate;
749 : }
750 : }
751 0 : if (state.afn->AirflowNetworkCompData(state.afn->AirflowNetworkLinkageData(k).CompNum).CompTypeNum == iComponentTypeNum::ELR) {
752 : // Calculate supply leak sensible losses
753 0 : Node1 = state.afn->AirflowNetworkLinkageData(k).NodeNums[0];
754 0 : Node2 = state.afn->AirflowNetworkLinkageData(k).NodeNums[1];
755 0 : if ((state.afn->AirflowNetworkNodeData(Node2).EPlusZoneNum > 0) && (state.afn->AirflowNetworkNodeData(Node1).EPlusNodeNum == 0) &&
756 0 : (state.afn->AirflowNetworkNodeData(Node1).AirLoopNum == AirLoopNum)) {
757 0 : SumFracSuppLeak +=
758 0 : state.afn->DisSysCompELRData(state.afn->AirflowNetworkCompData(state.afn->AirflowNetworkLinkageData(k).CompNum).TypeNum)
759 0 : .ELR;
760 : }
761 : }
762 : }
763 0 : F[0] = SumTermFlow / (1.0 - SumFracSuppLeak);
764 0 : state.afn->VAVTerminalRatio = 0.0;
765 0 : if (F[0] > MaxAirMassFlowRate) {
766 0 : state.afn->VAVTerminalRatio = MaxAirMassFlowRate / F[0];
767 0 : F[0] = MaxAirMassFlowRate;
768 : }
769 : }
770 163425 : DF[0] = 0.0;
771 163425 : return NF;
772 : }
773 :
774 0 : int DetailedFan::calculate(EnergyPlusData &state,
775 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
776 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
777 : int const i, // Linkage number
778 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
779 : const Real64 control, // Element control signal
780 : const AirState &propN, // Node 1 properties
781 : const AirState &propM, // Node 2 properties
782 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
783 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
784 : )
785 : {
786 :
787 : // SUBROUTINE INFORMATION:
788 : // AUTHOR George Walton
789 : // DATE WRITTEN Extracted from AIRNET
790 : // MODIFIED Lixing Gu, 2/1/04
791 : // Revised the subroutine to meet E+ needs
792 : // MODIFIED Lixing Gu, 6/8/05
793 : // RE-ENGINEERED na
794 :
795 : // PURPOSE OF THIS SUBROUTINE:
796 : // This subroutine solves airflow for a detailed fan component -- using standard interface.
797 :
798 : // SUBROUTINE PARAMETER DEFINITIONS:
799 0 : Real64 constexpr TOL(0.00001);
800 :
801 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
802 : int j;
803 : int k;
804 : int L;
805 : Real64 DPDF;
806 : Real64 PRISE; // pressure rise (negative of pressure drop) (Pa).
807 : Real64 BX;
808 : Real64 BY;
809 : Real64 CX;
810 : Real64 CY;
811 : Real64 CCY;
812 : Real64 DX;
813 : Real64 DY;
814 :
815 : // Formats
816 : // static gio::Fmt Format_901("(A5,I3,5E14.6)");
817 :
818 0 : int NumCur = n;
819 0 : if (control <= 0.0) {
820 : // Speed = 0; treat fan as resistance.
821 0 : generic_crack(FlowCoef, FlowExpo, LFLAG, PDROP, propN, propM, F, DF);
822 0 : return 1;
823 : }
824 : // Pressure rise at reference fan speed.
825 0 : if (control >= TranRat) {
826 0 : PRISE = -PDROP * (RhoAir / propN.density) / (control * control);
827 : } else {
828 0 : PRISE = -PDROP * (RhoAir / propN.density) / (TranRat * control);
829 : }
830 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " fan:" << i << PDROP << PRISE << AFECTL(i) << DisSysCompDetFanData(CompNum).TranRat;
831 0 : if (LFLAG) {
832 : // Initialization by linear approximation.
833 0 : F[0] = -Qfree * control * (1.0 - PRISE / Pshut);
834 0 : DPDF = -Pshut / Qfree;
835 : // if (LIST >= 4)
836 : // gio::write(Unit21, Format_901) << " fni:" << JA << DisSysCompDetFanData(CompNum).Qfree << DisSysCompDetFanData(CompNum).Pshut;
837 : } else {
838 : // Solution of the fan performance curve.
839 : // Determine curve fit range.
840 0 : j = 1;
841 0 : k = 5 * (j - 1) + 1;
842 0 : BX = Coeff(k);
843 0 : BY = Coeff(k + 1) + BX * (Coeff(k + 2) + BX * (Coeff(k + 3) + BX * Coeff(k + 4))) - PRISE;
844 0 : if (BY < 0.0) ShowFatalError(state, "Out of range, too low in an AirflowNetwork detailed Fan");
845 :
846 : while (true) {
847 0 : DX = Coeff(k + 5);
848 0 : DY = Coeff(k + 1) + DX * (Coeff(k + 2) + DX * (Coeff(k + 3) + DX * Coeff(k + 5))) - PRISE;
849 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " fp0:" << j << BX << BY << DX << DY;
850 0 : if (BY * DY <= 0.0) break;
851 0 : ++j;
852 0 : if (j > NumCur) ShowFatalError(state, "Out of range, too high (FAN) in ADS simulation");
853 0 : k += 5;
854 0 : BX = DX;
855 0 : BY = DY;
856 : }
857 : // Determine reference mass flow rate by false position method.
858 0 : L = 0;
859 0 : CY = 0.0;
860 0 : Label40:;
861 0 : ++L;
862 0 : if (L > 100) ShowFatalError(state, "Too many iterations (FAN) in AirflowNtework simulation");
863 0 : CCY = CY;
864 0 : CX = BX - BY * ((DX - BX) / (DY - BY));
865 0 : CY = Coeff(k + 1) + CX * (Coeff(k + 2) + CX * (Coeff(k + 3) + CX * Coeff(k + 4))) - PRISE;
866 0 : if (BY * CY == 0.0) goto Label90;
867 0 : if (BY * CY > 0.0) goto Label60;
868 0 : DX = CX;
869 0 : DY = CY;
870 0 : if (CY * CCY > 0.0) BY *= 0.5;
871 0 : goto Label70;
872 0 : Label60:;
873 0 : BX = CX;
874 0 : BY = CY;
875 0 : if (CY * CCY > 0.0) DY *= 0.5;
876 0 : Label70:;
877 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " fpi:" << j << BX << CX << DX << BY << DY;
878 0 : if (DX - BX < TOL * CX) goto Label80;
879 0 : if (DX - BX < TOL) goto Label80;
880 0 : goto Label40;
881 0 : Label80:;
882 0 : CX = 0.5 * (BX + DX);
883 0 : Label90:;
884 0 : F[0] = CX;
885 0 : DPDF = Coeff(k + 2) + CX * (2.0 * Coeff(k + 3) + CX * 3.0 * Coeff(k + 4));
886 : }
887 : // Convert to flow at given speed.
888 0 : F[0] *= (propN.density / RhoAir) * control;
889 : // Set derivative w/r pressure drop (-).
890 0 : if (control >= TranRat) {
891 0 : DF[0] = -control / DPDF;
892 : } else {
893 0 : DF[0] = -1.0 / DPDF;
894 : }
895 0 : return 1;
896 : }
897 :
898 0 : int DetailedFan::calculate(EnergyPlusData &state,
899 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
900 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
901 : const Real64 control, // Element control signal
902 : const AirState &propN, // Node 1 properties
903 : const AirState &propM, // Node 2 properties
904 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
905 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
906 : )
907 : {
908 :
909 : // SUBROUTINE INFORMATION:
910 : // AUTHOR George Walton
911 : // DATE WRITTEN Extracted from AIRNET
912 : // MODIFIED Lixing Gu, 2/1/04
913 : // Revised the subroutine to meet E+ needs
914 : // MODIFIED Lixing Gu, 6/8/05
915 : // RE-ENGINEERED na
916 :
917 : // PURPOSE OF THIS SUBROUTINE:
918 : // This subroutine solves airflow for a detailed fan component -- using standard interface.
919 :
920 : // SUBROUTINE PARAMETER DEFINITIONS:
921 0 : Real64 constexpr TOL(0.00001);
922 :
923 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
924 : int j;
925 : int k;
926 : int L;
927 : Real64 DPDF;
928 : Real64 PRISE; // pressure rise (negative of pressure drop) (Pa).
929 : Real64 BX;
930 : Real64 BY;
931 : Real64 CX;
932 : Real64 CY;
933 : Real64 CCY;
934 : Real64 DX;
935 : Real64 DY;
936 :
937 : // Formats
938 : // static gio::Fmt Format_901("(A5,I3,5E14.6)");
939 :
940 0 : int NumCur = n;
941 :
942 0 : if (control <= 0.0) {
943 : // Speed = 0; treat fan as resistance.
944 0 : generic_crack(FlowCoef, FlowExpo, false, PDROP, propN, propM, F, DF);
945 0 : return 1;
946 : }
947 : // Pressure rise at reference fan speed.
948 0 : if (control >= TranRat) {
949 0 : PRISE = -PDROP * (RhoAir / propN.density) / pow_2(control);
950 : } else {
951 0 : PRISE = -PDROP * (RhoAir / propN.density) / (TranRat * control);
952 : }
953 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " fan:" << i << PDROP << PRISE << AFECTL(i) << DisSysCompDetFanData(CompNum).TranRat;
954 : // if (LFLAG) {
955 : // // Initialization by linear approximation.
956 : // F[0] = -Qfree * control * (1.0 - PRISE / Pshut);
957 : // DPDF = -Pshut / Qfree;
958 : // // if (LIST >= 4)
959 : // // gio::write(Unit21, Format_901) << " fni:" << JA << DisSysCompDetFanData(CompNum).Qfree << DisSysCompDetFanData(CompNum).Pshut;
960 : //} else {
961 : // Solution of the fan performance curve.
962 : // Determine curve fit range.
963 0 : j = 1;
964 0 : k = 5 * (j - 1) + 1;
965 0 : BX = Coeff(k);
966 0 : BY = Coeff(k + 1) + BX * (Coeff(k + 2) + BX * (Coeff(k + 3) + BX * Coeff(k + 4))) - PRISE;
967 0 : if (BY < 0.0) ShowFatalError(state, "Out of range, too low in an AirflowNetwork detailed Fan");
968 :
969 : while (true) {
970 0 : DX = Coeff(k + 5);
971 0 : DY = Coeff(k + 1) + DX * (Coeff(k + 2) + DX * (Coeff(k + 3) + DX * Coeff(k + 5))) - PRISE;
972 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " fp0:" << j << BX << BY << DX << DY;
973 0 : if (BY * DY <= 0.0) break;
974 0 : ++j;
975 0 : if (j > NumCur) ShowFatalError(state, "Out of range, too high (FAN) in ADS simulation");
976 0 : k += 5;
977 0 : BX = DX;
978 0 : BY = DY;
979 : }
980 : // Determine reference mass flow rate by false position method.
981 0 : L = 0;
982 0 : CY = 0.0;
983 0 : Label40:;
984 0 : ++L;
985 0 : if (L > 100) ShowFatalError(state, "Too many iterations (FAN) in AirflowNtework simulation");
986 0 : CCY = CY;
987 0 : CX = BX - BY * ((DX - BX) / (DY - BY));
988 0 : CY = Coeff(k + 1) + CX * (Coeff(k + 2) + CX * (Coeff(k + 3) + CX * Coeff(k + 4))) - PRISE;
989 0 : if (BY * CY == 0.0) goto Label90;
990 0 : if (BY * CY > 0.0) goto Label60;
991 0 : DX = CX;
992 0 : DY = CY;
993 0 : if (CY * CCY > 0.0) BY *= 0.5;
994 0 : goto Label70;
995 0 : Label60:;
996 0 : BX = CX;
997 0 : BY = CY;
998 0 : if (CY * CCY > 0.0) DY *= 0.5;
999 0 : Label70:;
1000 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " fpi:" << j << BX << CX << DX << BY << DY;
1001 0 : if (DX - BX < TOL * CX) goto Label80;
1002 0 : if (DX - BX < TOL) goto Label80;
1003 0 : goto Label40;
1004 0 : Label80:;
1005 0 : CX = 0.5 * (BX + DX);
1006 0 : Label90:;
1007 0 : F[0] = CX;
1008 0 : DPDF = Coeff(k + 2) + CX * (2.0 * Coeff(k + 3) + CX * 3.0 * Coeff(k + 4));
1009 :
1010 : // Convert to flow at given speed.
1011 0 : F[0] *= (propN.density / RhoAir) * control;
1012 : // Set derivative w/r pressure drop (-).
1013 0 : if (control >= TranRat) {
1014 0 : DF[0] = -control / DPDF;
1015 : } else {
1016 0 : DF[0] = -1.0 / DPDF;
1017 : }
1018 0 : return 1;
1019 : }
1020 :
1021 0 : int Damper::calculate([[maybe_unused]] EnergyPlusData &state,
1022 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
1023 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1024 : int const i, // Linkage number
1025 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1026 : [[maybe_unused]] const Real64 control, // Element control signal
1027 : const AirState &propN, // Node 1 properties
1028 : const AirState &propM, // Node 2 properties
1029 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1030 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1031 : )
1032 : {
1033 :
1034 : // SUBROUTINE INFORMATION:
1035 : // AUTHOR George Walton
1036 : // DATE WRITTEN Extracted from AIRNET
1037 : // MODIFIED Lixing Gu, 2/1/04
1038 : // Revised the subroutine to meet E+ needs
1039 : // MODIFIED Lixing Gu, 6/8/05
1040 : // RE-ENGINEERED na
1041 :
1042 : // PURPOSE OF THIS SUBROUTINE:
1043 : // This subroutine solves airflow for a Controlled power law resistance airflow component (damper)
1044 :
1045 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1046 : Real64 C;
1047 :
1048 : // Formats
1049 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
1050 :
1051 0 : C = control;
1052 0 : if (C < FlowMin) C = FlowMin;
1053 0 : if (C > FlowMax) C = FlowMax;
1054 0 : C = A0 + C * (A1 + C * (A2 + C * A3));
1055 : // if (LIST >= 4)
1056 : // gio::write(Unit21, Format_901) << " Dmp:" << i << AFECTL(i) << DisSysCompDamperData(CompNum).FlowMin
1057 : // << DisSysCompDamperData(CompNum).FlowMax << C;
1058 0 : if (LFLAG || std::abs(PDROP) <= LTP) {
1059 : // Laminar flow.
1060 0 : if (PDROP >= 0.0) {
1061 0 : DF[0] = C * LamFlow * propN.density / propN.viscosity;
1062 : } else {
1063 0 : DF[0] = C * LamFlow * propM.density / propM.viscosity;
1064 : }
1065 0 : F[0] = DF[0] * PDROP;
1066 : } else {
1067 : // Turbulent flow.
1068 0 : if (PDROP >= 0.0) {
1069 0 : F[0] = C * TurFlow * propN.sqrt_density * std::pow(PDROP, FlowExpo);
1070 : } else {
1071 0 : F[0] = -C * TurFlow * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
1072 : }
1073 0 : DF[0] = F[0] * FlowExpo / PDROP;
1074 : }
1075 0 : return 1;
1076 : }
1077 :
1078 0 : int Damper::calculate([[maybe_unused]] EnergyPlusData &state,
1079 : const Real64 PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1080 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1081 : const Real64 control, // Element control signal
1082 : const AirState &propN, // Node 1 properties
1083 : const AirState &propM, // Node 2 properties
1084 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1085 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1086 : )
1087 : {
1088 :
1089 : // SUBROUTINE INFORMATION:
1090 : // AUTHOR George Walton
1091 : // DATE WRITTEN Extracted from AIRNET
1092 : // MODIFIED Lixing Gu, 2/1/04
1093 : // Revised the subroutine to meet E+ needs
1094 : // MODIFIED Lixing Gu, 6/8/05
1095 : // RE-ENGINEERED na
1096 :
1097 : // PURPOSE OF THIS SUBROUTINE:
1098 : // This subroutine solves airflow for a Controlled power law resistance airflow component (damper)
1099 :
1100 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1101 : Real64 C;
1102 :
1103 : // Formats
1104 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
1105 :
1106 0 : C = control;
1107 0 : if (C < FlowMin) C = FlowMin;
1108 0 : if (C > FlowMax) C = FlowMax;
1109 0 : C = A0 + C * (A1 + C * (A2 + C * A3));
1110 : // if (LIST >= 4)
1111 : // gio::write(Unit21, Format_901) << " Dmp:" << i << AFECTL(i) << DisSysCompDamperData(CompNum).FlowMin
1112 : // << DisSysCompDamperData(CompNum).FlowMax << C;
1113 0 : if (std::abs(PDROP) <= LTP) {
1114 : // Laminar flow.
1115 0 : if (PDROP >= 0.0) {
1116 0 : DF[0] = C * LamFlow * propN.density / propN.viscosity;
1117 : } else {
1118 0 : DF[0] = C * LamFlow * propM.density / propM.viscosity;
1119 : }
1120 0 : F[0] = DF[0] * PDROP;
1121 : } else {
1122 : // Turbulent flow.
1123 0 : if (PDROP >= 0.0) {
1124 0 : F[0] = C * TurFlow * propN.sqrt_density * std::pow(PDROP, FlowExpo);
1125 : } else {
1126 0 : F[0] = -C * TurFlow * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
1127 : }
1128 0 : DF[0] = F[0] * FlowExpo / PDROP;
1129 : }
1130 0 : return 1;
1131 : }
1132 :
1133 653732 : int EffectiveLeakageRatio::calculate([[maybe_unused]] EnergyPlusData &state,
1134 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
1135 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1136 : [[maybe_unused]] int const i, // Linkage number
1137 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1138 : [[maybe_unused]] const Real64 control, // Element control signal
1139 : const AirState &propN, // Node 1 properties
1140 : const AirState &propM, // Node 2 properties
1141 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1142 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1143 : )
1144 : {
1145 :
1146 : // SUBROUTINE INFORMATION:
1147 : // AUTHOR George Walton
1148 : // DATE WRITTEN Extracted from AIRNET
1149 : // MODIFIED Lixing Gu, 2/1/04
1150 : // Revised the subroutine to meet E+ needs
1151 : // MODIFIED Lixing Gu, 6/8/05
1152 : // RE-ENGINEERED na
1153 :
1154 : // PURPOSE OF THIS SUBROUTINE:
1155 : // This subroutine solves airflow for a Effective leakage ratio component
1156 :
1157 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1158 : Real64 CDM;
1159 : Real64 FL;
1160 : Real64 FT;
1161 : Real64 FlowCoef;
1162 :
1163 : // Formats
1164 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
1165 :
1166 : // Get component properties
1167 653732 : FlowCoef = ELR * FlowRate / propN.density * std::pow(RefPres, -FlowExpo);
1168 :
1169 653732 : if (LFLAG) {
1170 : // Initialization by linear relation.
1171 0 : if (PDROP >= 0.0) {
1172 0 : DF[0] = FlowCoef * propN.density / propN.viscosity;
1173 : } else {
1174 0 : DF[0] = FlowCoef * propM.density / propM.viscosity;
1175 : }
1176 0 : F[0] = -DF[0] * PDROP;
1177 : } else {
1178 : // Standard calculation.
1179 653732 : if (PDROP >= 0.0) {
1180 : // Flow in positive direction.
1181 : // Laminar flow.
1182 563553 : CDM = FlowCoef * propN.density / propN.viscosity;
1183 563553 : FL = CDM * PDROP;
1184 : // Turbulent flow.
1185 563553 : if (FlowExpo == 0.5) {
1186 0 : FT = FlowCoef * propN.sqrt_density * std::sqrt(PDROP);
1187 : } else {
1188 563553 : FT = FlowCoef * propN.sqrt_density * std::pow(PDROP, FlowExpo);
1189 : }
1190 : } else {
1191 : // Flow in negative direction.
1192 : // Laminar flow.
1193 90179 : CDM = FlowCoef * propM.density / propM.viscosity;
1194 90179 : FL = CDM * PDROP;
1195 : // Turbulent flow.
1196 90179 : if (FlowExpo == 0.5) {
1197 0 : FT = -FlowCoef * propM.sqrt_density * std::sqrt(-PDROP);
1198 : } else {
1199 90179 : FT = -FlowCoef * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
1200 : }
1201 : }
1202 : // Select laminar or turbulent flow.
1203 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " plr: " << i << PDROP << FL << FT;
1204 653732 : if (std::abs(FL) <= std::abs(FT)) {
1205 50 : F[0] = FL;
1206 50 : DF[0] = CDM;
1207 : } else {
1208 653682 : F[0] = FT;
1209 653682 : DF[0] = FT * FlowExpo / PDROP;
1210 : }
1211 : }
1212 653732 : return 1;
1213 : }
1214 :
1215 0 : int EffectiveLeakageRatio::calculate([[maybe_unused]] EnergyPlusData &state,
1216 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1217 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1218 : [[maybe_unused]] const Real64 control, // Element control signal
1219 : const AirState &propN, // Node 1 properties
1220 : const AirState &propM, // Node 2 properties
1221 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1222 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1223 : )
1224 : {
1225 :
1226 : // SUBROUTINE INFORMATION:
1227 : // AUTHOR George Walton
1228 : // DATE WRITTEN Extracted from AIRNET
1229 : // MODIFIED Lixing Gu, 2/1/04
1230 : // Revised the subroutine to meet E+ needs
1231 : // MODIFIED Lixing Gu, 6/8/05
1232 : // RE-ENGINEERED na
1233 :
1234 : // PURPOSE OF THIS SUBROUTINE:
1235 : // This subroutine solves airflow for a Effective leakage ratio component
1236 :
1237 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1238 : Real64 CDM;
1239 : Real64 FL;
1240 : Real64 FT;
1241 : Real64 FlowCoef;
1242 :
1243 : // Formats
1244 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
1245 :
1246 : // Get component properties
1247 0 : FlowCoef = ELR * FlowRate / propN.density * std::pow(RefPres, -FlowExpo);
1248 :
1249 : // Standard calculation.
1250 0 : if (PDROP >= 0.0) {
1251 : // Flow in positive direction.
1252 : // Laminar flow.
1253 0 : CDM = FlowCoef * propN.density / propN.viscosity;
1254 0 : FL = CDM * PDROP;
1255 : // Turbulent flow.
1256 0 : if (FlowExpo == 0.5) {
1257 0 : FT = FlowCoef * propN.sqrt_density * std::sqrt(PDROP);
1258 : } else {
1259 0 : FT = FlowCoef * propN.sqrt_density * std::pow(PDROP, FlowExpo);
1260 : }
1261 : } else {
1262 : // Flow in negative direction.
1263 : // Laminar flow.
1264 0 : CDM = FlowCoef * propM.density / propM.viscosity;
1265 0 : FL = CDM * PDROP;
1266 : // Turbulent flow.
1267 0 : if (FlowExpo == 0.5) {
1268 0 : FT = -FlowCoef * propM.sqrt_density * std::sqrt(-PDROP);
1269 : } else {
1270 0 : FT = -FlowCoef * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
1271 : }
1272 : }
1273 : // Select laminar or turbulent flow.
1274 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " plr: " << i << PDROP << FL << FT;
1275 0 : if (std::abs(FL) <= std::abs(FT)) {
1276 0 : F[0] = FL;
1277 0 : DF[0] = CDM;
1278 : } else {
1279 0 : F[0] = FT;
1280 0 : DF[0] = FT * FlowExpo / PDROP;
1281 : }
1282 0 : return 1;
1283 : }
1284 :
1285 14 : int DetailedOpening::calculate(EnergyPlusData &state,
1286 : [[maybe_unused]] bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
1287 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1288 : int const IL, // Linkage number
1289 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1290 : [[maybe_unused]] const Real64 control, // Element control signal
1291 : [[maybe_unused]] const AirState &propN, // Node 1 properties
1292 : [[maybe_unused]] const AirState &propM, // Node 2 properties
1293 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1294 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1295 : )
1296 : {
1297 :
1298 : // SUBROUTINE INFORMATION:
1299 : // AUTHOR Lixing Gu
1300 : // DATE WRITTEN Oct. 2005
1301 : // MODIFIED na
1302 : // RE-ENGINEERED This subroutine is revised based on a vertical large opening subroutine from COMIS
1303 :
1304 : // PURPOSE OF THIS SUBROUTINE:
1305 : // This subroutine simulates airflow and pressure of a detailed large opening component.
1306 :
1307 : // METHODOLOGY EMPLOYED:
1308 : // Purpose: This routine calculates the massflow and its derivative
1309 : // through a large opening in both flow directions. As input
1310 : // the density profiles RhoProfF/T are required aswell as the
1311 : // effective pressure difference profile DpProfNew, which is the
1312 : // sum of the stack pressure difference profile DpProf and the
1313 : // difference of the actual pressures at reference height. The
1314 : // profiles are calculated in the routine PresProfile.
1315 : // The massflow and its derivative are calculated for each
1316 : // interval representing a step of the pressure difference
1317 : // profile. The total flow and derivative are obtained by
1318 : // summation over the whole opening.
1319 : // The calculation is split into different cases representing
1320 : // different situations of the opening:
1321 : // - closed opening (opening factor = 0): summation of top and
1322 : // bottom crack (crack length = lwmax) plus "integration" over
1323 : // a vertically distributed crack of length (2*lhmax+lextra).
1324 : // - type 1: normal rectangular opening: "integration" over NrInt
1325 : // openings with width actlw and height actlh/NrInt
1326 : // - type 2: horizontally pivoted window: flow direction assumed
1327 : // strictly perpendicular to the plane of the opening
1328 : // -> "integration" over normal rectangular openings at top
1329 : // and bottom of LO plus a rectangular opening in series with two
1330 : // triangular openings in the middle of the LO (most general
1331 : // situation). The geometry is defined by the input parameters
1332 : // actlw(=lwmax), actlh, axisheight, opening angle.
1333 : // Assuming the massflow perpendicular to the opening plane in all
1334 : // cases the ownheightfactor has no influence on the massflow.
1335 :
1336 : // REFERENCES:
1337 : // Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
1338 : // Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990
1339 :
1340 : // USE STATEMENTS:
1341 : // Locals
1342 : // SUBROUTINE ARGUMENT DEFINITIONS:
1343 :
1344 : // SUBROUTINE PARAMETER DEFINITIONS:
1345 14 : Real64 constexpr RealMin(1e-37);
1346 : static Real64 const sqrt_1_2(std::sqrt(1.2));
1347 :
1348 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1349 : static Real64 const sqrt_2(std::sqrt(2.0));
1350 :
1351 : Real64 Width;
1352 : Real64 Height;
1353 :
1354 : Real64 fma12; // massflow in direction "from-to" [kg/s]
1355 : Real64 fma21; // massflow in direction "to-from" [kg/s]
1356 : Real64 dp1fma12; // derivative d fma12 / d Dp [kg/s/Pa]
1357 : Real64 dp1fma21; // derivative d fma21 / d Dp [kg/s/Pa]
1358 14 : Array1D<Real64> DpProfNew(NrInt + 2); // Differential pressure profile for Large Openings, taking into account fixed
1359 : // pressures and actual zone pressures at reference height
1360 : Real64 Fact; // Actual opening factor
1361 : Real64 DifLim; // Limit for the pressure difference where laminarization takes place [Pa]
1362 : Real64 Cfact;
1363 : Real64 FvVeloc;
1364 :
1365 : Real64 ActLh;
1366 : Real64 ActLw;
1367 : Real64 Lextra;
1368 : Real64 Axishght;
1369 : Real64 ActCD;
1370 : Real64 Cs;
1371 : Real64 expn;
1372 : Real64 Type;
1373 : Real64 Interval;
1374 : Real64 fmasum;
1375 : Real64 dfmasum;
1376 : Real64 Prefact;
1377 14 : Array1D<Real64> EvalHghts(NrInt + 2);
1378 : Real64 h2;
1379 : Real64 h4;
1380 : Real64 alpha;
1381 : Real64 rholink;
1382 : Real64 c1;
1383 : Real64 c2;
1384 : Real64 DpZeroOffset;
1385 : Real64 area;
1386 : Real64 WFact;
1387 : Real64 HFact;
1388 : int i;
1389 : int Loc;
1390 : int iNum;
1391 :
1392 : // Get component properties
1393 14 : DifLim = 1.0e-4;
1394 14 : Width = state.afn->MultizoneSurfaceData(IL).Width;
1395 14 : Height = state.afn->MultizoneSurfaceData(IL).Height;
1396 14 : Fact = state.afn->MultizoneSurfaceData(IL).OpenFactor;
1397 14 : Loc = (state.afn->AirflowNetworkLinkageData(IL).DetOpenNum - 1) * (NrInt + 2);
1398 14 : iNum = NumFac;
1399 14 : ActCD = 0.0;
1400 :
1401 14 : if (iNum == 2) {
1402 14 : if (Fact <= OpenFac2) {
1403 14 : WFact = WidthFac1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (WidthFac2 - WidthFac1);
1404 14 : HFact = HeightFac1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (HeightFac2 - HeightFac1);
1405 14 : Cfact = DischCoeff1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (DischCoeff2 - DischCoeff1);
1406 : } else {
1407 0 : ShowFatalError(
1408 : state,
1409 0 : "Open Factor is above the maximum input range for opening factors in AirflowNetwork:MultiZone:Component:DetailedOpening = " +
1410 0 : name);
1411 : }
1412 : }
1413 :
1414 14 : if (iNum == 3) {
1415 0 : if (Fact <= OpenFac2) {
1416 0 : WFact = WidthFac1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (WidthFac2 - WidthFac1);
1417 0 : HFact = HeightFac1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (HeightFac2 - HeightFac1);
1418 0 : Cfact = DischCoeff1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (DischCoeff2 - DischCoeff1);
1419 0 : } else if (Fact <= OpenFac3) {
1420 0 : WFact = WidthFac2 + (Fact - OpenFac2) / (OpenFac3 - OpenFac2) * (WidthFac3 - WidthFac2);
1421 0 : HFact = HeightFac2 + (Fact - OpenFac2) / (OpenFac3 - OpenFac2) * (HeightFac3 - HeightFac2);
1422 0 : Cfact = DischCoeff2 + (Fact - OpenFac2) / (OpenFac3 - OpenFac2) * (DischCoeff3 - DischCoeff2);
1423 : } else {
1424 0 : ShowFatalError(
1425 : state,
1426 0 : "Open Factor is above the maximum input range for opening factors in AirflowNetwork:MultiZone:Component:DetailedOpening = " +
1427 0 : name);
1428 : }
1429 : }
1430 :
1431 14 : if (iNum == 4) {
1432 0 : if (Fact <= OpenFac2) {
1433 0 : WFact = WidthFac1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (WidthFac2 - WidthFac1);
1434 0 : HFact = HeightFac1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (HeightFac2 - HeightFac1);
1435 0 : Cfact = DischCoeff1 + (Fact - OpenFac1) / (OpenFac2 - OpenFac1) * (DischCoeff2 - DischCoeff1);
1436 0 : } else if (Fact <= OpenFac3) {
1437 0 : WFact = WidthFac2 + (Fact - OpenFac2) / (OpenFac3 - OpenFac2) * (WidthFac3 - WidthFac2);
1438 0 : HFact = HeightFac2 + (Fact - OpenFac2) / (OpenFac3 - OpenFac2) * (HeightFac3 - HeightFac2);
1439 0 : Cfact = DischCoeff2 + (Fact - OpenFac2) / (OpenFac3 - OpenFac2) * (DischCoeff3 - DischCoeff2);
1440 0 : } else if (Fact <= OpenFac4) {
1441 0 : WFact = WidthFac3 + (Fact - OpenFac3) / (OpenFac4 - OpenFac3) * (WidthFac4 - WidthFac3);
1442 0 : HFact = HeightFac3 + (Fact - OpenFac3) / (OpenFac4 - OpenFac3) * (HeightFac4 - HeightFac3);
1443 0 : Cfact = DischCoeff3 + (Fact - OpenFac3) / (OpenFac4 - OpenFac3) * (DischCoeff4 - DischCoeff3);
1444 : } else {
1445 0 : ShowFatalError(
1446 : state,
1447 0 : "Open Factor is above the maximum input range for opening factors in AirflowNetwork:MultiZone:Component:DetailedOpening = " +
1448 0 : name);
1449 : }
1450 : }
1451 :
1452 : // calculate DpProfNew
1453 322 : for (i = 1; i <= NrInt + 2; ++i) {
1454 308 : DpProfNew(i) = PDROP + state.afn->dos.DpProf(Loc + i) - state.afn->dos.DpL(IL, 1);
1455 : }
1456 :
1457 : // Get opening data based on the opening factor
1458 14 : if (Fact == 0) {
1459 0 : ActLw = state.afn->MultizoneSurfaceData(IL).Width;
1460 0 : ActLh = state.afn->MultizoneSurfaceData(IL).Height;
1461 0 : Cfact = 0.0;
1462 : } else {
1463 14 : ActLw = state.afn->MultizoneSurfaceData(IL).Width * WFact;
1464 14 : ActLh = state.afn->MultizoneSurfaceData(IL).Height * HFact;
1465 14 : ActCD = Cfact;
1466 : }
1467 :
1468 14 : Cs = FlowCoef;
1469 14 : expn = FlowExpo;
1470 14 : Type = LVOType;
1471 14 : if (Type == 1) {
1472 14 : Lextra = LVOValue;
1473 14 : Axishght = 0.0;
1474 0 : } else if (Type == 2) {
1475 0 : Lextra = 0.0;
1476 0 : Axishght = LVOValue;
1477 0 : ActLw = state.afn->MultizoneSurfaceData(IL).Width;
1478 0 : ActLh = state.afn->MultizoneSurfaceData(IL).Height;
1479 : }
1480 :
1481 : // Add window multiplier with window close
1482 14 : if (state.afn->MultizoneSurfaceData(IL).Multiplier > 1.0) Cs *= state.afn->MultizoneSurfaceData(IL).Multiplier;
1483 : // Add window multiplier with window open
1484 14 : if (Fact > 0.0) {
1485 14 : if (state.afn->MultizoneSurfaceData(IL).Multiplier > 1.0) ActLw *= state.afn->MultizoneSurfaceData(IL).Multiplier;
1486 : // Add recurring warnings
1487 14 : if (ActLw == 0.0) {
1488 0 : ++WidthErrCount;
1489 0 : if (WidthErrCount < 2) {
1490 0 : ShowWarningError(state, "The actual width of the AirflowNetwork:MultiZone:Component:DetailedOpening of " + name + " is 0.");
1491 0 : ShowContinueError(state, "The actual width is set to 1.0E-6 m.");
1492 0 : ShowContinueErrorTimeStamp(state, "Occurrence info:");
1493 : } else {
1494 0 : ShowRecurringWarningErrorAtEnd(state,
1495 0 : "The actual width of the AirflowNetwork:MultiZone:Component:DetailedOpening of " + name +
1496 : " is 0 error continues.",
1497 0 : WidthErrIndex,
1498 : ActLw,
1499 : ActLw);
1500 : }
1501 0 : ActLw = 1.0e-6;
1502 : }
1503 14 : if (ActLh == 0.0) {
1504 0 : ++HeightErrCount;
1505 0 : if (HeightErrCount < 2) {
1506 0 : ShowWarningError(state, "The actual height of the AirflowNetwork:MultiZone:Component:DetailedOpening of " + name + " is 0.");
1507 0 : ShowContinueError(state, "The actual height is set to 1.0E-6 m.");
1508 0 : ShowContinueErrorTimeStamp(state, "Occurrence info:");
1509 : } else {
1510 0 : ShowRecurringWarningErrorAtEnd(state,
1511 0 : "The actual width of the AirflowNetwork:MultiZone:Component:DetailedOpening of " + name +
1512 : " is 0 error continues.",
1513 0 : HeightErrIndex,
1514 : ActLh,
1515 : ActLh);
1516 : }
1517 0 : ActLh = 1.0e-6;
1518 : }
1519 : }
1520 : // Initialization:
1521 14 : int NF(1);
1522 14 : Interval = ActLh / NrInt;
1523 14 : fma12 = 0.0;
1524 14 : fma21 = 0.0;
1525 14 : dp1fma12 = 0.0;
1526 14 : dp1fma21 = 0.0;
1527 :
1528 : // Closed LO
1529 14 : if (Cfact == 0) {
1530 0 : DpZeroOffset = DifLim;
1531 : // bottom crack
1532 0 : if (DpProfNew(1) > 0) {
1533 0 : if (std::abs(DpProfNew(1)) <= DpZeroOffset) {
1534 0 : dfmasum = Cs * ActLw * std::pow(DpZeroOffset, expn) / DpZeroOffset;
1535 0 : fmasum = DpProfNew(1) * dfmasum;
1536 : } else {
1537 0 : fmasum = Cs * ActLw * std::pow(DpProfNew(1), expn);
1538 0 : dfmasum = fmasum * expn / DpProfNew(1);
1539 : }
1540 0 : fma12 += fmasum;
1541 0 : dp1fma12 += dfmasum;
1542 : } else {
1543 0 : if (std::abs(DpProfNew(1)) <= DpZeroOffset) {
1544 0 : dfmasum = -Cs * ActLw * std::pow(DpZeroOffset, expn) / DpZeroOffset;
1545 0 : fmasum = DpProfNew(1) * dfmasum;
1546 : } else {
1547 0 : fmasum = Cs * ActLw * std::pow(-DpProfNew(1), expn);
1548 0 : dfmasum = fmasum * expn / DpProfNew(1);
1549 : }
1550 0 : fma21 += fmasum;
1551 0 : dp1fma21 += dfmasum;
1552 : }
1553 : // top crack
1554 0 : if (DpProfNew(NrInt + 2) > 0) {
1555 0 : if (std::abs(DpProfNew(NrInt + 2)) <= DpZeroOffset) {
1556 0 : dfmasum = Cs * ActLw * std::pow(DpZeroOffset, expn) / DpZeroOffset;
1557 0 : fmasum = DpProfNew(NrInt + 2) * dfmasum;
1558 : } else {
1559 0 : fmasum = Cs * ActLw * std::pow(DpProfNew(NrInt + 2), expn);
1560 0 : dfmasum = fmasum * expn / DpProfNew(NrInt + 2);
1561 : }
1562 0 : fma12 += fmasum;
1563 0 : dp1fma12 += dfmasum;
1564 : } else {
1565 0 : if (std::abs(DpProfNew(NrInt + 2)) <= DpZeroOffset) {
1566 0 : dfmasum = -Cs * ActLw * std::pow(DpZeroOffset, expn) / DpZeroOffset;
1567 0 : fmasum = DpProfNew(NrInt + 2) * dfmasum;
1568 : } else {
1569 0 : fmasum = Cs * ActLw * std::pow(-DpProfNew(NrInt + 2), expn);
1570 0 : dfmasum = fmasum * expn / DpProfNew(NrInt + 2);
1571 : }
1572 0 : fma21 += fmasum;
1573 0 : dp1fma21 += dfmasum;
1574 : }
1575 : // side and extra cracks
1576 0 : Prefact = Interval * (2 + Lextra / ActLh) * Cs;
1577 0 : for (i = 2; i <= NrInt + 1; ++i) {
1578 0 : if (DpProfNew(i) > 0) {
1579 0 : if (std::abs(DpProfNew(i)) <= DpZeroOffset) {
1580 0 : dfmasum = Prefact * std::pow(DpZeroOffset, expn) / DpZeroOffset;
1581 0 : fmasum = DpProfNew(i) * dfmasum;
1582 : } else {
1583 0 : fmasum = Prefact * std::pow(DpProfNew(i), expn);
1584 0 : dfmasum = fmasum * expn / DpProfNew(i);
1585 : }
1586 0 : fma12 += fmasum;
1587 0 : dp1fma12 += dfmasum;
1588 : } else {
1589 0 : if (std::abs(DpProfNew(i)) <= DpZeroOffset) {
1590 0 : dfmasum = -Prefact * std::pow(DpZeroOffset, expn) / DpZeroOffset;
1591 0 : fmasum = DpProfNew(i) * dfmasum;
1592 : } else {
1593 0 : fmasum = Prefact * std::pow(-DpProfNew(i), expn);
1594 0 : dfmasum = fmasum * expn / DpProfNew(i);
1595 : }
1596 0 : fma21 += fmasum;
1597 0 : dp1fma21 += dfmasum;
1598 : }
1599 : }
1600 : }
1601 :
1602 : // Open LO, type 1
1603 14 : if ((Cfact != 0) && (Type == 1)) {
1604 14 : DpZeroOffset = DifLim * 1e-3;
1605 14 : Prefact = ActLw * ActCD * Interval * sqrt_2;
1606 294 : for (i = 2; i <= NrInt + 1; ++i) {
1607 280 : if (DpProfNew(i) > 0) {
1608 160 : if (std::abs(DpProfNew(i)) <= DpZeroOffset) {
1609 0 : dfmasum = std::sqrt(state.afn->dos.RhoProfF(Loc + i) * DpZeroOffset) / DpZeroOffset;
1610 0 : fmasum = DpProfNew(i) * dfmasum;
1611 : } else {
1612 160 : fmasum = std::sqrt(state.afn->dos.RhoProfF(Loc + i) * DpProfNew(i));
1613 160 : dfmasum = 0.5 * fmasum / DpProfNew(i);
1614 : }
1615 160 : fma12 += fmasum;
1616 160 : dp1fma12 += dfmasum;
1617 : } else {
1618 120 : if (std::abs(DpProfNew(i)) <= DpZeroOffset) {
1619 0 : dfmasum = -std::sqrt(state.afn->dos.RhoProfT(Loc + i) * DpZeroOffset) / DpZeroOffset;
1620 0 : fmasum = DpProfNew(i) * dfmasum;
1621 : } else {
1622 120 : fmasum = std::sqrt(-state.afn->dos.RhoProfT(Loc + i) * DpProfNew(i));
1623 120 : dfmasum = 0.5 * fmasum / DpProfNew(i);
1624 : }
1625 120 : fma21 += fmasum;
1626 120 : dp1fma21 += dfmasum;
1627 : }
1628 : }
1629 :
1630 14 : fma12 *= Prefact;
1631 14 : fma21 *= Prefact;
1632 14 : dp1fma12 *= Prefact;
1633 14 : dp1fma21 *= Prefact;
1634 : }
1635 :
1636 : // Open LO, type 2
1637 14 : if ((Cfact != 0) && (Type == 2)) {
1638 : // Initialization
1639 0 : DpZeroOffset = DifLim * 1e-3;
1640 : // New definition for opening factors for LVO type 2: opening angle = 90 degrees --> opening factor = 1.0
1641 : // should be PIOvr2 in below?
1642 0 : alpha = Fact * Constant::PiOvr2;
1643 0 : Real64 const cos_alpha(std::cos(alpha));
1644 0 : Real64 const tan_alpha(std::tan(alpha));
1645 0 : h2 = Axishght * (1.0 - cos_alpha);
1646 0 : h4 = Axishght + (ActLh - Axishght) * cos_alpha;
1647 0 : EvalHghts(1) = 0.0;
1648 0 : EvalHghts(NrInt + 2) = ActLh;
1649 : // New definition for opening factors for LVO type 2: pening angle = 90 degrees --> opening factor = 1.0
1650 0 : if (Fact == 1.0) {
1651 0 : h2 = Axishght;
1652 0 : h4 = Axishght;
1653 : }
1654 :
1655 0 : for (i = 2; i <= NrInt + 1; ++i) {
1656 0 : EvalHghts(i) = Interval * (i - 1.5);
1657 : }
1658 :
1659 : // Calculation of massflow and its derivative
1660 0 : for (i = 2; i <= NrInt + 1; ++i) {
1661 0 : if (DpProfNew(i) > 0) {
1662 0 : rholink = state.afn->dos.RhoProfF(Loc + i);
1663 : } else {
1664 0 : rholink = state.afn->dos.RhoProfT(Loc + i);
1665 : }
1666 :
1667 0 : if ((EvalHghts(i) <= h2) || (EvalHghts(i) >= h4)) {
1668 0 : if (std::abs(DpProfNew(i)) <= DpZeroOffset) {
1669 0 : dfmasum = ActCD * ActLw * Interval * std::sqrt(2.0 * rholink * DpZeroOffset) / DpZeroOffset * sign(1, DpProfNew(i));
1670 0 : fmasum = DpProfNew(i) * dfmasum;
1671 : } else {
1672 0 : fmasum = ActCD * ActLw * Interval * std::sqrt(2.0 * rholink * std::abs(DpProfNew(i)));
1673 0 : dfmasum = 0.5 * fmasum / DpProfNew(i);
1674 : }
1675 : } else {
1676 : // triangular opening at the side of LO
1677 0 : c1 = ActCD * ActLw * Interval * std::sqrt(2.0 * rholink);
1678 0 : c2 = 2 * ActCD * std::abs(Axishght - EvalHghts(i)) * tan_alpha * Interval * std::sqrt(2.0 * rholink);
1679 0 : if ((c1 != 0) && (c2 != 0)) {
1680 0 : if (std::abs(DpProfNew(i)) <= DpZeroOffset) {
1681 0 : dfmasum = std::sqrt(DpZeroOffset / (1 / c1 / c1 + 1 / c2 / c2)) / DpZeroOffset * sign(1, DpProfNew(i));
1682 0 : fmasum = DpProfNew(i) * dfmasum;
1683 : } else {
1684 0 : fmasum = std::sqrt(std::abs(DpProfNew(i)) / (1 / c1 / c1 + 1 / c2 / c2));
1685 0 : dfmasum = 0.5 * fmasum / DpProfNew(i);
1686 : }
1687 : } else {
1688 0 : fmasum = 0.0;
1689 0 : dfmasum = 0.0;
1690 : }
1691 : }
1692 :
1693 0 : if (DpProfNew(i) > 0) {
1694 0 : fma12 += fmasum;
1695 0 : dp1fma12 += dfmasum;
1696 : } else {
1697 0 : fma21 += fmasum;
1698 0 : dp1fma21 += dfmasum;
1699 : }
1700 : }
1701 : }
1702 :
1703 : // Calculate some velocity in the large opening
1704 14 : area = ActLh * ActLw * ActCD;
1705 14 : if (area > (Cs + RealMin)) {
1706 14 : if (area > RealMin) {
1707 14 : FvVeloc = (fma21 + fma12) / area;
1708 : } else {
1709 0 : FvVeloc = 0.0;
1710 : }
1711 : } else {
1712 : // here the average velocity over the full area, may blow half in half out.
1713 : // velocity= Fva/Nett area=Fma/Rho/(Cm/( (2**N)* SQRT(1.2) ) )
1714 0 : if (Cs > 0.0) {
1715 : // get the average Rho for this closed window
1716 0 : for (i = 2; i <= NrInt + 1; ++i) {
1717 : // rholink = 0.0;
1718 : // if (DpProfNew(i) > 0) {
1719 : // rholink = state.afn->dos.RhoProfF(Loc + i);
1720 : //} else {
1721 : // rholink = state.afn->dos.RhoProfT(Loc + i);
1722 : //}
1723 : // rholink /= NrInt;
1724 0 : rholink = 1.2;
1725 : }
1726 0 : FvVeloc = (fma21 + fma12) * std::pow(2.0, expn) * sqrt_1_2 / (rholink * Cs);
1727 : } else {
1728 0 : FvVeloc = 0.0;
1729 : }
1730 : }
1731 :
1732 : // Output mass flow rates and associated derivatives
1733 14 : F[0] = fma12 - fma21;
1734 14 : DF[0] = dp1fma12 - dp1fma21;
1735 14 : F[1] = 0.0;
1736 14 : if (fma12 != 0.0 && fma21 != 0.0) {
1737 8 : F[1] = fma21;
1738 : }
1739 14 : DF[1] = 0.0;
1740 14 : return NF;
1741 14 : }
1742 :
1743 0 : int SimpleOpening::calculate(EnergyPlusData &state,
1744 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
1745 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1746 : int const i, // Linkage number
1747 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1748 : [[maybe_unused]] const Real64 control, // Element control signal
1749 : const AirState &propN, // Node 1 properties
1750 : const AirState &propM, // Node 2 properties
1751 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1752 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1753 : )
1754 : {
1755 :
1756 : // SUBROUTINE INFORMATION:
1757 : // AUTHOR George Walton
1758 : // DATE WRITTEN Extracted from AIRNET
1759 : // MODIFIED Lixing Gu, 2/1/04
1760 : // Revised the subroutine to meet E+ needs
1761 : // MODIFIED Lixing Gu, 6/8/05
1762 : // RE-ENGINEERED na
1763 :
1764 : // PURPOSE OF THIS SUBROUTINE:
1765 : // This subroutine solves airflow for a Doorway airflow component using standard interface.
1766 : // A doorway may have two-way airflows. Heights measured relative to the bottom of the door.
1767 :
1768 : // SUBROUTINE PARAMETER DEFINITIONS:
1769 : // Replace this with the value from numbers when we have C++20
1770 0 : Real64 constexpr SQRT2(1.414213562373095);
1771 :
1772 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1773 : Real64 DPMID; // pressure drop at mid-height of doorway.
1774 : Real64 C;
1775 : Real64 DF0; // derivative factor at the bottom of the door.
1776 : Real64 DFH; // derivative factor at the top of the door.
1777 : Real64 DRHO; // difference in air densities between rooms.
1778 : Real64 GDRHO;
1779 : Real64 F0; // flow factor at the bottom of the door.
1780 : Real64 FH; // flow factor at the top of the door.
1781 : Real64 Y; // height of neutral plane rel. to bottom of door (m).
1782 : Real64 coeff;
1783 : Real64 Width;
1784 : Real64 Height;
1785 : Real64 OpenFactor;
1786 :
1787 0 : int NF(1);
1788 :
1789 : // Formats
1790 : // static gio::Fmt Format_900("(A5,9X,4E16.7)");
1791 : // static gio::Fmt Format_903("(A5,3I3,4E16.7)");
1792 :
1793 0 : Width = state.afn->MultizoneSurfaceData(i).Width;
1794 0 : Height = state.afn->MultizoneSurfaceData(i).Height;
1795 0 : coeff = FlowCoef * 2.0 * (Width + Height);
1796 0 : OpenFactor = state.afn->MultizoneSurfaceData(i).OpenFactor;
1797 0 : if (OpenFactor > 0.0) {
1798 0 : Width *= OpenFactor;
1799 0 : if (state.dataSurface->Surface(state.afn->MultizoneSurfaceData(i).SurfNum).Tilt < 90.0) {
1800 0 : Height *= state.dataSurface->Surface(state.afn->MultizoneSurfaceData(i).SurfNum).SinTilt;
1801 : }
1802 : }
1803 :
1804 0 : if (PDROP >= 0.0) {
1805 0 : coeff /= propN.sqrt_density;
1806 : } else {
1807 0 : coeff /= propM.sqrt_density;
1808 : }
1809 :
1810 : // Add window multiplier with window close
1811 0 : if (state.afn->MultizoneSurfaceData(i).Multiplier > 1.0) coeff *= state.afn->MultizoneSurfaceData(i).Multiplier;
1812 : // Add window multiplier with window open
1813 0 : if (OpenFactor > 0.0) {
1814 0 : if (state.afn->MultizoneSurfaceData(i).Multiplier > 1.0) Width *= state.afn->MultizoneSurfaceData(i).Multiplier;
1815 : }
1816 :
1817 0 : DRHO = propN.density - propM.density;
1818 0 : GDRHO = 9.8 * DRHO;
1819 : // if (LIST >= 4) gio::write(Unit21, Format_903) << " DOR:" << i << n << m << PDROP << std::abs(DRHO) << MinRhoDiff;
1820 0 : if (OpenFactor == 0.0) {
1821 0 : generic_crack(coeff, FlowExpo, LFLAG, PDROP, propN, propM, F, DF);
1822 0 : return 1;
1823 : }
1824 0 : if (std::abs(DRHO) < MinRhoDiff || LFLAG) {
1825 0 : DPMID = PDROP - 0.5 * Height * GDRHO;
1826 : // Initialization or identical temps: treat as one-way flow.
1827 0 : NF = 1;
1828 0 : generic_crack(coeff, FlowExpo, LFLAG, DPMID, propN, propM, F, DF);
1829 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " Drs:" << DPMID << F[0] << DF[0];
1830 : } else {
1831 : // Possible two-way flow:
1832 0 : Y = PDROP / GDRHO;
1833 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " DrY:" << PDROP << GDRHO << Y;
1834 : // F0 = lower flow, FH = upper flow.
1835 0 : C = SQRT2 * Width * DischCoeff;
1836 0 : DF0 = C * std::sqrt(std::abs(PDROP)) / std::abs(GDRHO);
1837 : // F0 = 0.666667d0*C*SQRT(ABS(GDRHO*Y))*ABS(Y)
1838 0 : F0 = (2.0 / 3.0) * C * std::sqrt(std::abs(GDRHO * Y)) * std::abs(Y);
1839 0 : DFH = C * std::sqrt(std::abs((Height - Y) / GDRHO));
1840 : // FH = 0.666667d0*DFH*ABS(GDRHO*(Height-Y))
1841 0 : FH = (2.0 / 3.0) * DFH * std::abs(GDRHO * (Height - Y));
1842 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " DrF:" << F0 << DF0 << FH << DFH;
1843 0 : if (Y <= 0.0) {
1844 : // One-way flow (negative).
1845 0 : if (DRHO >= 0.0) {
1846 0 : F[0] = -propM.sqrt_density * std::abs(FH - F0);
1847 0 : DF[0] = propM.sqrt_density * std::abs(DFH - DF0);
1848 : } else {
1849 0 : F[0] = propN.sqrt_density * std::abs(FH - F0);
1850 0 : DF[0] = propN.sqrt_density * std::abs(DFH - DF0);
1851 : }
1852 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " Dr1:" << C << F[0] << DF[0];
1853 0 : } else if (Y >= Height) {
1854 : // One-way flow (positive).
1855 0 : if (DRHO >= 0.0) {
1856 0 : F[0] = propN.sqrt_density * std::abs(FH - F0);
1857 0 : DF[0] = propN.sqrt_density * std::abs(DFH - DF0);
1858 : } else {
1859 0 : F[0] = -propM.sqrt_density * std::abs(FH - F0);
1860 0 : DF[0] = propM.sqrt_density * std::abs(DFH - DF0);
1861 : }
1862 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " Dr2:" << C << F[0] << DF[0];
1863 : } else {
1864 : // Two-way flow.
1865 0 : NF = 2;
1866 0 : if (DRHO >= 0.0) {
1867 0 : F[0] = -propM.sqrt_density * FH;
1868 0 : DF[0] = propM.sqrt_density * DFH;
1869 0 : F[1] = propN.sqrt_density * F0;
1870 0 : DF[1] = propN.sqrt_density * DF0;
1871 : } else {
1872 0 : F[0] = propN.sqrt_density * FH;
1873 0 : DF[0] = propN.sqrt_density * DFH;
1874 0 : F[1] = -propM.sqrt_density * F0;
1875 0 : DF[1] = propM.sqrt_density * DF0;
1876 : }
1877 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " Dr3:" << C << F[0] << DF[0];
1878 : // if (LIST >= 4) gio::write(Unit21, Format_900) << " Dr4:" << C << F[1] << DF[1];
1879 : }
1880 : }
1881 0 : return NF;
1882 : }
1883 :
1884 132 : int ConstantPressureDrop::calculate([[maybe_unused]] EnergyPlusData &state,
1885 : [[maybe_unused]] bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
1886 : const Real64 PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1887 : int const i, // Linkage number
1888 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1889 : [[maybe_unused]] const Real64 control, // Element control signal
1890 : const AirState &propN, // Node 1 properties
1891 : [[maybe_unused]] const AirState &propM, // Node 2 properties
1892 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1893 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1894 : )
1895 : {
1896 :
1897 : // SUBROUTINE INFORMATION:
1898 : // AUTHOR George Walton
1899 : // DATE WRITTEN Extracted from AIRNET
1900 : // MODIFIED Lixing Gu, 2/1/04
1901 : // Revised the subroutine to meet E+ needs
1902 : // MODIFIED Lixing Gu, 6/8/05
1903 :
1904 : // PURPOSE OF THIS SUBROUTINE:
1905 : // This subroutine solves airflow for a Constant pressure drop component
1906 :
1907 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1908 : // Real64 Co;
1909 :
1910 132 : int n = state.afn->AirflowNetworkLinkageData(i).NodeNums[0];
1911 132 : int m = state.afn->AirflowNetworkLinkageData(i).NodeNums[1];
1912 :
1913 : // Get component properties
1914 : // A = Cross section area [m2]
1915 : // DP = Pressure difference across the element [Pa]
1916 :
1917 132 : if (PDROP == 0.0) {
1918 12 : F[0] = std::sqrt(2.0 * propN.density) * A * std::sqrt(DP);
1919 12 : DF[0] = 0.5 * F[0] / DP;
1920 : } else {
1921 2476 : for (int k = 1; k <= state.afn->ActualNumOfLinks; ++k) {
1922 2476 : if (state.afn->AirflowNetworkLinkageData(k).NodeNums[1] == n) {
1923 120 : F[0] = state.afn->AFLOW(k);
1924 120 : break;
1925 : }
1926 : }
1927 120 : state.afn->PZ(m) = state.afn->PZ(n) - DP;
1928 : // Co = F[0] / DP; // can be deleted since it is not used
1929 120 : DF[0] = 10.e10;
1930 : }
1931 132 : return 1;
1932 : }
1933 :
1934 3536826 : int EffectiveLeakageArea::calculate([[maybe_unused]] EnergyPlusData &state,
1935 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
1936 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
1937 : [[maybe_unused]] int const i, // Linkage number
1938 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
1939 : [[maybe_unused]] const Real64 control, // Element control signal
1940 : const AirState &propN, // Node 1 properties
1941 : const AirState &propM, // Node 2 properties
1942 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
1943 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
1944 : )
1945 : {
1946 :
1947 : // SUBROUTINE INFORMATION:
1948 : // AUTHOR George Walton
1949 : // DATE WRITTEN Extracted from AIRNET
1950 : // MODIFIED Lixing Gu, 2/1/04
1951 : // Revised the subroutine to meet E+ needs
1952 : // MODIFIED Lixing Gu, 6/8/05
1953 : // RE-ENGINEERED na
1954 :
1955 : // PURPOSE OF THIS SUBROUTINE:
1956 : // This subroutine solves airflow for a Surface effective leakage area component
1957 :
1958 : // SUBROUTINE PARAMETER DEFINITIONS:
1959 : static Real64 const sqrt_2(std::sqrt(2.0));
1960 :
1961 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1962 : Real64 CDM;
1963 : Real64 FL;
1964 : Real64 FT;
1965 : Real64 FlowCoef;
1966 :
1967 : // Formats
1968 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
1969 :
1970 : // Get component properties
1971 3536826 : FlowCoef = ELA * DischCoeff * sqrt_2 * std::pow(RefDeltaP, 0.5 - FlowExpo);
1972 :
1973 3536826 : if (LFLAG) {
1974 : // Initialization by linear relation.
1975 0 : if (PDROP >= 0.0) {
1976 0 : DF[0] = FlowCoef * propN.density / propN.viscosity;
1977 : } else {
1978 0 : DF[0] = FlowCoef * propM.density / propM.viscosity;
1979 : }
1980 0 : F[0] = -DF[0] * PDROP;
1981 : } else {
1982 : // Standard calculation.
1983 3536826 : if (PDROP >= 0.0) {
1984 : // Flow in positive direction.
1985 : // Laminar flow.
1986 1813488 : CDM = FlowCoef * propN.density / propN.viscosity;
1987 1813488 : FL = CDM * PDROP;
1988 : // Turbulent flow.
1989 1813488 : if (FlowExpo == 0.5) {
1990 0 : FT = FlowCoef * propN.sqrt_density * std::sqrt(PDROP);
1991 : } else {
1992 1813488 : FT = FlowCoef * propN.sqrt_density * std::pow(PDROP, FlowExpo);
1993 : }
1994 : } else {
1995 : // Flow in negative direction.
1996 : // Laminar flow.
1997 1723338 : CDM = FlowCoef * propM.density / propM.viscosity;
1998 1723338 : FL = CDM * PDROP;
1999 : // Turbulent flow.
2000 1723338 : if (FlowExpo == 0.5) {
2001 0 : FT = -FlowCoef * propM.sqrt_density * std::sqrt(-PDROP);
2002 : } else {
2003 1723338 : FT = -FlowCoef * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
2004 : }
2005 : }
2006 : // Select laminar or turbulent flow.
2007 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " plr: " << i << PDROP << FL << FT;
2008 3536826 : if (std::abs(FL) <= std::abs(FT)) {
2009 168 : F[0] = FL;
2010 168 : DF[0] = CDM;
2011 : } else {
2012 3536658 : F[0] = FT;
2013 3536658 : DF[0] = FT * FlowExpo / PDROP;
2014 : }
2015 : }
2016 3536826 : return 1;
2017 : }
2018 :
2019 0 : int EffectiveLeakageArea::calculate([[maybe_unused]] EnergyPlusData &state,
2020 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2021 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2022 : [[maybe_unused]] const Real64 control, // Element control signal
2023 : const AirState &propN, // Node 1 properties
2024 : const AirState &propM, // Node 2 properties
2025 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2026 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2027 : )
2028 : {
2029 :
2030 : // SUBROUTINE INFORMATION:
2031 : // AUTHOR George Walton
2032 : // DATE WRITTEN Extracted from AIRNET
2033 : // MODIFIED Lixing Gu, 2/1/04
2034 : // Revised the subroutine to meet E+ needs
2035 : // MODIFIED Lixing Gu, 6/8/05
2036 : // RE-ENGINEERED na
2037 :
2038 : // PURPOSE OF THIS SUBROUTINE:
2039 : // This subroutine solves airflow for a Surface effective leakage area component
2040 :
2041 : // SUBROUTINE PARAMETER DEFINITIONS:
2042 : static Real64 const sqrt_2(std::sqrt(2.0));
2043 :
2044 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2045 : Real64 CDM;
2046 : Real64 FL;
2047 : Real64 FT;
2048 : Real64 FlowCoef;
2049 :
2050 : // Formats
2051 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
2052 :
2053 : // Get component properties
2054 0 : FlowCoef = ELA * DischCoeff * sqrt_2 * std::pow(RefDeltaP, 0.5 - FlowExpo);
2055 :
2056 : // Standard calculation.
2057 0 : if (PDROP >= 0.0) {
2058 : // Flow in positive direction.
2059 : // Laminar flow.
2060 0 : CDM = FlowCoef * propN.density / propN.viscosity;
2061 0 : FL = CDM * PDROP;
2062 : // Turbulent flow.
2063 0 : if (FlowExpo == 0.5) {
2064 0 : FT = FlowCoef * propN.sqrt_density * std::sqrt(PDROP);
2065 : } else {
2066 0 : FT = FlowCoef * propN.sqrt_density * std::pow(PDROP, FlowExpo);
2067 : }
2068 : } else {
2069 : // Flow in negative direction.
2070 : // Laminar flow.
2071 0 : CDM = FlowCoef * propM.density / propM.viscosity;
2072 0 : FL = CDM * PDROP;
2073 : // Turbulent flow.
2074 0 : if (FlowExpo == 0.5) {
2075 0 : FT = -FlowCoef * propM.sqrt_density * std::sqrt(-PDROP);
2076 : } else {
2077 0 : FT = -FlowCoef * propM.sqrt_density * std::pow(-PDROP, FlowExpo);
2078 : }
2079 : }
2080 : // Select laminar or turbulent flow.
2081 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " plr: " << i << PDROP << FL << FT;
2082 0 : if (std::abs(FL) <= std::abs(FT)) {
2083 0 : F[0] = FL;
2084 0 : DF[0] = CDM;
2085 : } else {
2086 0 : F[0] = FT;
2087 0 : DF[0] = FT * FlowExpo / PDROP;
2088 : }
2089 :
2090 0 : return 1;
2091 : }
2092 :
2093 490145 : int DisSysCompCoilProp::calculate([[maybe_unused]] EnergyPlusData &state,
2094 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
2095 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2096 : [[maybe_unused]] int const i, // Linkage number
2097 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2098 : [[maybe_unused]] const Real64 control, // Element control signal
2099 : const AirState &propN, // Node 1 properties
2100 : const AirState &propM, // Node 2 properties
2101 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2102 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2103 : )
2104 : {
2105 :
2106 : // SUBROUTINE INFORMATION:
2107 : // AUTHOR George Walton
2108 : // DATE WRITTEN Extracted from AIRNET
2109 : // MODIFIED Lixing Gu, 2/1/04
2110 : // Revised the subroutine to meet E+ needs
2111 : // MODIFIED Lixing Gu, 6/8/05
2112 : // RE-ENGINEERED na
2113 :
2114 : // PURPOSE OF THIS SUBROUTINE:
2115 : // This subroutine solves airflow for a coil component
2116 :
2117 : // SUBROUTINE PARAMETER DEFINITIONS:
2118 490145 : Real64 constexpr C(0.868589);
2119 490145 : Real64 constexpr EPS(0.001);
2120 490145 : Real64 constexpr Rough(0.0001);
2121 490145 : Real64 constexpr InitLamCoef(128.0);
2122 490145 : Real64 constexpr LamDynCoef(64.0);
2123 490145 : Real64 constexpr LamFriCoef(0.0001);
2124 490145 : Real64 constexpr TurDynCoef(0.0001);
2125 :
2126 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2127 : Real64 A0;
2128 : Real64 A1;
2129 : Real64 A2;
2130 : Real64 B;
2131 : Real64 D;
2132 : Real64 S2;
2133 : Real64 CDM;
2134 : Real64 FL;
2135 : Real64 FT;
2136 : Real64 FTT;
2137 : Real64 RE;
2138 : Real64 ed;
2139 : Real64 ld;
2140 : Real64 g;
2141 : Real64 AA1;
2142 : Real64 area;
2143 :
2144 : // Formats
2145 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
2146 :
2147 : // Get component properties
2148 : // ed = Rough / DisSysCompCoilData(CompNum).hydraulicDiameter;
2149 490145 : ed = Rough / hydraulicDiameter;
2150 :
2151 490145 : area = square(hydraulicDiameter) * Constant::Pi;
2152 490145 : ld = L / hydraulicDiameter;
2153 490145 : g = 1.14 - 0.868589 * std::log(ed);
2154 490145 : AA1 = g;
2155 :
2156 490145 : if (LFLAG) {
2157 : // Initialization by linear relation.
2158 2 : if (PDROP >= 0.0) {
2159 1 : DF[0] = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * InitLamCoef * ld);
2160 : } else {
2161 1 : DF[0] = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * InitLamCoef * ld);
2162 : }
2163 2 : F[0] = -DF[0] * PDROP;
2164 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwi:" << i << InitLamCoef << F[0] << DF[0];
2165 : } else {
2166 : // Standard calculation.
2167 490143 : if (PDROP >= 0.0) {
2168 : // Flow in positive direction.
2169 : // Laminar flow coefficient !=0
2170 : if (LamFriCoef >= 0.001) {
2171 : A2 = LamFriCoef / (2.0 * propN.density * area * area);
2172 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * area * hydraulicDiameter);
2173 : A0 = -PDROP;
2174 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2175 : FL = (CDM - A1) / (2.0 * A2);
2176 : CDM = 1.0 / CDM;
2177 : } else {
2178 437295 : CDM = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
2179 437295 : FL = CDM * PDROP;
2180 : }
2181 437295 : RE = FL * hydraulicDiameter / (propN.viscosity * area);
2182 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
2183 : // Turbulent flow; test when Re>10.
2184 437295 : if (RE >= 10.0) {
2185 414072 : S2 = std::sqrt(2.0 * propN.density * PDROP) * area;
2186 414072 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2187 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
2188 : while (true) {
2189 1573913 : FT = FTT;
2190 1573913 : B = (9.3 * propN.viscosity * area) / (FT * Rough);
2191 1573913 : D = 1.0 + g * B;
2192 1573913 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2193 1573913 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2194 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
2195 1573913 : if (std::abs(FTT - FT) / FTT < EPS) break;
2196 : }
2197 414072 : FT = FTT;
2198 : } else {
2199 23223 : FT = FL;
2200 : }
2201 : } else {
2202 : // Flow in negative direction.
2203 : // Laminar flow coefficient !=0
2204 : if (LamFriCoef >= 0.001) {
2205 : A2 = LamFriCoef / (2.0 * propM.density * area * area);
2206 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * area * hydraulicDiameter);
2207 : A0 = PDROP;
2208 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2209 : FL = -(CDM - A1) / (2.0 * A2);
2210 : CDM = 1.0 / CDM;
2211 : } else {
2212 52848 : CDM = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
2213 52848 : FL = CDM * PDROP;
2214 : }
2215 52848 : RE = -FL * hydraulicDiameter / (propM.viscosity * area);
2216 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
2217 : // Turbulent flow; test when Re>10.
2218 52848 : if (RE >= 10.0) {
2219 52848 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * area;
2220 52848 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2221 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
2222 : while (true) {
2223 136306 : FT = FTT;
2224 136306 : B = (9.3 * propM.viscosity * area) / (FT * Rough);
2225 136306 : D = 1.0 + g * B;
2226 136306 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2227 136306 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2228 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
2229 136306 : if (std::abs(FTT - FT) / FTT < EPS) break;
2230 : }
2231 52848 : FT = -FTT;
2232 : } else {
2233 0 : FT = FL;
2234 : }
2235 : }
2236 : // Select laminar or turbulent flow.
2237 490143 : if (std::abs(FL) <= std::abs(FT)) {
2238 23224 : F[0] = FL;
2239 23224 : DF[0] = CDM;
2240 : } else {
2241 466919 : F[0] = FT;
2242 466919 : DF[0] = 0.5 * FT / PDROP;
2243 : }
2244 : }
2245 490145 : return 1;
2246 : }
2247 :
2248 0 : int DisSysCompCoilProp::calculate([[maybe_unused]] EnergyPlusData &state,
2249 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2250 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2251 : [[maybe_unused]] const Real64 control, // Element control signal
2252 : const AirState &propN, // Node 1 properties
2253 : const AirState &propM, // Node 2 properties
2254 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2255 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2256 : )
2257 : {
2258 :
2259 : // SUBROUTINE INFORMATION:
2260 : // AUTHOR George Walton
2261 : // DATE WRITTEN Extracted from AIRNET
2262 : // MODIFIED Lixing Gu, 2/1/04
2263 : // Revised the subroutine to meet E+ needs
2264 : // MODIFIED Lixing Gu, 6/8/05
2265 : // RE-ENGINEERED na
2266 :
2267 : // PURPOSE OF THIS SUBROUTINE:
2268 : // This subroutine solves airflow for a coil component
2269 :
2270 : // SUBROUTINE PARAMETER DEFINITIONS:
2271 0 : Real64 constexpr C(0.868589);
2272 0 : Real64 constexpr EPS(0.001);
2273 0 : Real64 constexpr Rough(0.0001);
2274 0 : Real64 constexpr LamDynCoef(64.0);
2275 0 : Real64 constexpr LamFriCoef(0.0001);
2276 0 : Real64 constexpr TurDynCoef(0.0001);
2277 :
2278 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2279 : Real64 A0;
2280 : Real64 A1;
2281 : Real64 A2;
2282 : Real64 B;
2283 : Real64 D;
2284 : Real64 S2;
2285 : Real64 CDM;
2286 : Real64 FL;
2287 : Real64 FT;
2288 : Real64 FTT;
2289 : Real64 RE;
2290 : Real64 ed;
2291 : Real64 ld;
2292 : Real64 g;
2293 : Real64 AA1;
2294 : Real64 area;
2295 :
2296 : // Formats
2297 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
2298 :
2299 : // Get component properties
2300 : // ed = Rough / DisSysCompCoilData(CompNum).hydraulicDiameter;
2301 0 : ed = Rough / hydraulicDiameter;
2302 :
2303 0 : area = square(hydraulicDiameter) * Constant::Pi;
2304 0 : ld = L / hydraulicDiameter;
2305 0 : g = 1.14 - 0.868589 * std::log(ed);
2306 0 : AA1 = g;
2307 :
2308 : // Standard calculation.
2309 0 : if (PDROP >= 0.0) {
2310 : // Flow in positive direction.
2311 : // Laminar flow coefficient !=0
2312 : if (LamFriCoef >= 0.001) {
2313 : A2 = LamFriCoef / (2.0 * propN.density * area * area);
2314 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * area * hydraulicDiameter);
2315 : A0 = -PDROP;
2316 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2317 : FL = (CDM - A1) / (2.0 * A2);
2318 : CDM = 1.0 / CDM;
2319 : } else {
2320 0 : CDM = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
2321 0 : FL = CDM * PDROP;
2322 : }
2323 0 : RE = FL * hydraulicDiameter / (propN.viscosity * area);
2324 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
2325 : // Turbulent flow; test when Re>10.
2326 0 : if (RE >= 10.0) {
2327 0 : S2 = std::sqrt(2.0 * propN.density * PDROP) * area;
2328 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2329 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
2330 : while (true) {
2331 0 : FT = FTT;
2332 0 : B = (9.3 * propN.viscosity * area) / (FT * Rough);
2333 0 : D = 1.0 + g * B;
2334 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2335 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2336 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
2337 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2338 : }
2339 0 : FT = FTT;
2340 : } else {
2341 0 : FT = FL;
2342 : }
2343 : } else {
2344 : // Flow in negative direction.
2345 : // Laminar flow coefficient !=0
2346 : if (LamFriCoef >= 0.001) {
2347 : A2 = LamFriCoef / (2.0 * propM.density * area * area);
2348 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * area * hydraulicDiameter);
2349 : A0 = PDROP;
2350 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2351 : FL = -(CDM - A1) / (2.0 * A2);
2352 : CDM = 1.0 / CDM;
2353 : } else {
2354 0 : CDM = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
2355 0 : FL = CDM * PDROP;
2356 : }
2357 0 : RE = -FL * hydraulicDiameter / (propM.viscosity * area);
2358 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
2359 : // Turbulent flow; test when Re>10.
2360 0 : if (RE >= 10.0) {
2361 0 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * area;
2362 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2363 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
2364 : while (true) {
2365 0 : FT = FTT;
2366 0 : B = (9.3 * propM.viscosity * area) / (FT * Rough);
2367 0 : D = 1.0 + g * B;
2368 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2369 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2370 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
2371 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2372 : }
2373 0 : FT = -FTT;
2374 : } else {
2375 0 : FT = FL;
2376 : }
2377 : }
2378 : // Select laminar or turbulent flow.
2379 0 : if (std::abs(FL) <= std::abs(FT)) {
2380 0 : F[0] = FL;
2381 0 : DF[0] = CDM;
2382 : } else {
2383 0 : F[0] = FT;
2384 0 : DF[0] = 0.5 * FT / PDROP;
2385 : }
2386 0 : return 1;
2387 : }
2388 :
2389 224 : int DisSysCompTermUnitProp::calculate([[maybe_unused]] EnergyPlusData &state,
2390 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
2391 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2392 : int const i, // Linkage number
2393 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2394 : [[maybe_unused]] const Real64 control, // Element control signal
2395 : const AirState &propN, // Node 1 properties
2396 : const AirState &propM, // Node 2 properties
2397 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2398 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2399 : )
2400 : {
2401 :
2402 : // SUBROUTINE INFORMATION:
2403 : // AUTHOR George Walton
2404 : // DATE WRITTEN Extracted from AIRNET
2405 : // MODIFIED Lixing Gu, 2/1/04
2406 : // Revised the subroutine to meet E+ needs
2407 : // MODIFIED Lixing Gu, 6/8/05
2408 : // RE-ENGINEERED na
2409 :
2410 : // PURPOSE OF THIS SUBROUTINE:
2411 : // This subroutine solves airflow for a terminal unit component
2412 :
2413 : // SUBROUTINE PARAMETER DEFINITIONS:
2414 224 : Real64 constexpr C(0.868589);
2415 224 : Real64 constexpr EPS(0.001);
2416 224 : Real64 constexpr Rough(0.0001);
2417 224 : Real64 constexpr InitLamCoef(128.0);
2418 224 : Real64 constexpr LamDynCoef(64.0);
2419 224 : Real64 constexpr LamFriCoef(0.0001);
2420 224 : Real64 constexpr TurDynCoef(0.0001);
2421 :
2422 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2423 : Real64 A0;
2424 : Real64 A1;
2425 : Real64 A2;
2426 : Real64 B;
2427 : Real64 D;
2428 : Real64 S2;
2429 : Real64 CDM;
2430 : Real64 FL;
2431 : Real64 FT;
2432 : Real64 FTT;
2433 : Real64 RE;
2434 : Real64 ed;
2435 : Real64 ld;
2436 : Real64 g;
2437 : Real64 AA1;
2438 : Real64 area;
2439 :
2440 : // Formats
2441 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
2442 :
2443 : // Get component properties
2444 224 : ed = Rough / hydraulicDiameter;
2445 224 : area = pow_2(hydraulicDiameter) * Constant::Pi;
2446 224 : ld = L / hydraulicDiameter;
2447 224 : g = 1.14 - 0.868589 * std::log(ed);
2448 224 : AA1 = g;
2449 :
2450 224 : if (LFLAG) {
2451 : // Initialization by linear relation.
2452 0 : if (PDROP >= 0.0) {
2453 0 : DF[0] = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * InitLamCoef * ld);
2454 : } else {
2455 0 : DF[0] = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * InitLamCoef * ld);
2456 : }
2457 0 : F[0] = -DF[0] * PDROP;
2458 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwi:" << i << InitLamCoef << F[0] << DF[0];
2459 : } else {
2460 : // Standard calculation.
2461 224 : if (PDROP >= 0.0) {
2462 : // Flow in positive direction.
2463 : // Laminar flow coefficient !=0
2464 : if (LamFriCoef >= 0.001) {
2465 : A2 = LamFriCoef / (2.0 * propN.density * area * area);
2466 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * area * hydraulicDiameter);
2467 : A0 = -PDROP;
2468 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2469 : FL = (CDM - A1) / (2.0 * A2);
2470 : CDM = 1.0 / CDM;
2471 : } else {
2472 224 : CDM = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
2473 224 : FL = CDM * PDROP;
2474 : }
2475 224 : RE = FL * hydraulicDiameter / (propN.viscosity * area);
2476 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
2477 : // Turbulent flow; test when Re>10.
2478 224 : if (RE >= 10.0) {
2479 202 : S2 = std::sqrt(2.0 * propN.density * PDROP) * area;
2480 202 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2481 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
2482 : while (true) {
2483 830 : FT = FTT;
2484 830 : B = (9.3 * propN.viscosity * area) / (FT * Rough);
2485 830 : D = 1.0 + g * B;
2486 830 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2487 830 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2488 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
2489 830 : if (std::abs(FTT - FT) / FTT < EPS) break;
2490 : }
2491 202 : FT = FTT;
2492 : } else {
2493 22 : FT = FL;
2494 : }
2495 : } else {
2496 : // Flow in negative direction.
2497 : // Laminar flow coefficient !=0
2498 : if (LamFriCoef >= 0.001) {
2499 : A2 = LamFriCoef / (2.0 * propM.density * area * area);
2500 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * area * hydraulicDiameter);
2501 : A0 = PDROP;
2502 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2503 : FL = -(CDM - A1) / (2.0 * A2);
2504 : CDM = 1.0 / CDM;
2505 : } else {
2506 0 : CDM = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
2507 0 : FL = CDM * PDROP;
2508 : }
2509 0 : RE = -FL * hydraulicDiameter / (propM.viscosity * area);
2510 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwl:" << i << PDROP << FL << CDM << RE;
2511 : // Turbulent flow; test when Re>10.
2512 0 : if (RE >= 10.0) {
2513 0 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * area;
2514 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2515 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << S2 << FTT << g;
2516 : while (true) {
2517 0 : FT = FTT;
2518 0 : B = (9.3 * propM.viscosity * area) / (FT * Rough);
2519 0 : D = 1.0 + g * B;
2520 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2521 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2522 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " dwt:" << i << B << FTT << g;
2523 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2524 : }
2525 0 : FT = -FTT;
2526 : } else {
2527 0 : FT = FL;
2528 : }
2529 : }
2530 : // Select laminar or turbulent flow.
2531 224 : if (std::abs(FL) <= std::abs(FT)) {
2532 44 : F[0] = FL;
2533 44 : DF[0] = CDM;
2534 : } else {
2535 180 : F[0] = FT;
2536 180 : DF[0] = 0.5 * FT / PDROP;
2537 : }
2538 : }
2539 : // If damper, setup the airflows from nodal values calculated from terminal
2540 224 : if (state.afn->AirflowNetworkLinkageData(i).VAVTermDamper) {
2541 0 : F[0] = state.dataLoopNodes->Node(DamperInletNode).MassFlowRate;
2542 0 : if (state.afn->VAVTerminalRatio > 0.0) {
2543 0 : F[0] *= state.afn->VAVTerminalRatio;
2544 : }
2545 0 : DF[0] = 0.0;
2546 : }
2547 224 : return 1;
2548 : }
2549 :
2550 0 : int DisSysCompHXProp::calculate([[maybe_unused]] EnergyPlusData &state,
2551 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
2552 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2553 : [[maybe_unused]] int const i, // Linkage number
2554 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2555 : [[maybe_unused]] const Real64 control, // Element control signal
2556 : const AirState &propN, // Node 1 properties
2557 : const AirState &propM, // Node 2 properties
2558 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2559 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2560 : )
2561 : {
2562 :
2563 : // SUBROUTINE INFORMATION:
2564 : // AUTHOR George Walton
2565 : // DATE WRITTEN Extracted from AIRNET
2566 : // MODIFIED Lixing Gu, 2/1/04
2567 : // Revised the subroutine to meet E+ needs
2568 : // MODIFIED Lixing Gu, 1/18/09
2569 : // RE-ENGINEERED na
2570 :
2571 : // PURPOSE OF THIS SUBROUTINE:
2572 : // This subroutine solves airflow for a heat exchanger component
2573 :
2574 : // SUBROUTINE PARAMETER DEFINITIONS:
2575 0 : Real64 constexpr C(0.868589);
2576 0 : Real64 constexpr EPS(0.001);
2577 0 : Real64 constexpr Rough(0.0001);
2578 0 : Real64 constexpr InitLamCoef(128.0);
2579 0 : Real64 constexpr LamDynCoef(64.0);
2580 0 : Real64 constexpr LamFriCoef(0.0001);
2581 0 : Real64 constexpr TurDynCoef(0.0001);
2582 :
2583 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2584 : Real64 A0;
2585 : Real64 A1;
2586 : Real64 A2;
2587 : Real64 B;
2588 : Real64 D;
2589 : Real64 S2;
2590 : Real64 CDM;
2591 : Real64 FL;
2592 : Real64 FT;
2593 : Real64 FTT;
2594 : Real64 RE;
2595 : Real64 ed;
2596 : Real64 ld;
2597 : Real64 g;
2598 : Real64 AA1;
2599 : Real64 area;
2600 :
2601 : // Get component properties
2602 0 : ed = Rough / hydraulicDiameter;
2603 0 : area = pow_2(hydraulicDiameter) * Constant::Pi;
2604 0 : ld = L / hydraulicDiameter;
2605 0 : g = 1.14 - 0.868589 * std::log(ed);
2606 0 : AA1 = g;
2607 :
2608 0 : if (LFLAG) {
2609 : // Initialization by linear relation.
2610 0 : if (PDROP >= 0.0) {
2611 0 : DF[0] = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * InitLamCoef * ld);
2612 : } else {
2613 0 : DF[0] = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * InitLamCoef * ld);
2614 : }
2615 0 : F[0] = -DF[0] * PDROP;
2616 : } else {
2617 : // Standard calculation.
2618 0 : if (PDROP >= 0.0) {
2619 : // Flow in positive direction.
2620 : // Laminar flow coefficient !=0
2621 : if (LamFriCoef >= 0.001) {
2622 : A2 = LamFriCoef / (2.0 * propN.density * area * area);
2623 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * area * hydraulicDiameter);
2624 : A0 = -PDROP;
2625 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2626 : FL = (CDM - A1) / (2.0 * A2);
2627 : CDM = 1.0 / CDM;
2628 : } else {
2629 0 : CDM = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
2630 0 : FL = CDM * PDROP;
2631 : }
2632 0 : RE = FL * hydraulicDiameter / (propN.viscosity * area);
2633 : // Turbulent flow; test when Re>10.
2634 0 : if (RE >= 10.0) {
2635 0 : S2 = std::sqrt(2.0 * propN.density * PDROP) * area;
2636 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2637 : while (true) {
2638 0 : FT = FTT;
2639 0 : B = (9.3 * propN.viscosity * area) / (FT * Rough);
2640 0 : D = 1.0 + g * B;
2641 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2642 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2643 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2644 : }
2645 0 : FT = FTT;
2646 : } else {
2647 0 : FT = FL;
2648 : }
2649 : } else {
2650 : // Flow in negative direction.
2651 : // Laminar flow coefficient !=0
2652 : if (LamFriCoef >= 0.001) {
2653 : A2 = LamFriCoef / (2.0 * propM.density * area * area);
2654 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * area * hydraulicDiameter);
2655 : A0 = PDROP;
2656 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2657 : FL = -(CDM - A1) / (2.0 * A2);
2658 : CDM = 1.0 / CDM;
2659 : } else {
2660 0 : CDM = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
2661 0 : FL = CDM * PDROP;
2662 : }
2663 0 : RE = -FL * hydraulicDiameter / (propM.viscosity * area);
2664 : // Turbulent flow; test when Re>10.
2665 0 : if (RE >= 10.0) {
2666 0 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * area;
2667 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2668 : while (true) {
2669 0 : FT = FTT;
2670 0 : B = (9.3 * propM.viscosity * area) / (FT * Rough);
2671 0 : D = 1.0 + g * B;
2672 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2673 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2674 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2675 : }
2676 0 : FT = -FTT;
2677 : } else {
2678 0 : FT = FL;
2679 : }
2680 : }
2681 : // Select laminar or turbulent flow.
2682 0 : if (std::abs(FL) <= std::abs(FT)) {
2683 0 : F[0] = FL;
2684 0 : DF[0] = CDM;
2685 : } else {
2686 0 : F[0] = FT;
2687 0 : DF[0] = 0.5 * FT / PDROP;
2688 : }
2689 : }
2690 0 : return 1;
2691 : }
2692 :
2693 0 : int DisSysCompHXProp::calculate([[maybe_unused]] EnergyPlusData &state,
2694 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2695 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2696 : [[maybe_unused]] const Real64 control, // Element control signal
2697 : const AirState &propN, // Node 1 properties
2698 : const AirState &propM, // Node 2 properties
2699 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2700 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2701 : )
2702 : {
2703 :
2704 : // SUBROUTINE INFORMATION:
2705 : // AUTHOR George Walton
2706 : // DATE WRITTEN Extracted from AIRNET
2707 : // MODIFIED Lixing Gu, 2/1/04
2708 : // Revised the subroutine to meet E+ needs
2709 : // MODIFIED Lixing Gu, 1/18/09
2710 : // RE-ENGINEERED na
2711 :
2712 : // PURPOSE OF THIS SUBROUTINE:
2713 : // This subroutine solves airflow for a heat exchanger component
2714 :
2715 : // SUBROUTINE PARAMETER DEFINITIONS:
2716 0 : Real64 constexpr C(0.868589);
2717 0 : Real64 constexpr EPS(0.001);
2718 0 : Real64 constexpr Rough(0.0001);
2719 0 : Real64 constexpr LamDynCoef(64.0);
2720 0 : Real64 constexpr LamFriCoef(0.0001);
2721 0 : Real64 constexpr TurDynCoef(0.0001);
2722 :
2723 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2724 : Real64 A0;
2725 : Real64 A1;
2726 : Real64 A2;
2727 : Real64 B;
2728 : Real64 D;
2729 : Real64 S2;
2730 : Real64 CDM;
2731 : Real64 FL;
2732 : Real64 FT;
2733 : Real64 FTT;
2734 : Real64 RE;
2735 : Real64 ed;
2736 : Real64 ld;
2737 : Real64 g;
2738 : Real64 AA1;
2739 : Real64 area;
2740 :
2741 : // Get component properties
2742 0 : ed = Rough / hydraulicDiameter;
2743 0 : area = pow_2(hydraulicDiameter) * Constant::Pi;
2744 0 : ld = L / hydraulicDiameter;
2745 0 : g = 1.14 - 0.868589 * std::log(ed);
2746 0 : AA1 = g;
2747 :
2748 : // Standard calculation.
2749 0 : if (PDROP >= 0.0) {
2750 : // Flow in positive direction.
2751 : // Laminar flow coefficient !=0
2752 : if (LamFriCoef >= 0.001) {
2753 : A2 = LamFriCoef / (2.0 * propN.density * area * area);
2754 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * area * hydraulicDiameter);
2755 : A0 = -PDROP;
2756 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2757 : FL = (CDM - A1) / (2.0 * A2);
2758 : CDM = 1.0 / CDM;
2759 : } else {
2760 0 : CDM = (2.0 * propN.density * area * hydraulicDiameter) / (propN.viscosity * LamDynCoef * ld);
2761 0 : FL = CDM * PDROP;
2762 : }
2763 0 : RE = FL * hydraulicDiameter / (propN.viscosity * area);
2764 : // Turbulent flow; test when Re>10.
2765 0 : if (RE >= 10.0) {
2766 0 : S2 = std::sqrt(2.0 * propN.density * PDROP) * area;
2767 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2768 : while (true) {
2769 0 : FT = FTT;
2770 0 : B = (9.3 * propN.viscosity * area) / (FT * Rough);
2771 0 : D = 1.0 + g * B;
2772 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2773 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2774 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2775 : }
2776 0 : FT = FTT;
2777 : } else {
2778 0 : FT = FL;
2779 : }
2780 : } else {
2781 : // Flow in negative direction.
2782 : // Laminar flow coefficient !=0
2783 : if (LamFriCoef >= 0.001) {
2784 : A2 = LamFriCoef / (2.0 * propM.density * area * area);
2785 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * area * hydraulicDiameter);
2786 : A0 = PDROP;
2787 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
2788 : FL = -(CDM - A1) / (2.0 * A2);
2789 : CDM = 1.0 / CDM;
2790 : } else {
2791 0 : CDM = (2.0 * propM.density * area * hydraulicDiameter) / (propM.viscosity * LamDynCoef * ld);
2792 0 : FL = CDM * PDROP;
2793 : }
2794 0 : RE = -FL * hydraulicDiameter / (propM.viscosity * area);
2795 : // Turbulent flow; test when Re>10.
2796 0 : if (RE >= 10.0) {
2797 0 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * area;
2798 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2799 : while (true) {
2800 0 : FT = FTT;
2801 0 : B = (9.3 * propM.viscosity * area) / (FT * Rough);
2802 0 : D = 1.0 + g * B;
2803 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
2804 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
2805 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
2806 : }
2807 0 : FT = -FTT;
2808 : } else {
2809 0 : FT = FL;
2810 : }
2811 : }
2812 : // Select laminar or turbulent flow.
2813 0 : if (std::abs(FL) <= std::abs(FT)) {
2814 0 : F[0] = FL;
2815 0 : DF[0] = CDM;
2816 : } else {
2817 0 : F[0] = FT;
2818 0 : DF[0] = 0.5 * FT / PDROP;
2819 : }
2820 0 : return 1;
2821 : }
2822 :
2823 196615 : int ZoneExhaustFan::calculate(EnergyPlusData &state,
2824 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
2825 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2826 : int const i, // Linkage number
2827 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2828 : [[maybe_unused]] const Real64 control, // Element control signal
2829 : const AirState &propN, // Node 1 properties
2830 : const AirState &propM, // Node 2 properties
2831 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2832 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2833 : )
2834 : {
2835 : // SUBROUTINE INFORMATION:
2836 : // AUTHOR George Walton
2837 : // DATE WRITTEN Extracted from AIRNET
2838 : // MODIFIED Lixing Gu, 12/17/06
2839 : // Revised for zone exhaust fan
2840 : // RE-ENGINEERED na
2841 :
2842 : // PURPOSE OF THIS SUBROUTINE:
2843 : // This subroutine solves airflow for a surface crack component
2844 :
2845 : // Using/Aliasing
2846 : using HVAC::VerySmallMassFlow;
2847 :
2848 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2849 : Real64 CDM;
2850 : Real64 FL;
2851 : Real64 FT;
2852 : Real64 RhozNorm;
2853 : Real64 VisczNorm;
2854 : Real64 expn;
2855 : Real64 Ctl;
2856 : Real64 coef;
2857 : Real64 Corr;
2858 : Real64 VisAve;
2859 : Real64 Tave;
2860 : Real64 RhoCor;
2861 : // int InletNode;
2862 :
2863 : // Formats
2864 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
2865 :
2866 196615 : if (state.dataLoopNodes->Node(InletNode).MassFlowRate > VerySmallMassFlow) {
2867 : // Treat the component as an exhaust fan
2868 180191 : if (state.afn->PressureSetFlag == PressureCtrlExhaust) {
2869 0 : F[0] = state.afn->ExhaustFanMassFlowRate;
2870 : } else {
2871 180191 : F[0] = state.dataLoopNodes->Node(InletNode).MassFlowRate;
2872 : }
2873 180191 : DF[0] = 0.0;
2874 180191 : return 1;
2875 : } else {
2876 : // Treat the component as a surface crack
2877 : // Crack standard condition from given inputs
2878 16424 : Corr = state.afn->MultizoneSurfaceData(i).Factor;
2879 16424 : RhozNorm = state.afn->properties.density(StandardP, StandardT, StandardW);
2880 16424 : VisczNorm = 1.71432e-5 + 4.828e-8 * StandardT;
2881 :
2882 16424 : expn = FlowExpo;
2883 16424 : VisAve = (propN.viscosity + propM.viscosity) / 2.0;
2884 16424 : Tave = (propN.temperature + propM.temperature) / 2.0;
2885 16424 : if (PDROP >= 0.0) {
2886 12837 : coef = FlowCoef / propN.sqrt_density * Corr;
2887 : } else {
2888 3587 : coef = FlowCoef / propM.sqrt_density * Corr;
2889 : }
2890 :
2891 16424 : if (LFLAG) {
2892 : // Initialization by linear relation.
2893 0 : if (PDROP >= 0.0) {
2894 0 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
2895 0 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
2896 0 : DF[0] = coef * propN.density / propN.viscosity * Ctl;
2897 : } else {
2898 0 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
2899 0 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
2900 0 : DF[0] = coef * propM.density / propM.viscosity * Ctl;
2901 : }
2902 0 : F[0] = -DF[0] * PDROP;
2903 : } else {
2904 : // Standard calculation.
2905 16424 : if (PDROP >= 0.0) {
2906 : // Flow in positive direction.
2907 : // Laminar flow.
2908 12837 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
2909 12837 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
2910 12837 : CDM = coef * propN.density / propN.viscosity * Ctl;
2911 12837 : FL = CDM * PDROP;
2912 : // Turbulent flow.
2913 12837 : if (expn == 0.5) {
2914 0 : FT = coef * propN.sqrt_density * std::sqrt(PDROP) * Ctl;
2915 : } else {
2916 12837 : FT = coef * propN.sqrt_density * std::pow(PDROP, expn) * Ctl;
2917 : }
2918 : } else {
2919 : // Flow in negative direction.
2920 : // Laminar flow.
2921 3587 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
2922 3587 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
2923 3587 : CDM = coef * propM.density / propM.viscosity * Ctl;
2924 3587 : FL = CDM * PDROP;
2925 : // Turbulent flow.
2926 3587 : if (expn == 0.5) {
2927 0 : FT = -coef * propM.sqrt_density * std::sqrt(-PDROP) * Ctl;
2928 : } else {
2929 3587 : FT = -coef * propM.sqrt_density * std::pow(-PDROP, expn) * Ctl;
2930 : }
2931 : }
2932 : // Select laminar or turbulent flow.
2933 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " scr: " << i << PDROP << FL << FT;
2934 16424 : if (std::abs(FL) <= std::abs(FT)) {
2935 0 : F[0] = FL;
2936 0 : DF[0] = CDM;
2937 : } else {
2938 16424 : F[0] = FT;
2939 16424 : DF[0] = FT * expn / PDROP;
2940 : }
2941 : }
2942 : }
2943 16424 : return 1;
2944 : }
2945 :
2946 0 : int ZoneExhaustFan::calculate(EnergyPlusData &state,
2947 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
2948 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
2949 : const Real64 control, // Element control signal
2950 : const AirState &propN, // Node 1 properties
2951 : const AirState &propM, // Node 2 properties
2952 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
2953 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
2954 : )
2955 : {
2956 : // SUBROUTINE INFORMATION:
2957 : // AUTHOR George Walton
2958 : // DATE WRITTEN Extracted from AIRNET
2959 : // MODIFIED Lixing Gu, 12/17/06
2960 : // Revised for zone exhaust fan
2961 : // RE-ENGINEERED na
2962 :
2963 : // PURPOSE OF THIS SUBROUTINE:
2964 : // This subroutine solves airflow for a surface crack component
2965 :
2966 : // Using/Aliasing
2967 : using HVAC::VerySmallMassFlow;
2968 :
2969 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2970 : Real64 CDM;
2971 : Real64 FL;
2972 : Real64 FT;
2973 : Real64 RhozNorm;
2974 : Real64 VisczNorm;
2975 : Real64 expn;
2976 : Real64 Ctl;
2977 : Real64 coef;
2978 : Real64 VisAve;
2979 : Real64 Tave;
2980 : Real64 RhoCor;
2981 : // int InletNode;
2982 :
2983 : // Formats
2984 : // static gio::Fmt Format_901("(A5,I3,6X,4E16.7)");
2985 :
2986 0 : if (state.dataLoopNodes->Node(InletNode).MassFlowRate > VerySmallMassFlow) {
2987 : // Treat the component as an exhaust fan
2988 0 : if (state.afn->PressureSetFlag == PressureCtrlExhaust) {
2989 0 : F[0] = state.afn->ExhaustFanMassFlowRate;
2990 : } else {
2991 0 : F[0] = state.dataLoopNodes->Node(InletNode).MassFlowRate;
2992 : }
2993 0 : DF[0] = 0.0;
2994 0 : return 1;
2995 : } else {
2996 : // Treat the component as a surface crack
2997 : // Crack standard condition from given inputs
2998 0 : RhozNorm = state.afn->properties.density(StandardP, StandardT, StandardW);
2999 0 : VisczNorm = 1.71432e-5 + 4.828e-8 * StandardT;
3000 :
3001 0 : expn = FlowExpo;
3002 0 : VisAve = (propN.viscosity + propM.viscosity) / 2.0;
3003 0 : Tave = (propN.temperature + propM.temperature) / 2.0;
3004 0 : if (PDROP >= 0.0) {
3005 0 : coef = control * FlowCoef / propN.sqrt_density;
3006 : } else {
3007 0 : coef = control * FlowCoef / propM.sqrt_density;
3008 : }
3009 :
3010 : // Standard calculation.
3011 0 : if (PDROP >= 0.0) {
3012 : // Flow in positive direction.
3013 : // Laminar flow.
3014 0 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
3015 0 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3016 0 : CDM = coef * propN.density / propN.viscosity * Ctl;
3017 0 : FL = CDM * PDROP;
3018 : // Turbulent flow.
3019 0 : if (expn == 0.5) {
3020 0 : FT = coef * propN.sqrt_density * std::sqrt(PDROP) * Ctl;
3021 : } else {
3022 0 : FT = coef * propN.sqrt_density * std::pow(PDROP, expn) * Ctl;
3023 : }
3024 : } else {
3025 : // Flow in negative direction.
3026 : // Laminar flow.
3027 0 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
3028 0 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3029 0 : CDM = coef * propM.density / propM.viscosity * Ctl;
3030 0 : FL = CDM * PDROP;
3031 : // Turbulent flow.
3032 0 : if (expn == 0.5) {
3033 0 : FT = -coef * propM.sqrt_density * std::sqrt(-PDROP) * Ctl;
3034 : } else {
3035 0 : FT = -coef * propM.sqrt_density * std::pow(-PDROP, expn) * Ctl;
3036 : }
3037 : }
3038 : // Select laminar or turbulent flow.
3039 : // if (LIST >= 4) gio::write(Unit21, Format_901) << " scr: " << i << PDROP << FL << FT;
3040 0 : if (std::abs(FL) <= std::abs(FT)) {
3041 0 : F[0] = FL;
3042 0 : DF[0] = CDM;
3043 : } else {
3044 0 : F[0] = FT;
3045 0 : DF[0] = FT * expn / PDROP;
3046 : }
3047 : }
3048 0 : return 1;
3049 : }
3050 :
3051 2 : int HorizontalOpening::calculate(EnergyPlusData &state,
3052 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
3053 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
3054 : int const i, // Linkage number
3055 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
3056 : [[maybe_unused]] const Real64 control, // Element control signal
3057 : const AirState &propN, // Node 1 properties
3058 : const AirState &propM, // Node 2 properties
3059 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3060 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3061 : )
3062 : {
3063 : // SUBROUTINE INFORMATION:
3064 : // AUTHOR Lixing Gu
3065 : // DATE WRITTEN Apr. 2009
3066 : // MODIFIED na
3067 : // MODIFIED na
3068 : // RE-ENGINEERED na
3069 :
3070 : // PURPOSE OF THIS SUBROUTINE:
3071 : // This subroutine solves airflow for a horizontal opening component. The subroutine was
3072 : // developed based on the subroutine AFEPLR of AIRNET.
3073 :
3074 : // METHODOLOGY EMPLOYED:
3075 : // Combine forced and buyancy airflows together with a cap
3076 :
3077 : // REFERENCES:
3078 : // Cooper, L., 1989, "Calculation of the Flow Through a Horizontal Ceiling/Floor Vent,"
3079 : // NISTIR 89-4052, National Institute of Standards and Technology, Gaithersburg, MD
3080 :
3081 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3082 : Real64 RhozAver;
3083 : Real64 expn;
3084 : Real64 coef;
3085 : Real64 Width; // Opening width
3086 : Real64 Height; // Opening height
3087 : Real64 Fact; // Opening factor
3088 : // Real64 Slope; // Opening slope
3089 : // Real64 DischCoeff; // Discharge coefficient
3090 : Real64 fma12; // massflow in direction "from-to" [kg/s]
3091 : Real64 fma21; // massflow in direction "to-from" [kg/s]
3092 : Real64 dp1fma12; // derivative d fma12 / d Dp [kg/s/Pa]
3093 : Real64 dp1fma21; // derivative d fma21 / d Dp [kg/s/Pa]
3094 : Real64 PurgedP; // Purge pressure [Pa]
3095 : Real64 BuoFlow; // Buoyancy flow rate [Pa]
3096 : Real64 BuoFlowMax; // Maximum buoyancy flow rate [Pa]
3097 : Real64 dPBuoFlow; // Derivative of buoyancy flow rate [kg/s/Pa]
3098 : Real64 DH; // Hydraulic diameter [m]
3099 : Real64 Cshape; // Shape factor [dimensionless]
3100 : Real64 OpenArea; // Opening area [m2]
3101 :
3102 : // Get information on the horizontal opening
3103 2 : RhozAver = (propN.density + propM.density) / 2.0;
3104 2 : Width = state.afn->MultizoneSurfaceData(i).Width;
3105 2 : Height = state.afn->MultizoneSurfaceData(i).Height;
3106 2 : Fact = state.afn->MultizoneSurfaceData(i).OpenFactor;
3107 2 : expn = FlowExpo;
3108 2 : coef = FlowCoef;
3109 : // Slope = MultizoneCompHorOpeningData(CompNum).Slope;
3110 : // DischCoeff = MultizoneCompHorOpeningData(CompNum).DischCoeff;
3111 2 : Cshape = 0.942 * Width / Height;
3112 2 : OpenArea = Width * Height * Fact * std::sin(Slope * Constant::Pi / 180.0) * (1.0 + std::cos(Slope * Constant::Pi / 180.0));
3113 2 : DH = 4.0 * (Width * Height) / 2.0 / (Width + Height) * Fact;
3114 :
3115 : // Check which zone is higher
3116 2 : if (Fact == 0.0) {
3117 0 : generic_crack(coef, expn, LFLAG, PDROP, propN, propM, F, DF);
3118 0 : return 1;
3119 : }
3120 :
3121 2 : fma12 = 0.0;
3122 2 : fma21 = 0.0;
3123 2 : dp1fma12 = 0.0;
3124 2 : dp1fma21 = 0.0;
3125 2 : BuoFlow = 0.0;
3126 2 : dPBuoFlow = 0.0;
3127 :
3128 2 : if (state.afn->AirflowNetworkLinkageData(i).NodeHeights[0] > state.afn->AirflowNetworkLinkageData(i).NodeHeights[1]) {
3129 : // Node N is upper zone
3130 2 : if (propN.density > propM.density) {
3131 2 : BuoFlowMax = RhozAver * 0.055 * std::sqrt(9.81 * std::abs(propN.density - propM.density) * pow_5(DH) / RhozAver);
3132 2 : PurgedP = Cshape * Cshape * 9.81 * std::abs(propN.density - propM.density) * pow_5(DH) / (2.0 * pow_2(OpenArea));
3133 2 : if (std::abs(PDROP) <= PurgedP) {
3134 2 : BuoFlow = BuoFlowMax * (1.0 - std::abs(PDROP) / PurgedP);
3135 2 : dPBuoFlow = BuoFlowMax / PurgedP;
3136 : }
3137 : }
3138 : } else {
3139 : // Node M is upper zone
3140 0 : if (propN.density < propM.density) {
3141 0 : BuoFlowMax = RhozAver * 0.055 * std::sqrt(9.81 * std::abs(propN.density - propM.density) * pow_5(DH) / RhozAver);
3142 0 : PurgedP = Cshape * Cshape * 9.81 * std::abs(propN.density - propM.density) * pow_5(DH) / (2.0 * pow_2(OpenArea));
3143 0 : if (std::abs(PDROP) <= PurgedP) {
3144 0 : BuoFlow = BuoFlowMax * (1.0 - std::abs(PDROP) / PurgedP);
3145 0 : dPBuoFlow = BuoFlowMax / PurgedP;
3146 : }
3147 : }
3148 : }
3149 :
3150 2 : if (PDROP == 0.0) {
3151 0 : fma12 = BuoFlow;
3152 0 : fma21 = BuoFlow;
3153 0 : dp1fma12 = 0.0;
3154 0 : dp1fma21 = 0.0;
3155 2 : } else if (PDROP > 0.0) {
3156 1 : fma12 = propN.density * OpenArea * Fact * DischCoeff * std::sqrt(2.0 * PDROP / RhozAver) + BuoFlow;
3157 1 : dp1fma12 = propN.density * OpenArea * DischCoeff / std::sqrt(2.0 * PDROP * RhozAver) + dPBuoFlow;
3158 1 : if (BuoFlow > 0.0) {
3159 1 : fma21 = BuoFlow;
3160 1 : dp1fma21 = dPBuoFlow;
3161 : }
3162 : } else { // PDROP.LT.0.0
3163 1 : fma21 = propM.density * OpenArea * Fact * DischCoeff * std::sqrt(2.0 * std::abs(PDROP) / RhozAver) + BuoFlow;
3164 1 : dp1fma21 = -propM.density * OpenArea * DischCoeff / std::sqrt(2.0 * std::abs(PDROP) * RhozAver) + dPBuoFlow;
3165 1 : if (BuoFlow > 0.0) {
3166 1 : fma12 = BuoFlow;
3167 1 : dp1fma12 = dPBuoFlow;
3168 : }
3169 : }
3170 :
3171 2 : F[0] = fma12 - fma21;
3172 2 : DF[0] = dp1fma12 - dp1fma21;
3173 2 : F[1] = 0.0;
3174 2 : if (fma12 != 0.0 && fma21 != 0.0) {
3175 2 : F[1] = BuoFlow;
3176 : }
3177 2 : DF[1] = 0.0;
3178 :
3179 2 : return 1;
3180 : }
3181 :
3182 4 : int SpecifiedMassFlow::calculate([[maybe_unused]] EnergyPlusData &state,
3183 : [[maybe_unused]] bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
3184 : [[maybe_unused]] Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
3185 : [[maybe_unused]] int const i, // Linkage number
3186 : const Real64 multiplier, // Element multiplier
3187 : const Real64 control, // Element control signal
3188 : [[maybe_unused]] const AirState &propN, // Node 1 properties
3189 : [[maybe_unused]] const AirState &propM, // Node 2 properties
3190 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3191 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3192 : )
3193 : {
3194 : // SUBROUTINE INFORMATION:
3195 : // AUTHOR Jason DeGraw and Prateek Shrestha
3196 : // DATE WRITTEN June 2021
3197 : // MODIFIED na
3198 : // MODIFIED na
3199 : // RE-ENGINEERED na
3200 :
3201 : // PURPOSE OF THIS SUBROUTINE:
3202 : // This subroutine computes airflow for a specified mass flow element.
3203 :
3204 4 : F[0] = mass_flow * control * multiplier;
3205 4 : DF[0] = 0.0;
3206 4 : F[1] = 0.0;
3207 4 : DF[1] = 0.0;
3208 :
3209 4 : return 1;
3210 : }
3211 :
3212 4 : int SpecifiedVolumeFlow::calculate([[maybe_unused]] EnergyPlusData &state,
3213 : [[maybe_unused]] bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
3214 : [[maybe_unused]] Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
3215 : [[maybe_unused]] int const i, // Linkage number
3216 : const Real64 multiplier, // Element multiplier
3217 : const Real64 control, // Element control signal
3218 : const AirState &propN, // Node 1 properties
3219 : const AirState &propM, // Node 2 properties
3220 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3221 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3222 : )
3223 : {
3224 : // SUBROUTINE INFORMATION:
3225 : // AUTHOR Jason DeGraw and Prateek Shrestha
3226 : // DATE WRITTEN June 2021
3227 : // MODIFIED na
3228 : // MODIFIED na
3229 : // RE-ENGINEERED na
3230 :
3231 : // PURPOSE OF THIS SUBROUTINE:
3232 : // This subroutine computes airflow for a specified volume flow element.
3233 :
3234 4 : Real64 flow = volume_flow * control * multiplier;
3235 :
3236 4 : Real64 upwind_density{propN.density};
3237 :
3238 4 : if (flow < 0.0) {
3239 0 : upwind_density = propM.density;
3240 : }
3241 :
3242 4 : F[0] = flow * upwind_density;
3243 4 : DF[0] = 0.0;
3244 4 : F[1] = 0.0;
3245 4 : DF[1] = 0.0;
3246 :
3247 4 : return 1;
3248 : }
3249 :
3250 172 : int OutdoorAirFan::calculate(EnergyPlusData &state,
3251 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
3252 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
3253 : int const i, // Linkage number
3254 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
3255 : [[maybe_unused]] const Real64 control, // Element control signal
3256 : const AirState &propN, // Node 1 properties
3257 : const AirState &propM, // Node 2 properties
3258 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3259 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3260 : )
3261 : {
3262 :
3263 : // PURPOSE OF THIS SUBROUTINE:
3264 : // This subroutine solves airflow for a constant flow rate airflow component -- using standard interface.
3265 :
3266 : // Using/Aliasing
3267 : using HVAC::VerySmallMassFlow;
3268 :
3269 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3270 : Real64 expn;
3271 : Real64 Ctl;
3272 : Real64 coef;
3273 : Real64 Corr;
3274 : Real64 VisAve;
3275 : Real64 Tave;
3276 : Real64 RhoCor;
3277 : // int InletNode;
3278 : Real64 RhozNorm;
3279 : Real64 VisczNorm;
3280 : Real64 CDM;
3281 : Real64 FL;
3282 : Real64 FT;
3283 :
3284 172 : int AirLoopNum = state.afn->AirflowNetworkLinkageData(i).AirLoopNum;
3285 :
3286 172 : if (state.dataLoopNodes->Node(InletNode).MassFlowRate > VerySmallMassFlow) {
3287 : // Treat the component as an exhaust fan
3288 172 : F[0] = state.dataLoopNodes->Node(InletNode).MassFlowRate;
3289 172 : DF[0] = 0.0;
3290 212 : if (state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopFanOperationMode == HVAC::FanOp::Cycling &&
3291 40 : state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopOnOffFanPartLoadRatio > 0.0) {
3292 40 : F[0] = F[0] / state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopOnOffFanPartLoadRatio;
3293 : }
3294 172 : return 1;
3295 : } else {
3296 : // Treat the component as a surface crack
3297 : // Crack standard condition from given inputs
3298 0 : Corr = 1.0;
3299 0 : RhozNorm = state.afn->properties.density(StandardP, StandardT, StandardW);
3300 0 : VisczNorm = 1.71432e-5 + 4.828e-8 * StandardT;
3301 :
3302 0 : expn = FlowExpo;
3303 0 : VisAve = (propN.viscosity + propM.viscosity) / 2.0;
3304 0 : Tave = (propN.temperature + propM.temperature) / 2.0;
3305 0 : if (PDROP >= 0.0) {
3306 0 : coef = FlowCoef / propN.sqrt_density * Corr;
3307 : } else {
3308 0 : coef = FlowCoef / propM.sqrt_density * Corr;
3309 : }
3310 :
3311 0 : if (LFLAG) {
3312 : // Initialization by linear relation.
3313 0 : if (PDROP >= 0.0) {
3314 0 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
3315 0 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3316 0 : DF[0] = coef * propN.density / propN.viscosity * Ctl;
3317 : } else {
3318 0 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
3319 0 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3320 0 : DF[0] = coef * propM.density / propM.viscosity * Ctl;
3321 : }
3322 0 : F[0] = -DF[0] * PDROP;
3323 : } else {
3324 : // Standard calculation.
3325 0 : if (PDROP >= 0.0) {
3326 : // Flow in positive direction.
3327 : // Laminar flow.
3328 0 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
3329 0 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3330 0 : CDM = coef * propN.density / propN.viscosity * Ctl;
3331 0 : FL = CDM * PDROP;
3332 : // Turbulent flow.
3333 0 : if (expn == 0.5) {
3334 0 : FT = coef * propN.sqrt_density * std::sqrt(PDROP) * Ctl;
3335 : } else {
3336 0 : FT = coef * propN.sqrt_density * std::pow(PDROP, expn) * Ctl;
3337 : }
3338 : } else {
3339 : // Flow in negative direction.
3340 : // Laminar flow.
3341 0 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
3342 0 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3343 0 : CDM = coef * propM.density / propM.viscosity * Ctl;
3344 0 : FL = CDM * PDROP;
3345 : // Turbulent flow.
3346 0 : if (expn == 0.5) {
3347 0 : FT = -coef * propM.sqrt_density * std::sqrt(-PDROP) * Ctl;
3348 : } else {
3349 0 : FT = -coef * propM.sqrt_density * std::pow(-PDROP, expn) * Ctl;
3350 : }
3351 : }
3352 : // Select laminar or turbulent flow.
3353 0 : if (std::abs(FL) <= std::abs(FT)) {
3354 0 : F[0] = FL;
3355 0 : DF[0] = CDM;
3356 : } else {
3357 0 : F[0] = FT;
3358 0 : DF[0] = FT * expn / PDROP;
3359 : }
3360 : }
3361 : }
3362 0 : return 1;
3363 : }
3364 :
3365 132 : int ReliefFlow::calculate(EnergyPlusData &state,
3366 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
3367 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
3368 : int const i, // Linkage number
3369 : [[maybe_unused]] const Real64 multiplier, // Element multiplier
3370 : [[maybe_unused]] const Real64 control, // Element control signal
3371 : const AirState &propN, // Node 1 properties
3372 : const AirState &propM, // Node 2 properties
3373 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3374 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3375 : )
3376 : {
3377 :
3378 : // PURPOSE OF THIS SUBROUTINE:
3379 : // This subroutine solves airflow for a constant flow rate airflow component -- using standard interface.
3380 :
3381 : // Using/Aliasing
3382 : using HVAC::VerySmallMassFlow;
3383 :
3384 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3385 : Real64 expn;
3386 : Real64 Ctl;
3387 : Real64 coef;
3388 : Real64 Corr;
3389 : Real64 VisAve;
3390 : Real64 Tave;
3391 : Real64 RhoCor;
3392 : Real64 RhozNorm;
3393 : Real64 VisczNorm;
3394 : Real64 CDM;
3395 : Real64 FL;
3396 : Real64 FT;
3397 :
3398 132 : int AirLoopNum = state.afn->AirflowNetworkLinkageData(i).AirLoopNum;
3399 :
3400 132 : if (state.dataLoopNodes->Node(OutletNode).MassFlowRate > VerySmallMassFlow) {
3401 : // Treat the component as an exhaust fan
3402 132 : DF[0] = 0.0;
3403 132 : if (state.afn->PressureSetFlag == PressureCtrlRelief) {
3404 92 : F[0] = state.afn->ReliefMassFlowRate;
3405 : } else {
3406 40 : F[0] = state.dataLoopNodes->Node(OutletNode).MassFlowRate;
3407 40 : if (state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopFanOperationMode == HVAC::FanOp::Cycling &&
3408 0 : state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopOnOffFanPartLoadRatio > 0.0) {
3409 0 : F[0] = F[0] / state.dataAirLoop->AirLoopAFNInfo(AirLoopNum).LoopOnOffFanPartLoadRatio;
3410 : }
3411 : }
3412 132 : return 1;
3413 : } else {
3414 : // Treat the component as a surface crack
3415 : // Crack standard condition from given inputs
3416 0 : Corr = 1.0;
3417 0 : RhozNorm = state.afn->properties.density(StandardP, StandardT, StandardW);
3418 0 : VisczNorm = 1.71432e-5 + 4.828e-8 * StandardT;
3419 :
3420 0 : expn = FlowExpo;
3421 0 : VisAve = (propN.viscosity + propM.viscosity) / 2.0;
3422 0 : Tave = (propN.temperature + propM.temperature) / 2.0;
3423 0 : if (PDROP >= 0.0) {
3424 0 : coef = FlowCoef / propN.sqrt_density * Corr;
3425 : } else {
3426 0 : coef = FlowCoef / propM.sqrt_density * Corr;
3427 : }
3428 :
3429 0 : if (LFLAG) {
3430 : // Initialization by linear relation.
3431 0 : if (PDROP >= 0.0) {
3432 0 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
3433 0 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3434 0 : DF[0] = coef * propN.density / propN.viscosity * Ctl;
3435 : } else {
3436 0 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
3437 0 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3438 0 : DF[0] = coef * propM.density / propM.viscosity * Ctl;
3439 : }
3440 0 : F[0] = -DF[0] * PDROP;
3441 : } else {
3442 : // Standard calculation.
3443 0 : if (PDROP >= 0.0) {
3444 : // Flow in positive direction.
3445 : // Laminar flow.
3446 0 : RhoCor = TOKELVIN(propN.temperature) / TOKELVIN(Tave);
3447 0 : Ctl = std::pow(RhozNorm / propN.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3448 0 : CDM = coef * propN.density / propN.viscosity * Ctl;
3449 0 : FL = CDM * PDROP;
3450 : // Turbulent flow.
3451 0 : if (expn == 0.5) {
3452 0 : FT = coef * propN.sqrt_density * std::sqrt(PDROP) * Ctl;
3453 : } else {
3454 0 : FT = coef * propN.sqrt_density * std::pow(PDROP, expn) * Ctl;
3455 : }
3456 : } else {
3457 : // Flow in negative direction.
3458 : // Laminar flow.
3459 0 : RhoCor = TOKELVIN(propM.temperature) / TOKELVIN(Tave);
3460 0 : Ctl = std::pow(RhozNorm / propM.density / RhoCor, expn - 1.0) * std::pow(VisczNorm / VisAve, 2.0 * expn - 1.0);
3461 0 : CDM = coef * propM.density / propM.viscosity * Ctl;
3462 0 : FL = CDM * PDROP;
3463 : // Turbulent flow.
3464 0 : if (expn == 0.5) {
3465 0 : FT = -coef * propM.sqrt_density * std::sqrt(-PDROP) * Ctl;
3466 : } else {
3467 0 : FT = -coef * propM.sqrt_density * std::pow(-PDROP, expn) * Ctl;
3468 : }
3469 : }
3470 : // Select laminar or turbulent flow.
3471 0 : if (std::abs(FL) <= std::abs(FT)) {
3472 0 : F[0] = FL;
3473 0 : DF[0] = CDM;
3474 : } else {
3475 0 : F[0] = FT;
3476 0 : DF[0] = FT * expn / PDROP;
3477 : }
3478 : }
3479 : }
3480 0 : return 1;
3481 : }
3482 :
3483 4 : void generic_crack(Real64 const coefficient, // Flow coefficient
3484 : Real64 const exponent, // Flow exponent
3485 : bool const linear, // Initialization flag. If true, use linear relationship
3486 : Real64 const pdrop, // Total pressure drop across a component (P1 - P2) [Pa]
3487 : const AirState &propN, // Node 1 properties
3488 : const AirState &propM, // Node 2 properties
3489 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3490 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3491 : )
3492 : {
3493 : // SUBROUTINE INFORMATION:
3494 : // AUTHOR George Walton
3495 : // DATE WRITTEN Extracted from AIRNET
3496 : // MODIFIED Lixing Gu, 2/1/04
3497 : // Revised the subroutine to meet E+ needs
3498 : // MODIFIED Lixing Gu, 6/8/05
3499 : // RE-ENGINEERED This subroutine is revised from AFEPLR developed by George Walton, NIST
3500 : // Jason DeGraw
3501 :
3502 : // PURPOSE OF THIS SUBROUTINE:
3503 : // This subroutine solves airflow for a power law component
3504 :
3505 : // METHODOLOGY EMPLOYED:
3506 : // Using Q=C(dP)^n
3507 :
3508 : // FLOW:
3509 : // Calculate normal density and viscocity at reference conditions
3510 4 : constexpr Real64 reference_density = AIRDENSITY_CONSTEXPR(101325.0, 20.0, 0.0);
3511 4 : constexpr Real64 reference_viscosity = 1.71432e-5 + 4.828e-8 * 20.0;
3512 :
3513 4 : Real64 VisAve{0.5 * (propN.viscosity + propM.viscosity)};
3514 4 : Real64 Tave{0.5 * (propN.temperature + propM.temperature)};
3515 :
3516 4 : Real64 sign{1.0};
3517 4 : Real64 upwind_temperature{propN.temperature};
3518 4 : Real64 upwind_density{propN.density};
3519 4 : Real64 upwind_viscosity{propN.viscosity};
3520 4 : Real64 upwind_sqrt_density{propN.sqrt_density};
3521 4 : Real64 abs_pdrop = pdrop;
3522 :
3523 4 : if (pdrop < 0.0) {
3524 2 : sign = -1.0;
3525 2 : upwind_temperature = propM.temperature;
3526 2 : upwind_density = propM.density;
3527 2 : upwind_viscosity = propM.viscosity;
3528 2 : upwind_sqrt_density = propM.sqrt_density;
3529 2 : abs_pdrop = -pdrop;
3530 : }
3531 :
3532 4 : Real64 coef = coefficient / upwind_sqrt_density;
3533 :
3534 : // Laminar calculation
3535 4 : Real64 RhoCor{TOKELVIN(upwind_temperature) / TOKELVIN(Tave)};
3536 4 : Real64 Ctl{std::pow(reference_density / upwind_density / RhoCor, exponent - 1.0) *
3537 4 : std::pow(reference_viscosity / VisAve, 2.0 * exponent - 1.0)};
3538 4 : Real64 CDM{coef * upwind_density / upwind_viscosity * Ctl};
3539 4 : Real64 FL{CDM * pdrop};
3540 :
3541 4 : if (linear) {
3542 2 : DF[0] = CDM;
3543 2 : F[0] = FL;
3544 : } else {
3545 : // Turbulent flow.
3546 : Real64 abs_FT;
3547 2 : if (exponent == 0.5) {
3548 0 : abs_FT = coef * upwind_sqrt_density * std::sqrt(abs_pdrop) * Ctl;
3549 : } else {
3550 2 : abs_FT = coef * upwind_sqrt_density * std::pow(abs_pdrop, exponent) * Ctl;
3551 : }
3552 : // Select linear or nonlinear flow.
3553 2 : if (std::abs(FL) <= abs_FT) {
3554 0 : F[0] = FL;
3555 0 : DF[0] = CDM;
3556 : } else {
3557 2 : F[0] = sign * abs_FT;
3558 2 : DF[0] = F[0] * exponent / pdrop;
3559 : }
3560 : }
3561 4 : }
3562 :
3563 0 : int GenericDuct(Real64 const Length, // Duct length
3564 : Real64 const Diameter, // Duct diameter
3565 : bool const LFLAG, // Initialization flag.If = 1, use laminar relationship
3566 : Real64 const PDROP, // Total pressure drop across a component (P1 - P2) [Pa]
3567 : const AirState &propN, // Node 1 properties
3568 : const AirState &propM, // Node 2 properties
3569 : std::array<Real64, 2> &F, // Airflow through the component [kg/s]
3570 : std::array<Real64, 2> &DF // Partial derivative: DF/DP
3571 : )
3572 : {
3573 :
3574 : // This subroutine solve air flow as a duct if fan has zero flow rate
3575 :
3576 : // SUBROUTINE PARAMETER DEFINITIONS:
3577 0 : Real64 constexpr C(0.868589);
3578 0 : Real64 constexpr EPS(0.001);
3579 0 : Real64 constexpr Rough(0.0001);
3580 0 : Real64 constexpr InitLamCoef(128.0);
3581 0 : Real64 constexpr LamDynCoef(64.0);
3582 0 : Real64 constexpr LamFriCoef(0.0001);
3583 0 : Real64 constexpr TurDynCoef(0.0001);
3584 :
3585 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3586 : Real64 A0;
3587 : Real64 A1;
3588 : Real64 A2;
3589 : Real64 B;
3590 : Real64 D;
3591 : Real64 S2;
3592 : Real64 CDM;
3593 : Real64 FL;
3594 : Real64 FT;
3595 : Real64 FTT;
3596 : Real64 RE;
3597 :
3598 : // Get component properties
3599 0 : Real64 ed = Rough / Diameter;
3600 0 : Real64 area = Diameter * Diameter * Constant::Pi / 4.0;
3601 0 : Real64 ld = Length / Diameter;
3602 0 : Real64 g = 1.14 - 0.868589 * std::log(ed);
3603 0 : Real64 AA1 = g;
3604 :
3605 0 : if (LFLAG) {
3606 : // Initialization by linear relation.
3607 0 : if (PDROP >= 0.0) {
3608 0 : DF[0] = (2.0 * propN.density * area * Diameter) / (propN.viscosity * InitLamCoef * ld);
3609 : } else {
3610 0 : DF[0] = (2.0 * propM.density * area * Diameter) / (propM.viscosity * InitLamCoef * ld);
3611 : }
3612 0 : F[0] = -DF[0] * PDROP;
3613 : } else {
3614 : // Standard calculation.
3615 0 : if (PDROP >= 0.0) {
3616 : // Flow in positive direction.
3617 : // Laminar flow coefficient !=0
3618 : if (LamFriCoef >= 0.001) {
3619 : A2 = LamFriCoef / (2.0 * propN.density * area * area);
3620 : A1 = (propN.viscosity * LamDynCoef * ld) / (2.0 * propN.density * area * Diameter);
3621 : A0 = -PDROP;
3622 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
3623 : FL = (CDM - A1) / (2.0 * A2);
3624 : CDM = 1.0 / CDM;
3625 : } else {
3626 0 : CDM = (2.0 * propN.density * area * Diameter) / (propN.viscosity * LamDynCoef * ld);
3627 0 : FL = CDM * PDROP;
3628 : }
3629 0 : RE = FL * Diameter / (propN.viscosity * area);
3630 : // Turbulent flow; test when Re>10.
3631 0 : if (RE >= 10.0) {
3632 0 : S2 = std::sqrt(2.0 * propN.density * PDROP) * area;
3633 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
3634 : while (true) {
3635 0 : FT = FTT;
3636 0 : B = (9.3 * propN.viscosity * area) / (FT * Rough);
3637 0 : D = 1.0 + g * B;
3638 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
3639 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
3640 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
3641 : }
3642 0 : FT = FTT;
3643 : } else {
3644 0 : FT = FL;
3645 : }
3646 : } else {
3647 : // Flow in negative direction.
3648 : // Laminar flow coefficient !=0
3649 : if (LamFriCoef >= 0.001) {
3650 : A2 = LamFriCoef / (2.0 * propM.density * area * area);
3651 : A1 = (propM.viscosity * LamDynCoef * ld) / (2.0 * propM.density * area * Diameter);
3652 : A0 = PDROP;
3653 : CDM = std::sqrt(A1 * A1 - 4.0 * A2 * A0);
3654 : FL = -(CDM - A1) / (2.0 * A2);
3655 : CDM = 1.0 / CDM;
3656 : } else {
3657 0 : CDM = (2.0 * propM.density * area * Diameter) / (propM.viscosity * LamDynCoef * ld);
3658 0 : FL = CDM * PDROP;
3659 : }
3660 0 : RE = -FL * Diameter / (propM.viscosity * area);
3661 : // Turbulent flow; test when Re>10.
3662 0 : if (RE >= 10.0) {
3663 0 : S2 = std::sqrt(-2.0 * propM.density * PDROP) * area;
3664 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
3665 : while (true) {
3666 0 : FT = FTT;
3667 0 : B = (9.3 * propM.viscosity * area) / (FT * Rough);
3668 0 : D = 1.0 + g * B;
3669 0 : g -= (g - AA1 + C * std::log(D)) / (1.0 + C * B / D);
3670 0 : FTT = S2 / std::sqrt(ld / pow_2(g) + TurDynCoef);
3671 0 : if (std::abs(FTT - FT) / FTT < EPS) break;
3672 : }
3673 0 : FT = -FTT;
3674 : } else {
3675 0 : FT = FL;
3676 : }
3677 : }
3678 : // Select laminar or turbulent flow.
3679 0 : if (std::abs(FL) <= std::abs(FT)) {
3680 0 : F[0] = FL;
3681 0 : DF[0] = CDM;
3682 : } else {
3683 0 : F[0] = FT;
3684 0 : DF[0] = 0.5 * FT / PDROP;
3685 : }
3686 : }
3687 0 : return 1;
3688 : }
3689 :
3690 2 : void DetailedOpeningSolver::presprofile(EnergyPlusData &state,
3691 : int const il, // Linkage number
3692 : int const Pprof, // Opening number
3693 : Real64 const G, // gravitation field strength [N/kg]
3694 : const Array1D<Real64> &DpF, // Stack pressures at start heights of Layers
3695 : const Array1D<Real64> &DpT, // Stack pressures at start heights of Layers
3696 : const Array1D<Real64> &BetaF, // Density gradients in the FROM zone (starting at linkheight) [Kg/m3/m]
3697 : const Array1D<Real64> &BetaT, // Density gradients in the TO zone (starting at linkheight) [Kg/m3/m]
3698 : const Array1D<Real64> &RhoStF, // Density at the start heights of Layers in the FROM zone
3699 : const Array1D<Real64> &RhoStT, // Density at the start heights of Layers in the TO zone
3700 : int const From, // Number of FROM zone
3701 : int const To, // Number of To zone
3702 : Real64 const ActLh, // Actual height of opening [m]
3703 : Real64 const OwnHeightFactor // Cosine of deviation angle of the opening plane from the vertical direction
3704 : )
3705 : {
3706 :
3707 : // SUBROUTINE INFORMATION:
3708 : // AUTHOR Lixing Gu
3709 : // DATE WRITTEN Oct. 2005
3710 : // MODIFIED na
3711 : // RE-ENGINEERED This subroutine is revised based on PresProfile subroutine from COMIS
3712 :
3713 : // PURPOSE OF THIS SUBROUTINE:
3714 : // This subroutine calculates for a large opening profiles of stack pressure difference and
3715 : // densities in the zones linked by the a detailed opening cmponent.
3716 :
3717 : // METHODOLOGY EMPLOYED:
3718 : // The profiles are obtained in the following
3719 : // way: - the opening is divided into NrInt vertical intervals
3720 : // - the stack pressure difference and densities in From-
3721 : // and To-zone are calculated at the centre of each
3722 : // interval aswell as at the top and bottom of the LO
3723 : // - these values are stored in the (NrInt+2)-dimensional
3724 : // arrays DpProf, RhoProfF, RhoProfT.
3725 : // The calculation of stack pressure and density in the two zones
3726 : // is based on the arrays DpF/T, RhoStF/T, BetaF/T. These arrays
3727 : // are calculated in the COMIS routine Lclimb. They contain the
3728 : // values of stack pressure and density at the startheight of the
3729 : // opening and at startheights of all layers lying inside the
3730 : // opening, and the density gradients across the layers.
3731 : // The effective startheight zl(1/2) in the From/To zone and the
3732 : // effective length actLh of the LO take into account the
3733 : // startheightfactor, heightfactor and ownheightfactor. Thus for
3734 : // slanted windows the range of the profiles is the vertical
3735 : // projection of the actual opening.
3736 :
3737 : // REFERENCES:
3738 : // Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
3739 : // Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990
3740 :
3741 : // Argument array dimensioning
3742 2 : EP_SIZE_CHECK(DpF, 2);
3743 2 : EP_SIZE_CHECK(DpT, 2);
3744 2 : EP_SIZE_CHECK(BetaF, 2);
3745 2 : EP_SIZE_CHECK(BetaT, 2);
3746 2 : EP_SIZE_CHECK(RhoStF, 2);
3747 2 : EP_SIZE_CHECK(RhoStT, 2);
3748 :
3749 : // Locals
3750 : // SUBROUTINE ARGUMENT DEFINITIONS:
3751 : // in the FROM zone (starting at linkheight) [Pa]
3752 : // (starting at linkheight) [Kg/m3]
3753 : // (starting at linkheight) [Kg/m3]
3754 :
3755 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3756 2 : Array1D<Real64> zF(2); // Startheights of layers in FROM-, TO-zone
3757 2 : Array1D<Real64> zT(2);
3758 2 : Array1D<Real64> zStF(2); // Startheights of layers within the LO, starting with the actual startheight of the LO.
3759 2 : Array1D<Real64> zStT(2);
3760 : // The values in the arrays DpF, DpT, BetaF, BetaT, RhoStF, RhoStT are calculated at these heights.
3761 : Real64 hghtsFR;
3762 : Real64 hghtsTR;
3763 2 : Array1D<Real64> hghtsF(NrInt + 2); // Heights of evaluation points for pressure and density profiles
3764 2 : Array1D<Real64> hghtsT(NrInt + 2);
3765 : Real64 Interval; // Distance between two evaluation points
3766 : Real64 delzF; // Interval between actual evaluation point and startheight of actual layer in FROM-, TO-zone
3767 : Real64 delzT;
3768 : int AnzLayF; // Number of layers in FROM-, TO-zone
3769 : int AnzLayT;
3770 : int lF; // Actual index for DpF/T, BetaF/T, RhoStF/T, zStF/T
3771 : int lT;
3772 : int n;
3773 : int i;
3774 : int k;
3775 :
3776 : // Initialization
3777 2 : delzF = 0.0;
3778 2 : delzT = 0.0;
3779 2 : Interval = ActLh * OwnHeightFactor / NrInt;
3780 :
3781 42 : for (n = 1; n <= NrInt; ++n) {
3782 40 : hghtsF(n + 1) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[0] + Interval * (n - 0.5);
3783 40 : hghtsT(n + 1) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[1] + Interval * (n - 0.5);
3784 : }
3785 2 : hghtsF(1) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[0];
3786 2 : hghtsT(1) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[1];
3787 2 : hghtsF(NrInt + 2) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[0] + ActLh * OwnHeightFactor;
3788 2 : hghtsT(NrInt + 2) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[1] + ActLh * OwnHeightFactor;
3789 :
3790 2 : lF = 1;
3791 2 : lT = 1;
3792 2 : if (From == 0) {
3793 0 : AnzLayF = 1;
3794 : } else {
3795 2 : AnzLayF = 0;
3796 : }
3797 2 : if (To == 0) {
3798 2 : AnzLayT = 1;
3799 : } else {
3800 0 : AnzLayT = 0;
3801 : }
3802 :
3803 2 : if (AnzLayF > 0) {
3804 0 : for (n = 1; n <= AnzLayF; ++n) {
3805 0 : zF(n) = 0.0;
3806 0 : if (hghtsF(1) < 0.0) zF(n) = hghtsF(1);
3807 : }
3808 : }
3809 :
3810 2 : if (AnzLayT > 0) {
3811 4 : for (n = 1; n <= AnzLayT; ++n) {
3812 2 : zT(n) = 0.0;
3813 2 : if (hghtsT(1) < 0.0) zT(n) = hghtsT(1);
3814 : }
3815 : }
3816 :
3817 2 : zStF(1) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[0];
3818 2 : i = 2;
3819 2 : k = 1;
3820 :
3821 2 : while (k <= AnzLayF) {
3822 0 : if (zF(k) > zStF(1)) break;
3823 0 : ++k;
3824 : }
3825 :
3826 2 : while (k <= AnzLayF) {
3827 0 : if (zF(k) > hghtsF(NrInt)) break;
3828 0 : zStF(i) = zF(k); // Autodesk:BoundsViolation zStF(i) @ i>2 and zF(k) @ k>2
3829 0 : ++i;
3830 0 : ++k;
3831 : }
3832 :
3833 2 : zStF(i) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[0] + ActLh * OwnHeightFactor; // Autodesk:BoundsViolation zStF(i) @ i>2
3834 2 : zStT(1) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[1];
3835 2 : i = 2;
3836 2 : k = 1;
3837 :
3838 4 : while (k <= AnzLayT) {
3839 2 : if (zT(k) > zStT(1)) break;
3840 2 : ++k;
3841 : }
3842 :
3843 2 : while (k <= AnzLayT) {
3844 0 : if (zT(k) > hghtsT(NrInt)) break; // Autodesk:BoundsViolation zT(k) @ k>2
3845 0 : zStT(i) = zT(k); // Autodesk:BoundsViolation zStF(i) @ i>2 and zT(k) @ k>2
3846 0 : ++i;
3847 0 : ++k;
3848 : }
3849 :
3850 2 : zStT(i) = state.afn->AirflowNetworkLinkageData(il).NodeHeights[1] + ActLh * OwnHeightFactor; // Autodesk:BoundsViolation zStT(i) @ i>2
3851 :
3852 : // Calculation of DpProf, RhoProfF, RhoProfT
3853 46 : for (i = 1; i <= NrInt + 2; ++i) {
3854 44 : hghtsFR = hghtsF(i);
3855 44 : hghtsTR = hghtsT(i);
3856 :
3857 : while (true) {
3858 44 : if (hghtsFR > zStF(lF + 1)) {
3859 0 : if (lF > 2) break;
3860 0 : ++lF;
3861 : }
3862 44 : if (hghtsFR <= zStF(lF + 1)) break;
3863 : }
3864 :
3865 : while (true) {
3866 44 : if (hghtsTR > zStT(lT + 1)) {
3867 0 : ++lT;
3868 : }
3869 44 : if (hghtsTR <= zStT(lT + 1)) break;
3870 : }
3871 :
3872 44 : delzF = hghtsF(i) - zStF(lF);
3873 44 : delzT = hghtsT(i) - zStT(lT);
3874 :
3875 44 : RhoProfF(i + Pprof) = RhoStF(lF) + BetaF(lF) * delzF;
3876 44 : RhoProfT(i + Pprof) = RhoStT(lT) + BetaT(lT) * delzT;
3877 :
3878 88 : DpProf(i + Pprof) = DpF(lF) - DpT(lT) - G * (RhoStF(lF) * delzF + BetaF(lF) * pow_2(delzF) / 2.0) +
3879 44 : G * (RhoStT(lT) * delzT + BetaT(lT) * pow_2(delzT) / 2.0);
3880 : }
3881 2 : }
3882 :
3883 16873 : void DetailedOpeningSolver::pstack(EnergyPlusData &state, std::vector<AirflowNetwork::AirState> &props, Array1D<Real64> &pz)
3884 : {
3885 :
3886 : // SUBROUTINE INFORMATION:
3887 : // AUTHOR Lixing Gu
3888 : // DATE WRITTEN Oct. 2005
3889 : // MODIFIED na
3890 : // RE-ENGINEERED This subroutine is revised based on PresProfile subroutine from COMIS
3891 :
3892 : // PURPOSE OF THIS SUBROUTINE:
3893 : // This subroutine calculates the stack pressures for a link between two zones
3894 :
3895 : // REFERENCES:
3896 : // Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
3897 : // Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990
3898 :
3899 : // USE STATEMENTS:
3900 : // Locals
3901 : // SUBROUTINE ARGUMENT DEFINITIONS:
3902 : // na
3903 :
3904 : // SUBROUTINE PARAMETER DEFINITIONS:
3905 16873 : Real64 constexpr PSea(101325.0);
3906 :
3907 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3908 : // REAL(r64) RhoOut ! air density outside [kg/m3]
3909 : Real64 G; // gravity field strength [N/kg]
3910 : Real64 Pbz; // Pbarom at entrance level [Pa]
3911 16873 : Array2D<Real64> RhoDrL(state.afn->NumOfLinksMultiZone, 2); // dry air density on both sides of the link [kg/m3]
3912 : Real64 TempL1; // Temp in From and To zone at link level [C]
3913 : Real64 TempL2;
3914 : // REAL(r64) Tout ! outside temperature [C]
3915 : Real64 Xhl1; // Humidity in From and To zone at link level [kg/kg]
3916 : Real64 Xhl2;
3917 : // REAL(r64) Xhout ! outside humidity [kg/kg]
3918 16873 : Array1D<Real64> Hfl(state.afn->NumOfLinksMultiZone); // Own height factor for large (slanted) openings
3919 : int Nl; // number of links
3920 :
3921 16873 : Array1D<Real64> DpF(2);
3922 : Real64 DpP;
3923 16873 : Array1D<Real64> DpT(2);
3924 : Real64 H;
3925 16873 : Array1D<Real64> RhoStF(2);
3926 16873 : Array1D<Real64> RhoStT(2);
3927 : Real64 RhoDrDummi;
3928 16873 : Array1D<Real64> BetaStF(2);
3929 16873 : Array1D<Real64> BetaStT(2);
3930 : Real64 T;
3931 : Real64 X;
3932 16873 : Array1D<Real64> HSt(2);
3933 : Real64 TzFrom;
3934 : Real64 XhzFrom;
3935 : Real64 TzTo;
3936 : Real64 XhzTo;
3937 : Real64 ActLh;
3938 : Real64 ActLOwnh;
3939 : Real64 Pref;
3940 : Real64 PzFrom;
3941 : Real64 PzTo;
3942 16873 : Array1D<Real64> RhoLd(2);
3943 : Real64 RhoStd;
3944 : int Fromz;
3945 : int Toz;
3946 : iComponentTypeNum Ltyp;
3947 : int i;
3948 : int ll;
3949 : int Pprof;
3950 : int OpenNum;
3951 :
3952 : Real64 RhoREF;
3953 : Real64 CONV;
3954 :
3955 16873 : RhoREF = state.afn->properties.density(PSea, state.dataEnvrn->OutDryBulbTemp, state.dataEnvrn->OutHumRat);
3956 :
3957 16873 : CONV = state.dataEnvrn->Latitude * 2.0 * Constant::Pi / 360.0;
3958 16873 : G = 9.780373 * (1.0 + 0.0052891 * pow_2(std::sin(CONV)) - 0.0000059 * pow_2(std::sin(2.0 * CONV)));
3959 :
3960 16873 : Hfl = 1.0;
3961 16873 : Pbz = state.dataEnvrn->OutBaroPress;
3962 16873 : Nl = state.afn->NumOfLinksMultiZone;
3963 16873 : OpenNum = 0;
3964 16873 : RhoLd(1) = 1.2;
3965 16873 : RhoLd(2) = 1.2;
3966 16873 : RhoStd = 1.2;
3967 :
3968 337332 : for (i = 1; i <= Nl; ++i) {
3969 : // Check surface tilt
3970 320459 : if (i <= Nl - state.afn->NumOfLinksIntraZone) { // Revised by L.Gu, on 9 / 29 / 10
3971 320461 : if (state.afn->AirflowNetworkLinkageData(i).DetOpenNum > 0 &&
3972 2 : state.dataSurface->Surface(state.afn->MultizoneSurfaceData(i).SurfNum).Tilt < 90) {
3973 0 : Hfl(i) = state.dataSurface->Surface(state.afn->MultizoneSurfaceData(i).SurfNum).SinTilt;
3974 : }
3975 : }
3976 : // Initialisation
3977 320459 : int From = state.afn->AirflowNetworkLinkageData(i).NodeNums[0];
3978 320459 : int To = state.afn->AirflowNetworkLinkageData(i).NodeNums[1];
3979 320459 : if (state.afn->AirflowNetworkNodeData(From).EPlusZoneNum > 0 && state.afn->AirflowNetworkNodeData(To).EPlusZoneNum > 0) {
3980 33810 : ll = 0;
3981 286649 : } else if (state.afn->AirflowNetworkNodeData(From).EPlusZoneNum == 0 && state.afn->AirflowNetworkNodeData(To).EPlusZoneNum > 0) {
3982 0 : ll = 1;
3983 : } else {
3984 286649 : ll = 3;
3985 : }
3986 :
3987 320459 : Ltyp = state.afn->AirflowNetworkCompData(state.afn->AirflowNetworkLinkageData(i).CompNum).CompTypeNum;
3988 320459 : if (Ltyp == iComponentTypeNum::DOP) {
3989 2 : ActLh = state.afn->MultizoneSurfaceData(i).Height;
3990 2 : ActLOwnh = ActLh * 1.0;
3991 : } else {
3992 320457 : ActLh = 0.0;
3993 320457 : ActLOwnh = 0.0;
3994 : }
3995 :
3996 320459 : TempL1 = props[From].temperature;
3997 320459 : Xhl1 = props[From].humidity_ratio;
3998 320459 : TzFrom = props[From].temperature;
3999 320459 : XhzFrom = props[From].humidity_ratio;
4000 : // RhoL1 = props[From].density;
4001 320459 : if (ll == 0 || ll == 3) {
4002 320459 : PzFrom = pz(From);
4003 : } else {
4004 0 : PzFrom = 0.0;
4005 0 : From = 0;
4006 : }
4007 :
4008 320459 : int ilayptr = 0;
4009 320459 : if (From == 0) ilayptr = 1;
4010 320459 : if (ilayptr == 0) {
4011 320459 : Fromz = 0;
4012 : } else {
4013 0 : Fromz = From;
4014 : }
4015 :
4016 320459 : TempL2 = props[To].temperature;
4017 320459 : Xhl2 = props[To].humidity_ratio;
4018 320459 : TzTo = props[To].temperature;
4019 320459 : XhzTo = props[To].humidity_ratio;
4020 : // RhoL2 = props[To].density;
4021 :
4022 320459 : if (ll < 3) {
4023 33810 : PzTo = pz(To);
4024 : } else {
4025 286649 : PzTo = 0.0;
4026 286649 : To = 0;
4027 : }
4028 320459 : ilayptr = 0;
4029 320459 : if (To == 0) ilayptr = 1;
4030 320459 : if (ilayptr == 0) {
4031 33810 : Toz = 0;
4032 : } else {
4033 286649 : Toz = To;
4034 : }
4035 :
4036 : // RhoDrL is Rho at link level without pollutant but with humidity
4037 320459 : RhoDrL(i, 1) = state.afn->properties.density(state.dataEnvrn->OutBaroPress + PzFrom, TempL1, Xhl1);
4038 320459 : RhoDrL(i, 2) = state.afn->properties.density(state.dataEnvrn->OutBaroPress + PzTo, TempL2, Xhl2);
4039 :
4040 : // End initialisation
4041 :
4042 : // calculate DpF the difference between Pz and P at Node 1 height
4043 320459 : ilayptr = 0;
4044 320459 : if (Fromz == 0) ilayptr = 1;
4045 320459 : int j = ilayptr;
4046 320459 : int k = 1;
4047 320459 : lclimb(state, G, RhoLd(1), state.afn->AirflowNetworkLinkageData(i).NodeHeights[0], TempL1, Xhl1, DpF(k), Toz, PzTo, Pbz, RhoDrL(i, 1));
4048 320459 : Real64 RhoL1 = RhoLd(1);
4049 : // For large openings calculate the stack pressure difference profile and the
4050 : // density profile within the the top- and the bottom- height of the large opening
4051 320459 : if (ActLOwnh > 0.0) {
4052 2 : HSt(k) = state.afn->AirflowNetworkLinkageData(i).NodeHeights[0];
4053 2 : RhoStF(k) = RhoL1;
4054 2 : ++k;
4055 2 : HSt(k) = 0.0;
4056 2 : if (HSt(k - 1) < 0.0) HSt(k) = HSt(k - 1);
4057 :
4058 : // Search for the first startheight of a layer which is within the top- and the
4059 : // bottom- height of the large opening.
4060 : while (true) {
4061 4 : ilayptr = 0;
4062 4 : if (Fromz == 0) ilayptr = 9;
4063 4 : if ((j > ilayptr) || (HSt(k) > state.afn->AirflowNetworkLinkageData(i).NodeHeights[0])) break;
4064 2 : j += 9;
4065 2 : HSt(k) = 0.0;
4066 2 : if (HSt(k - 1) < 0.0) HSt(k) = HSt(k - 1);
4067 : }
4068 :
4069 : // Calculate Rho and stack pressure for every StartHeight of a layer which is
4070 : // within the top- and the bottom-height of the large opening.
4071 : while (true) {
4072 2 : ilayptr = 0;
4073 2 : if (Fromz == 0) ilayptr = 9;
4074 2 : if ((j > ilayptr) || (HSt(k) >= (state.afn->AirflowNetworkLinkageData(i).NodeHeights[0] + ActLOwnh)))
4075 2 : break; // Autodesk:BoundsViolation HSt(k) @ k>2
4076 0 : T = TzFrom;
4077 0 : X = XhzFrom;
4078 0 : lclimb(
4079 0 : state, G, RhoStd, HSt(k), T, X, DpF(k), Fromz, PzFrom, Pbz, RhoDrDummi); // Autodesk:BoundsViolation HSt(k) and DpF(k) @ k>2
4080 0 : RhoStF(k) = RhoStd; // Autodesk:BoundsViolation RhoStF(k) @ k>2
4081 0 : j += 9;
4082 0 : ++k; // Autodesk:Note k>2 now
4083 0 : HSt(k) = 0.0; // Autodesk:BoundsViolation @ k>2
4084 0 : if (HSt(k - 1) < 0.0) HSt(k) = HSt(k - 1); // Autodesk:BoundsViolation @ k>2
4085 : }
4086 : // Stack pressure difference and rho for top-height of the large opening
4087 2 : HSt(k) = state.afn->AirflowNetworkLinkageData(i).NodeHeights[0] + ActLOwnh; // Autodesk:BoundsViolation k>2 poss
4088 2 : T = TzFrom;
4089 2 : X = XhzFrom;
4090 2 : lclimb(state, G, RhoStd, HSt(k), T, X, DpF(k), Fromz, PzFrom, Pbz, RhoDrDummi); // Autodesk:BoundsViolation k>2 poss
4091 2 : RhoStF(k) = RhoStd; // Autodesk:BoundsViolation k >= 3 poss
4092 :
4093 4 : for (j = 1; j <= (k - 1); ++j) {
4094 2 : BetaStF(j) = (RhoStF(j + 1) - RhoStF(j)) / (HSt(j + 1) - HSt(j));
4095 : }
4096 : }
4097 :
4098 : // repeat procedure for the "To" node, DpT
4099 320459 : ilayptr = 0;
4100 320459 : if (Toz == 0) ilayptr = 1;
4101 320459 : j = ilayptr;
4102 : // Calculate Rho at link height only if we have large openings or layered zones.
4103 320459 : k = 1;
4104 320459 : lclimb(state, G, RhoLd(2), state.afn->AirflowNetworkLinkageData(i).NodeHeights[1], TempL2, Xhl2, DpT(k), Toz, PzTo, Pbz, RhoDrL(i, 2));
4105 320459 : Real64 RhoL2 = RhoLd(2);
4106 :
4107 : // For large openings calculate the stack pressure difference profile and the
4108 : // density profile within the the top- and the bottom- height of the large opening
4109 320459 : if (ActLOwnh > 0.0) {
4110 2 : HSt(k) = state.afn->AirflowNetworkLinkageData(i).NodeHeights[1];
4111 2 : RhoStT(k) = RhoL2;
4112 2 : ++k;
4113 2 : HSt(k) = 0.0;
4114 2 : if (HSt(k - 1) < 0.0) HSt(k) = HSt(k - 1);
4115 : while (true) {
4116 4 : ilayptr = 0;
4117 4 : if (Toz == 0) ilayptr = 9;
4118 4 : if ((j > ilayptr) || (HSt(k) > state.afn->AirflowNetworkLinkageData(i).NodeHeights[1])) break;
4119 2 : j += 9;
4120 2 : HSt(k) = 0.0;
4121 2 : if (HSt(k - 1) < 0.0) HSt(k) = HSt(k - 1);
4122 : }
4123 : // Calculate Rho and stack pressure for every StartHeight of a layer which is
4124 : // within the top- and the bottom-height of the large opening.
4125 : while (true) {
4126 2 : ilayptr = 0;
4127 2 : if (Toz == 0) ilayptr = 9;
4128 2 : if ((j > ilayptr) || (HSt(k) >= (state.afn->AirflowNetworkLinkageData(i).NodeHeights[1] + ActLOwnh)))
4129 2 : break; // Autodesk:BoundsViolation Hst(k) @ k>2
4130 0 : T = TzTo;
4131 0 : X = XhzTo;
4132 0 : lclimb(state, G, RhoStd, HSt(k), T, X, DpT(k), Toz, PzTo, Pbz, RhoDrDummi); // Autodesk:BoundsViolation HSt(k) and DpT(k) @ k>2
4133 0 : RhoStT(k) = RhoStd; // Autodesk:BoundsViolation RhoStT(k) @ k>2
4134 0 : j += 9;
4135 0 : ++k; // Autodesk:Note k>2 now
4136 0 : HSt(k) = 0.0; // Autodesk:BoundsViolation @ k>2
4137 0 : if (HSt(k - 1) < 0.0) HSt(k) = HSt(k - 1); // Autodesk:BoundsViolation @ k>2
4138 : }
4139 : // Stack pressure difference and rho for top-height of the large opening
4140 2 : HSt(k) = state.afn->AirflowNetworkLinkageData(i).NodeHeights[1] + ActLOwnh; // Autodesk:BoundsViolation k>2 poss
4141 2 : T = TzTo;
4142 2 : X = XhzTo;
4143 2 : lclimb(state, G, RhoStd, HSt(k), T, X, DpT(k), Toz, PzTo, Pbz, RhoDrDummi); // Autodesk:BoundsViolation k>2 poss
4144 2 : RhoStT(k) = RhoStd; // Autodesk:BoundsViolation k>2 poss
4145 :
4146 4 : for (j = 1; j <= (k - 1); ++j) {
4147 2 : BetaStT(j) = (RhoStT(j + 1) - RhoStT(j)) / (HSt(j + 1) - HSt(j));
4148 : }
4149 : }
4150 :
4151 : // CALCULATE STACK PRESSURE FOR THE PATH ITSELF for different flow directions
4152 320459 : H = double(state.afn->AirflowNetworkLinkageData(i).NodeHeights[1]) - double(state.afn->AirflowNetworkLinkageData(i).NodeHeights[0]);
4153 320459 : if (ll == 0 || ll == 3 || ll == 6) {
4154 320459 : H -= state.afn->AirflowNetworkNodeData(From).NodeHeight;
4155 : }
4156 320459 : if (ll < 3) {
4157 33810 : H += state.afn->AirflowNetworkNodeData(To).NodeHeight;
4158 : }
4159 :
4160 : // IF AIR FLOWS from "From" to "To"
4161 320459 : Pref = Pbz + PzFrom + DpF(1);
4162 320459 : DpP = psz(Pref, RhoLd(1), 0.0, 0.0, H, G);
4163 320459 : DpL(i, 1) = (DpF(1) - DpT(1) + DpP);
4164 :
4165 : // IF AIR FLOWS from "To" to "From"
4166 320459 : Pref = Pbz + PzTo + DpT(1);
4167 320459 : DpP = -psz(Pref, RhoLd(2), 0.0, 0.0, -H, G);
4168 320459 : DpL(i, 2) = (DpF(1) - DpT(1) + DpP);
4169 :
4170 320459 : if (Ltyp == iComponentTypeNum::DOP) {
4171 2 : Pprof = OpenNum * (NrInt + 2);
4172 2 : presprofile(state, i, Pprof, G, DpF, DpT, BetaStF, BetaStT, RhoStF, RhoStT, From, To, ActLh, Hfl(i));
4173 2 : ++OpenNum;
4174 : }
4175 : }
4176 16873 : }
4177 :
4178 1281840 : Real64 DetailedOpeningSolver::psz(Real64 const Pz0, // Pressure at altitude z0 [Pa]
4179 : Real64 const Rho0, // density at altitude z0 [kg/m3]
4180 : Real64 const beta, // density gradient [kg/m4]
4181 : Real64 const z0, // reference altitude [m]
4182 : Real64 const z, // altitude[m]
4183 : Real64 const g // gravity field strength [N/kg]
4184 : )
4185 : {
4186 :
4187 : // FUNCTION INFORMATION:
4188 : // AUTHOR Lixing Gu
4189 : // DATE WRITTEN Oct. 2005
4190 : // MODIFIED na
4191 : // RE-ENGINEERED This subroutine is revised based on psz function from COMIS
4192 :
4193 : // PURPOSE OF THIS SUBROUTINE:
4194 : // This function determines the pressure due to buoyancy in a stratified density environment
4195 :
4196 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
4197 : Real64 dz;
4198 : Real64 rho;
4199 :
4200 1281840 : dz = z - z0;
4201 1281840 : rho = (Rho0 + beta * dz / 2.0);
4202 1281840 : return -Pz0 * (1.0 - std::exp(-dz * rho * g / Pz0)); // Differential pressure from z to z0 [Pa]
4203 : }
4204 :
4205 640922 : void DetailedOpeningSolver::lclimb(EnergyPlusData &state,
4206 : Real64 const G, // gravity field strength [N/kg]
4207 : Real64 &Rho, // Density link level (initialized with rho zone) [kg/m3]
4208 : Real64 const Z, // Height of the link above the zone reference [m]
4209 : Real64 &T, // temperature at link level [C]
4210 : Real64 &X, // absolute humidity at link level [kg/kg]
4211 : Real64 &Dp, // Stackpressure to the linklevel [Pa]
4212 : int const zone, // Zone number
4213 : Real64 const PZ, // Zone Pressure (reflevel) [Pa]
4214 : Real64 const Pbz, // Barometric pressure at entrance level [Pa]
4215 : Real64 &RhoDr // Air density of dry air on the link level used
4216 : )
4217 : {
4218 :
4219 : // FUNCTION INFORMATION:
4220 : // AUTHOR Lixing Gu
4221 : // DATE WRITTEN Oct. 2005
4222 : // MODIFIED na
4223 : // RE-ENGINEERED This subroutine is revised based on subroutine IClimb from COMIS
4224 :
4225 : // PURPOSE OF THIS SUBROUTINE:
4226 : // This function the differential pressure from the reflevel in a zone To Z, the level of a link
4227 :
4228 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
4229 : Real64 H; // Start Height of the layer
4230 : Real64 BetaT; // Temperature gradient of this layer
4231 : Real64 BetaXfct; // Humidity gradient factor of this layer
4232 : Real64 BetaCfct; // Concentration 1 gradient factor of this layer
4233 : Real64 X0;
4234 : Real64 P;
4235 : Real64 Htop;
4236 : Real64 Hbot;
4237 : Real64 Rho0;
4238 : Real64 Rho1;
4239 : Real64 BetaRho;
4240 : static int L(0);
4241 : static int ilayptr(0);
4242 :
4243 640922 : Dp = 0.0;
4244 640922 : Rho0 = Rho;
4245 640922 : X0 = X;
4246 640922 : if (Z > 0.0) {
4247 : // initialize start values
4248 640922 : H = 0.0;
4249 640922 : BetaT = 0.0;
4250 640922 : BetaXfct = 0.0;
4251 640922 : BetaCfct = 0.0;
4252 640922 : BetaRho = 0.0;
4253 640922 : Hbot = 0.0;
4254 :
4255 640922 : while (H < 0.0) {
4256 : // loop until H>0 ; The start of the layer is above 0
4257 0 : BetaT = 0.0;
4258 0 : BetaXfct = 0.0;
4259 0 : BetaCfct = 0.0;
4260 0 : L += 9;
4261 0 : ilayptr = 0;
4262 0 : if (zone == 0) ilayptr = 9;
4263 0 : if (L >= ilayptr) {
4264 0 : H = Z + 1.0;
4265 : } else {
4266 0 : H = 0.0;
4267 : }
4268 : }
4269 :
4270 : // The link is in this layer also if it is on the top of it.
4271 :
4272 : while (true) {
4273 1281844 : if (H >= Z) {
4274 : // the link ends in this layer , we reached the final Dp and BetaRho
4275 640922 : Htop = Z;
4276 640922 : P = PZ + Dp;
4277 640922 : if (Htop != Hbot) {
4278 640922 : Rho0 = state.afn->properties.density(Pbz + P, T, X);
4279 640922 : T += (Htop - Hbot) * BetaT;
4280 640922 : X += (Htop - Hbot) * BetaXfct * X0;
4281 640922 : Rho1 = state.afn->properties.density(Pbz + P, T, X);
4282 640922 : BetaRho = (Rho1 - Rho0) / (Htop - Hbot);
4283 640922 : Dp += psz(Pbz + P, Rho0, BetaRho, Hbot, Htop, G);
4284 : }
4285 640922 : RhoDr = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4286 640922 : Rho = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4287 640922 : return;
4288 :
4289 : } else {
4290 : // bottom of the layer is below Z (Z above ref)
4291 640922 : Htop = H;
4292 : // P is the pressure up to the start height of the layer we just reached
4293 640922 : P = PZ + Dp;
4294 640922 : if (Htop != Hbot) {
4295 0 : Rho0 = state.afn->properties.density(Pbz + P, T, X);
4296 0 : T += (Htop - Hbot) * BetaT;
4297 0 : X += (Htop - Hbot) * BetaXfct * X0;
4298 0 : Rho1 = state.afn->properties.density(Pbz + P, T, X);
4299 0 : BetaRho = (Rho1 - Rho0) / (Htop - Hbot);
4300 0 : Dp += psz(Pbz + P, Rho0, BetaRho, Hbot, Htop, G);
4301 : }
4302 :
4303 640922 : RhoDr = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4304 640922 : Rho = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4305 :
4306 : // place current values Hbot and Beta's
4307 640922 : Hbot = H;
4308 640922 : BetaT = 0.0;
4309 640922 : BetaXfct = 0.0;
4310 640922 : BetaCfct = 0.0;
4311 640922 : L += 9;
4312 640922 : ilayptr = 0;
4313 640922 : if (zone == 0) ilayptr = 9;
4314 640922 : if (L >= ilayptr) {
4315 640922 : H = Z + 1.0;
4316 : } else {
4317 0 : H = 0.0;
4318 : }
4319 : }
4320 : }
4321 :
4322 : } else {
4323 : // This is the ELSE for negative linkheights Z below the refplane
4324 0 : H = 0.0;
4325 0 : BetaT = 0.0;
4326 0 : BetaXfct = 0.0;
4327 0 : BetaCfct = 0.0;
4328 0 : BetaRho = 0.0;
4329 0 : Htop = 0.0;
4330 0 : while (H > 0.0) {
4331 : // loop until H<0 ; The start of the layer is below the zone refplane
4332 0 : L -= 9;
4333 0 : ilayptr = 0;
4334 0 : if (zone == 0) ilayptr = 1;
4335 0 : if (L < ilayptr) {
4336 : // with H=Z (negative) this loop will exit, no data for interval Z-refplane
4337 0 : H = Z;
4338 0 : BetaT = 0.0;
4339 0 : BetaXfct = 0.0;
4340 0 : BetaCfct = 0.0;
4341 0 : BetaRho = 0.0;
4342 : } else {
4343 0 : H = 0.0;
4344 0 : BetaT = 0.0;
4345 0 : BetaXfct = 0.0;
4346 0 : BetaCfct = 0.0;
4347 : }
4348 : }
4349 :
4350 : // The link is in this layer also if it is on the bottom of it.
4351 : while (true) {
4352 0 : if (H <= Z) {
4353 0 : Hbot = Z;
4354 0 : P = PZ + Dp;
4355 0 : if (Htop != Hbot) {
4356 0 : Rho1 = state.afn->properties.density(Pbz + P, T, X);
4357 0 : T += (Hbot - Htop) * BetaT;
4358 0 : X += (Hbot - Htop) * BetaXfct * X0;
4359 0 : Rho0 = state.afn->properties.density(Pbz + P, T, X);
4360 0 : BetaRho = (Rho1 - Rho0) / (Htop - Hbot);
4361 0 : Dp -= psz(Pbz + P, Rho0, BetaRho, Hbot, Htop, G);
4362 : }
4363 0 : RhoDr = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4364 0 : Rho = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4365 0 : return;
4366 : } else {
4367 : // bottom of the layer is below Z (Z below ref)
4368 0 : Hbot = H;
4369 0 : P = PZ + Dp;
4370 0 : if (Htop != Hbot) {
4371 0 : Rho1 = state.afn->properties.density(Pbz + P, T, X);
4372 : // T,X,C calculated for the lower height
4373 0 : T += (Hbot - Htop) * BetaT;
4374 0 : X += (Hbot - Htop) * BetaXfct * X0;
4375 0 : Rho0 = state.afn->properties.density(Pbz + P, T, X);
4376 0 : BetaRho = (Rho1 - Rho0) / (Htop - Hbot);
4377 0 : Dp -= psz(Pbz + P, Rho0, BetaRho, Hbot, Htop, G);
4378 : }
4379 0 : RhoDr = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4380 0 : Rho = state.afn->properties.density(Pbz + PZ + Dp, T, X);
4381 :
4382 : // place current values Hbot and Beta's
4383 0 : Htop = H;
4384 0 : L -= 9;
4385 0 : ilayptr = 0;
4386 0 : if (zone == 0) ilayptr = 1;
4387 0 : if (L < ilayptr) {
4388 0 : H = Z - 1.0;
4389 0 : BetaT = 0.0;
4390 0 : BetaXfct = 0.0;
4391 0 : BetaCfct = 0.0;
4392 : } else {
4393 0 : H = 0.0;
4394 0 : BetaT = 0.0;
4395 0 : BetaXfct = 0.0;
4396 0 : BetaCfct = 0.0;
4397 : }
4398 : }
4399 : // ENDIF H<Z
4400 : }
4401 : }
4402 : }
4403 :
4404 : } // namespace AirflowNetwork
4405 :
4406 : } // namespace EnergyPlus
|
