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