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