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 : // C++ Headers
49 : #include <cmath>
50 :
51 : // ObjexxFCL Headers
52 : #include <ObjexxFCL/Array.functions.hh>
53 : #include <ObjexxFCL/Fmath.hh>
54 :
55 : // EnergyPlus Headers
56 : #include <EnergyPlus/Data/EnergyPlusData.hh>
57 : #include <EnergyPlus/DataEnvironment.hh>
58 : #include <EnergyPlus/DataHVACGlobals.hh>
59 : #include <EnergyPlus/DataHeatBalance.hh>
60 : #include <EnergyPlus/DataIPShortCuts.hh>
61 : #include <EnergyPlus/EarthTube.hh>
62 : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
63 : #include <EnergyPlus/OutputProcessor.hh>
64 : #include <EnergyPlus/Psychrometrics.hh>
65 : #include <EnergyPlus/ScheduleManager.hh>
66 : #include <EnergyPlus/UtilityRoutines.hh>
67 : #include <EnergyPlus/ZoneTempPredictorCorrector.hh>
68 :
69 : namespace EnergyPlus::EarthTube {
70 : // Module containing the data for Earth Tube system
71 :
72 : // MODULE INFORMATION:
73 : // AUTHOR Kwang Ho Lee
74 : // DATE WRITTEN November 2005
75 :
76 : // PURPOSE OF THIS MODULE:
77 : // To encapsulate the data and algorithyms required to manage the EarthTube System Component
78 :
79 : // REFERENCES:
80 : // 1. M. Krarti, "Analytical Model to Predict Annual Soil Surface Temperature Variation",
81 : // Journal of Solar Energy Engineering 117, 1995, pp 91-99
82 : // 2. K. Labs In: J. Cook, editor, "Passive Cooling",
83 : // Cambridge Massachusetts, MIT Press, 1989, pp 206-212
84 :
85 : // This is an interesting one. The actual members of the enum are never explicitly used
86 : // The enum is used in a getEnumValue call to determine what was found in GetInput
87 : // The value is then used as an array index to lookup thermal conductivity and such from some std::arrays
88 : // So the IDE thinks these are unused, and I'm not sure the best way to hint that they sorta aren't
89 : enum class SoilType
90 : {
91 : Invalid = -1,
92 : HeavyAndSat,
93 : HeavyAndDamp,
94 : HeavyAndDry,
95 : LightAndDry,
96 : Num
97 : };
98 :
99 : int totEarthTube = 0;
100 :
101 : constexpr std::array<std::string_view, static_cast<int>(Ventilation::Num)> ventilationNamesUC = {"NATURAL", "INTAKE", "EXHAUST"};
102 : constexpr std::array<std::string_view, static_cast<int>(SoilType::Num)> soilTypeNamesUC = {
103 : "HEAVYANDSATURATED", "HEAVYANDDAMP", "HEAVYANDDRY", "LIGHTANDDRY"};
104 : constexpr std::array<std::string_view, static_cast<int>(EarthTubeModelType::Num)> solutionTypeNamesUC = {"BASIC", "VERTICAL"};
105 :
106 3757673 : void ManageEarthTube(EnergyPlusData &state)
107 : {
108 :
109 : // SUBROUTINE INFORMATION:
110 : // AUTHOR Kwang Ho Lee
111 : // DATE WRITTEN November 2005
112 :
113 : // PURPOSE OF THIS SUBROUTINE:
114 : // This subroutine manages the simulation of EarthTube unit.
115 : // This driver manages the calls to all of
116 : // the other drivers and simulation algorithms.
117 :
118 : // Obtains and Allocates heat balance related parameters from input file
119 3757673 : if (state.dataEarthTube->GetInputFlag) {
120 768 : bool ErrorsFound = false;
121 768 : GetEarthTube(state, ErrorsFound);
122 768 : state.dataEarthTube->GetInputFlag = false;
123 : }
124 :
125 3757673 : if (state.dataEarthTube->EarthTubeSys.empty()) {
126 3753377 : return;
127 : }
128 :
129 4296 : initEarthTubeVertical(state);
130 :
131 4296 : CalcEarthTube(state);
132 :
133 4296 : ReportEarthTube(state);
134 : }
135 :
136 768 : void GetEarthTube(EnergyPlusData &state, bool &ErrorsFound) // If errors found in input
137 : {
138 :
139 : // SUBROUTINE INFORMATION:
140 : // AUTHOR Kwang Ho Lee
141 : // DATE WRITTEN November 2005
142 :
143 : // PURPOSE OF THIS SUBROUTINE:
144 : // This subroutine obtains input data for EarthTube units and
145 : // stores it in the EarthTube data structure.
146 :
147 : static constexpr std::string_view routineName = "GetEarthTube";
148 :
149 : // SUBROUTINE PARAMETER DEFINITIONS:
150 768 : Real64 constexpr EarthTubeTempLimit(100.0); // degrees Celsius
151 :
152 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
153 : int NumAlpha;
154 : int NumNumber;
155 : int IOStat;
156 : int Loop;
157 768 : Array1D_bool RepVarSet;
158 :
159 768 : RepVarSet.dimension(state.dataGlobal->NumOfZones, true);
160 :
161 : // Following used for reporting
162 768 : state.dataEarthTube->ZnRptET.allocate(state.dataGlobal->NumOfZones);
163 :
164 768 : auto &s_ipsc = state.dataIPShortCut;
165 :
166 768 : s_ipsc->cCurrentModuleObject = "ZoneEarthtube:Parameters";
167 :
168 768 : int totEarthTubePars = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, s_ipsc->cCurrentModuleObject);
169 :
170 768 : state.dataEarthTube->EarthTubePars.allocate(totEarthTubePars);
171 :
172 770 : for (Loop = 1; Loop <= totEarthTubePars; ++Loop) {
173 2 : auto &thisEarthTubePars = state.dataEarthTube->EarthTubePars(Loop);
174 4 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
175 2 : s_ipsc->cCurrentModuleObject,
176 : Loop,
177 2 : s_ipsc->cAlphaArgs,
178 : NumAlpha,
179 2 : s_ipsc->rNumericArgs,
180 : NumNumber,
181 : IOStat,
182 2 : s_ipsc->lNumericFieldBlanks,
183 2 : s_ipsc->lAlphaFieldBlanks,
184 2 : s_ipsc->cAlphaFieldNames,
185 2 : s_ipsc->cNumericFieldNames);
186 :
187 2 : thisEarthTubePars.nameParameters = s_ipsc->cAlphaArgs(1);
188 : // Check to make sure name is unique
189 3 : for (int otherParams = 1; otherParams < Loop; ++otherParams) {
190 1 : if (Util::SameString(thisEarthTubePars.nameParameters, state.dataEarthTube->EarthTubePars(otherParams).nameParameters)) {
191 0 : ShowSevereError(
192 : state,
193 0 : format("{}: {} = {} is not a unique name.", s_ipsc->cCurrentModuleObject, s_ipsc->cAlphaFieldNames(1), s_ipsc->cAlphaArgs(1)));
194 0 : ShowContinueError(state, format("Check the other {} names for a duplicate.", s_ipsc->cCurrentModuleObject));
195 0 : ErrorsFound = true;
196 : }
197 : }
198 :
199 2 : thisEarthTubePars.numNodesAbove = s_ipsc->rNumericArgs(1);
200 2 : thisEarthTubePars.numNodesBelow = s_ipsc->rNumericArgs(2);
201 2 : thisEarthTubePars.dimBoundAbove = s_ipsc->rNumericArgs(3);
202 2 : thisEarthTubePars.dimBoundBelow = s_ipsc->rNumericArgs(4);
203 2 : thisEarthTubePars.width = s_ipsc->rNumericArgs(5);
204 : }
205 :
206 768 : s_ipsc->cCurrentModuleObject = "ZoneEarthtube";
207 768 : totEarthTube = state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, s_ipsc->cCurrentModuleObject);
208 :
209 768 : state.dataEarthTube->EarthTubeSys.allocate(totEarthTube);
210 :
211 774 : for (Loop = 1; Loop <= totEarthTube; ++Loop) {
212 6 : auto &thisEarthTube = state.dataEarthTube->EarthTubeSys(Loop);
213 12 : state.dataInputProcessing->inputProcessor->getObjectItem(state,
214 6 : s_ipsc->cCurrentModuleObject,
215 : Loop,
216 6 : s_ipsc->cAlphaArgs,
217 : NumAlpha,
218 6 : s_ipsc->rNumericArgs,
219 : NumNumber,
220 : IOStat,
221 6 : s_ipsc->lNumericFieldBlanks,
222 6 : s_ipsc->lAlphaFieldBlanks,
223 6 : s_ipsc->cAlphaFieldNames,
224 6 : s_ipsc->cNumericFieldNames);
225 :
226 6 : ErrorObjectHeader eoh{routineName, s_ipsc->cCurrentModuleObject, s_ipsc->cAlphaArgs(1)};
227 :
228 : // First Alpha is Zone Name
229 6 : thisEarthTube.ZonePtr = Util::FindItemInList(s_ipsc->cAlphaArgs(1), state.dataHeatBal->Zone);
230 6 : if (thisEarthTube.ZonePtr == 0) {
231 0 : ShowSevereError(state, format("{}: {} not found={}", s_ipsc->cCurrentModuleObject, s_ipsc->cAlphaFieldNames(1), s_ipsc->cAlphaArgs(1)));
232 0 : ErrorsFound = true;
233 : }
234 :
235 : // Second Alpha is Schedule Name
236 6 : if (s_ipsc->lAlphaFieldBlanks(2)) {
237 0 : ShowSevereEmptyField(state, eoh, s_ipsc->cAlphaFieldNames(2));
238 0 : ErrorsFound = true;
239 6 : } else if ((thisEarthTube.availSched = Sched::GetSchedule(state, s_ipsc->cAlphaArgs(2))) == nullptr) {
240 0 : ShowSevereItemNotFound(state, eoh, s_ipsc->cAlphaFieldNames(2), s_ipsc->cAlphaArgs(2));
241 0 : ErrorsFound = true;
242 : }
243 :
244 : // Overall parameters and their limits
245 6 : thisEarthTube.DesignLevel = s_ipsc->rNumericArgs(1);
246 :
247 6 : thisEarthTube.MinTemperature = s_ipsc->rNumericArgs(2);
248 6 : if ((thisEarthTube.MinTemperature < -EarthTubeTempLimit) || (thisEarthTube.MinTemperature > EarthTubeTempLimit)) {
249 0 : ShowSevereError(state,
250 0 : format("{}: {}={} must have a minimum temperature between -{:.0R}C and {:.0R}C",
251 0 : s_ipsc->cCurrentModuleObject,
252 0 : s_ipsc->cAlphaFieldNames(1),
253 0 : s_ipsc->cAlphaArgs(1),
254 : EarthTubeTempLimit,
255 : EarthTubeTempLimit));
256 0 : ShowContinueError(state, format("Entered value={:.0R}", thisEarthTube.MinTemperature));
257 0 : ErrorsFound = true;
258 : }
259 :
260 6 : thisEarthTube.MaxTemperature = s_ipsc->rNumericArgs(3);
261 6 : if ((thisEarthTube.MaxTemperature < -EarthTubeTempLimit) || (thisEarthTube.MaxTemperature > EarthTubeTempLimit)) {
262 0 : ShowSevereError(state,
263 0 : format("{}: {}={} must have a maximum temperature between -{:.0R}C and {:.0R}C",
264 0 : s_ipsc->cCurrentModuleObject,
265 0 : s_ipsc->cAlphaFieldNames(1),
266 0 : s_ipsc->cAlphaArgs(1),
267 : EarthTubeTempLimit,
268 : EarthTubeTempLimit));
269 0 : ShowContinueError(state, format("Entered value={:.0R}", thisEarthTube.MaxTemperature));
270 0 : ErrorsFound = true;
271 : }
272 :
273 6 : thisEarthTube.DelTemperature = s_ipsc->rNumericArgs(4); // 3/12/03 Negative del temp now allowed COP
274 :
275 : // if we have a blank, then just set it to the Natural type, otherwise, search on it
276 6 : if (s_ipsc->cAlphaArgs(3).empty()) {
277 0 : thisEarthTube.FanType = Ventilation::Natural;
278 : } else {
279 6 : thisEarthTube.FanType = static_cast<Ventilation>(getEnumValue(ventilationNamesUC, s_ipsc->cAlphaArgs(3)));
280 6 : if (thisEarthTube.FanType == Ventilation::Invalid) {
281 0 : ShowSevereInvalidKey(state, eoh, s_ipsc->cAlphaFieldNames(3), s_ipsc->cAlphaArgs(3));
282 0 : ErrorsFound = true;
283 : }
284 : }
285 :
286 6 : thisEarthTube.FanPressure = s_ipsc->rNumericArgs(5);
287 6 : if (thisEarthTube.FanPressure < 0.0) {
288 0 : ShowSevereError(state,
289 0 : format("{}: {}={}, {} must be positive, entered value={:.2R}",
290 0 : s_ipsc->cCurrentModuleObject,
291 0 : s_ipsc->cAlphaFieldNames(1),
292 0 : s_ipsc->cAlphaArgs(1),
293 0 : s_ipsc->cNumericFieldNames(5),
294 0 : thisEarthTube.FanPressure));
295 0 : ErrorsFound = true;
296 : }
297 :
298 6 : thisEarthTube.FanEfficiency = s_ipsc->rNumericArgs(6);
299 6 : if ((thisEarthTube.FanEfficiency <= 0.0) || (thisEarthTube.FanEfficiency > 1.0)) {
300 0 : ShowSevereError(state,
301 0 : format("{}: {}={}, {} must be greater than zero and less than or equal to one, entered value={:.2R}",
302 0 : s_ipsc->cCurrentModuleObject,
303 0 : s_ipsc->cAlphaFieldNames(1),
304 0 : s_ipsc->cAlphaArgs(1),
305 0 : s_ipsc->cNumericFieldNames(6),
306 0 : thisEarthTube.FanEfficiency));
307 0 : ErrorsFound = true;
308 : }
309 :
310 6 : thisEarthTube.r1 = s_ipsc->rNumericArgs(7);
311 6 : if (thisEarthTube.r1 <= 0.0) {
312 0 : ShowSevereError(state,
313 0 : format("{}: {}={}, {} must be positive, entered value={:.2R}",
314 0 : s_ipsc->cCurrentModuleObject,
315 0 : s_ipsc->cAlphaFieldNames(1),
316 0 : s_ipsc->cAlphaArgs(1),
317 0 : s_ipsc->cNumericFieldNames(7),
318 0 : thisEarthTube.r1));
319 0 : ErrorsFound = true;
320 : }
321 :
322 6 : thisEarthTube.r2 = s_ipsc->rNumericArgs(8);
323 6 : if (thisEarthTube.r2 <= 0.0) {
324 0 : ShowSevereError(state,
325 0 : format("{}: {}={}, {} must be positive, entered value={:.2R}",
326 0 : s_ipsc->cCurrentModuleObject,
327 0 : s_ipsc->cAlphaFieldNames(1),
328 0 : s_ipsc->cAlphaArgs(1),
329 0 : s_ipsc->cNumericFieldNames(8),
330 0 : thisEarthTube.r2));
331 0 : ErrorsFound = true;
332 : }
333 :
334 6 : thisEarthTube.r3 = 2.0 * thisEarthTube.r1;
335 :
336 6 : thisEarthTube.PipeLength = s_ipsc->rNumericArgs(9);
337 6 : if (thisEarthTube.PipeLength <= 0.0) {
338 0 : ShowSevereError(state,
339 0 : format("{}: {}={}, {} must be positive, entered value={:.2R}",
340 0 : s_ipsc->cCurrentModuleObject,
341 0 : s_ipsc->cAlphaFieldNames(1),
342 0 : s_ipsc->cAlphaArgs(1),
343 0 : s_ipsc->cNumericFieldNames(9),
344 0 : thisEarthTube.PipeLength));
345 0 : ErrorsFound = true;
346 : }
347 :
348 6 : thisEarthTube.PipeThermCond = s_ipsc->rNumericArgs(10);
349 6 : if (thisEarthTube.PipeThermCond <= 0.0) {
350 0 : ShowSevereError(state,
351 0 : format("{}: {}={}, {} must be positive, entered value={:.2R}",
352 0 : s_ipsc->cCurrentModuleObject,
353 0 : s_ipsc->cAlphaFieldNames(1),
354 0 : s_ipsc->cAlphaArgs(1),
355 0 : s_ipsc->cNumericFieldNames(10),
356 0 : thisEarthTube.PipeThermCond));
357 0 : ErrorsFound = true;
358 : }
359 :
360 6 : thisEarthTube.z = s_ipsc->rNumericArgs(11);
361 6 : if (thisEarthTube.z <= 0.0) {
362 0 : ShowSevereError(state,
363 0 : format("{}: {}={}, {} must be positive, entered value={:.2R}",
364 0 : s_ipsc->cCurrentModuleObject,
365 0 : s_ipsc->cAlphaFieldNames(1),
366 0 : s_ipsc->cAlphaArgs(1),
367 0 : s_ipsc->cNumericFieldNames(11),
368 0 : thisEarthTube.z));
369 0 : ErrorsFound = true;
370 : }
371 6 : if (thisEarthTube.z <= (thisEarthTube.r1 + thisEarthTube.r2 + thisEarthTube.r3)) {
372 : // Note that code in initEarthTubeVertical assumes that this check remains in place--if this ever gets changed,
373 : // code in initEarthTubeVertical must be modified
374 0 : ShowSevereError(state,
375 0 : format("{}: {}={}, {} must be greater than 3*{} + {} entered value={:.2R} ref sum={:.2R}",
376 0 : s_ipsc->cCurrentModuleObject,
377 0 : s_ipsc->cAlphaFieldNames(1),
378 0 : s_ipsc->cAlphaArgs(1),
379 0 : s_ipsc->cNumericFieldNames(11),
380 0 : s_ipsc->cNumericFieldNames(7),
381 0 : s_ipsc->cNumericFieldNames(8),
382 0 : thisEarthTube.z,
383 0 : thisEarthTube.r1 + thisEarthTube.r2 + thisEarthTube.r3));
384 0 : ErrorsFound = true;
385 : }
386 :
387 6 : SoilType soilType = static_cast<SoilType>(getEnumValue(soilTypeNamesUC, s_ipsc->cAlphaArgs(4)));
388 6 : constexpr std::array<Real64, static_cast<int>(SoilType::Num)> thermalDiffusivity = {0.0781056, 0.055728, 0.0445824, 0.024192};
389 6 : constexpr std::array<Real64, static_cast<int>(SoilType::Num)> thermalConductivity = {2.42, 1.3, 0.865, 0.346};
390 6 : if (soilType == SoilType::Invalid) {
391 0 : ShowSevereInvalidKey(state, eoh, s_ipsc->cAlphaFieldNames(4), s_ipsc->cAlphaArgs(4));
392 0 : ErrorsFound = true;
393 : } else {
394 6 : thisEarthTube.SoilThermDiff = thermalDiffusivity[static_cast<int>(soilType)];
395 6 : thisEarthTube.SoilThermCond = thermalConductivity[static_cast<int>(soilType)];
396 : }
397 :
398 6 : thisEarthTube.AverSoilSurTemp = s_ipsc->rNumericArgs(12);
399 6 : thisEarthTube.ApmlSoilSurTemp = s_ipsc->rNumericArgs(13);
400 6 : thisEarthTube.SoilSurPhaseConst = int(s_ipsc->rNumericArgs(14));
401 :
402 : // Override any user input for cases where natural ventilation is being used
403 6 : if (thisEarthTube.FanType == Ventilation::Natural) {
404 2 : thisEarthTube.FanPressure = 0.0;
405 2 : thisEarthTube.FanEfficiency = 1.0;
406 : }
407 :
408 6 : thisEarthTube.ConstantTermCoef = s_ipsc->rNumericArgs(15);
409 6 : thisEarthTube.TemperatureTermCoef = s_ipsc->rNumericArgs(16);
410 6 : thisEarthTube.VelocityTermCoef = s_ipsc->rNumericArgs(17);
411 6 : thisEarthTube.VelocitySQTermCoef = s_ipsc->rNumericArgs(18);
412 :
413 : // cAlphaArgs(5)--Model type: basic or vertical
414 : // only process cAlphaArgs(6) if cAlphaArgs(5) is "Vertical"
415 6 : if (s_ipsc->cAlphaArgs(5).empty()) {
416 0 : thisEarthTube.ModelType = EarthTubeModelType::Basic;
417 : } else {
418 6 : thisEarthTube.ModelType = static_cast<EarthTubeModelType>(getEnumValue(solutionTypeNamesUC, s_ipsc->cAlphaArgs(5)));
419 6 : if (thisEarthTube.ModelType == EarthTubeModelType::Invalid) {
420 0 : ShowSevereInvalidKey(state, eoh, s_ipsc->cAlphaFieldNames(5), s_ipsc->cAlphaArgs(5));
421 0 : ErrorsFound = true;
422 : }
423 : }
424 :
425 6 : if (thisEarthTube.ModelType == EarthTubeModelType::Vertical) {
426 2 : thisEarthTube.r3 = 0.0; // Vertical model does not use this parameter--reset to zero (keep because r3=0 necessary so Rs=0 in calc routine)
427 : // Process the parameters based on the name (link via index)
428 2 : thisEarthTube.vertParametersPtr = 0;
429 3 : for (int parIndex = 1; parIndex <= totEarthTubePars; ++parIndex) {
430 3 : if (Util::SameString(s_ipsc->cAlphaArgs(6), state.dataEarthTube->EarthTubePars(parIndex).nameParameters)) {
431 2 : thisEarthTube.vertParametersPtr = parIndex;
432 2 : break;
433 : }
434 : }
435 2 : if (thisEarthTube.vertParametersPtr == 0) { // didn't find a match
436 0 : ShowSevereItemNotFound(state, eoh, s_ipsc->cAlphaFieldNames(6), s_ipsc->cAlphaArgs(6));
437 0 : ErrorsFound = true;
438 : }
439 : }
440 :
441 6 : if (thisEarthTube.ZonePtr > 0) {
442 6 : if (RepVarSet(thisEarthTube.ZonePtr)) {
443 6 : RepVarSet(thisEarthTube.ZonePtr) = false;
444 6 : auto &zone = state.dataHeatBal->Zone(thisEarthTube.ZonePtr);
445 6 : auto &thisZnRptET = state.dataEarthTube->ZnRptET(thisEarthTube.ZonePtr);
446 :
447 12 : SetupOutputVariable(state,
448 : "Earth Tube Zone Sensible Cooling Energy",
449 : Constant::Units::J,
450 6 : thisZnRptET.EarthTubeHeatLoss,
451 : OutputProcessor::TimeStepType::System,
452 : OutputProcessor::StoreType::Sum,
453 6 : zone.Name);
454 12 : SetupOutputVariable(state,
455 : "Earth Tube Zone Sensible Cooling Rate",
456 : Constant::Units::W,
457 6 : thisZnRptET.EarthTubeHeatLossRate,
458 : OutputProcessor::TimeStepType::System,
459 : OutputProcessor::StoreType::Average,
460 6 : zone.Name);
461 12 : SetupOutputVariable(state,
462 : "Earth Tube Zone Sensible Heating Energy",
463 : Constant::Units::J,
464 6 : thisZnRptET.EarthTubeHeatGain,
465 : OutputProcessor::TimeStepType::System,
466 : OutputProcessor::StoreType::Sum,
467 6 : zone.Name);
468 12 : SetupOutputVariable(state,
469 : "Earth Tube Zone Sensible Heating Rate",
470 : Constant::Units::W,
471 6 : thisZnRptET.EarthTubeHeatGainRate,
472 : OutputProcessor::TimeStepType::System,
473 : OutputProcessor::StoreType::Average,
474 6 : zone.Name);
475 12 : SetupOutputVariable(state,
476 : "Earth Tube Air Flow Volume",
477 : Constant::Units::m3,
478 6 : thisZnRptET.EarthTubeVolume,
479 : OutputProcessor::TimeStepType::System,
480 : OutputProcessor::StoreType::Sum,
481 6 : zone.Name);
482 12 : SetupOutputVariable(state,
483 : "Earth Tube Current Density Air Volume Flow Rate",
484 : Constant::Units::m3_s,
485 6 : thisZnRptET.EarthTubeVolFlowRate,
486 : OutputProcessor::TimeStepType::System,
487 : OutputProcessor::StoreType::Average,
488 6 : zone.Name);
489 12 : SetupOutputVariable(state,
490 : "Earth Tube Standard Density Air Volume Flow Rate",
491 : Constant::Units::m3_s,
492 6 : thisZnRptET.EarthTubeVolFlowRateStd,
493 : OutputProcessor::TimeStepType::System,
494 : OutputProcessor::StoreType::Average,
495 6 : zone.Name);
496 12 : SetupOutputVariable(state,
497 : "Earth Tube Air Flow Mass",
498 : Constant::Units::kg,
499 6 : thisZnRptET.EarthTubeMass,
500 : OutputProcessor::TimeStepType::System,
501 : OutputProcessor::StoreType::Sum,
502 6 : zone.Name);
503 12 : SetupOutputVariable(state,
504 : "Earth Tube Air Mass Flow Rate",
505 : Constant::Units::kg_s,
506 6 : thisZnRptET.EarthTubeMassFlowRate,
507 : OutputProcessor::TimeStepType::System,
508 : OutputProcessor::StoreType::Average,
509 6 : zone.Name);
510 12 : SetupOutputVariable(state,
511 : "Earth Tube Water Mass Flow Rate",
512 : Constant::Units::kg_s,
513 6 : thisZnRptET.EarthTubeWaterMassFlowRate,
514 : OutputProcessor::TimeStepType::System,
515 : OutputProcessor::StoreType::Average,
516 6 : zone.Name);
517 12 : SetupOutputVariable(state,
518 : "Earth Tube Fan Electricity Energy",
519 : Constant::Units::J,
520 6 : thisZnRptET.EarthTubeFanElec,
521 : OutputProcessor::TimeStepType::System,
522 : OutputProcessor::StoreType::Sum,
523 6 : zone.Name,
524 : Constant::eResource::Electricity,
525 : OutputProcessor::Group::Building);
526 12 : SetupOutputVariable(state,
527 : "Earth Tube Fan Electricity Rate",
528 : Constant::Units::W,
529 6 : thisZnRptET.EarthTubeFanElecPower,
530 : OutputProcessor::TimeStepType::System,
531 : OutputProcessor::StoreType::Average,
532 6 : zone.Name);
533 12 : SetupOutputVariable(state,
534 : "Earth Tube Zone Inlet Air Temperature",
535 : Constant::Units::C,
536 6 : thisZnRptET.EarthTubeAirTemp,
537 : OutputProcessor::TimeStepType::System,
538 : OutputProcessor::StoreType::Average,
539 6 : zone.Name);
540 12 : SetupOutputVariable(state,
541 : "Earth Tube Ground Interface Temperature",
542 : Constant::Units::C,
543 6 : thisEarthTube.GroundTempt,
544 : OutputProcessor::TimeStepType::System,
545 : OutputProcessor::StoreType::Average,
546 6 : zone.Name);
547 12 : SetupOutputVariable(state,
548 : "Earth Tube Outdoor Air Heat Transfer Rate",
549 : Constant::Units::W,
550 6 : thisZnRptET.EarthTubeOATreatmentPower,
551 : OutputProcessor::TimeStepType::System,
552 : OutputProcessor::StoreType::Average,
553 6 : zone.Name);
554 12 : SetupOutputVariable(state,
555 : "Earth Tube Zone Inlet Wet Bulb Temperature",
556 : Constant::Units::C,
557 6 : thisZnRptET.EarthTubeWetBulbTemp,
558 : OutputProcessor::TimeStepType::System,
559 : OutputProcessor::StoreType::Average,
560 6 : zone.Name);
561 12 : SetupOutputVariable(state,
562 : "Earth Tube Zone Inlet Humidity Ratio",
563 : Constant::Units::kgWater_kgDryAir,
564 6 : thisZnRptET.EarthTubeHumRat,
565 : OutputProcessor::TimeStepType::System,
566 : OutputProcessor::StoreType::Average,
567 6 : zone.Name);
568 : }
569 : }
570 : }
571 :
572 768 : CheckEarthTubesInZones(state, s_ipsc->cAlphaArgs(1), s_ipsc->cCurrentModuleObject, ErrorsFound);
573 :
574 768 : if (ErrorsFound) {
575 0 : ShowFatalError(state, format("{}: Errors getting input. Program terminates.", s_ipsc->cCurrentModuleObject));
576 : }
577 768 : }
578 :
579 768 : void CheckEarthTubesInZones(EnergyPlusData &state,
580 : std::string const &ZoneName, // name of zone for error reporting
581 : std::string_view FieldName, // name of earth tube in input
582 : bool &ErrorsFound // Found a problem
583 : )
584 : {
585 : // Check to make sure there is only one earth tube statement per zone
586 768 : int numEarthTubes = (int)state.dataEarthTube->EarthTubeSys.size();
587 772 : for (int Loop = 1; Loop <= numEarthTubes - 1; ++Loop) {
588 10 : for (int Loop1 = Loop + 1; Loop1 <= numEarthTubes; ++Loop1) {
589 6 : if (state.dataEarthTube->EarthTubeSys(Loop).ZonePtr == state.dataEarthTube->EarthTubeSys(Loop1).ZonePtr) {
590 0 : ShowSevereError(state, format("{} has more than one {} associated with it.", ZoneName, FieldName));
591 0 : ShowContinueError(state, format("Only one {} is allowed per zone. Check the definitions of {}", FieldName, FieldName));
592 0 : ShowContinueError(state, "in your input file and make sure that there is only one defined for each zone.");
593 0 : ErrorsFound = true;
594 : }
595 : }
596 : }
597 768 : }
598 :
599 4296 : void initEarthTubeVertical(EnergyPlusData &state)
600 : {
601 4296 : if (state.dataEarthTube->initFirstTime) {
602 2 : state.dataEarthTube->initFirstTime = false;
603 8 : for (int etNum = 1; etNum <= totEarthTube; ++etNum) {
604 6 : auto &thisEarthTube = state.dataEarthTube->EarthTubeSys(etNum);
605 6 : if (thisEarthTube.ModelType != EarthTubeModelType::Vertical) {
606 4 : continue; // Skip earth tubes that do not use vertical solution
607 : }
608 2 : auto &thisEarthTubeParams = state.dataEarthTube->EarthTubePars(thisEarthTube.vertParametersPtr);
609 2 : thisEarthTube.totNodes = thisEarthTubeParams.numNodesAbove + thisEarthTubeParams.numNodesBelow + 1;
610 2 : thisEarthTube.aCoeff.resize(thisEarthTube.totNodes);
611 2 : thisEarthTube.bCoeff.resize(thisEarthTube.totNodes);
612 2 : thisEarthTube.cCoeff.resize(thisEarthTube.totNodes);
613 2 : thisEarthTube.cCoeff0.resize(thisEarthTube.totNodes);
614 2 : thisEarthTube.dCoeff.resize(thisEarthTube.totNodes);
615 2 : thisEarthTube.cPrime.resize(thisEarthTube.totNodes);
616 2 : thisEarthTube.dPrime.resize(thisEarthTube.totNodes);
617 2 : thisEarthTube.cPrime0.resize(thisEarthTube.totNodes);
618 2 : thisEarthTube.tCurrent.resize(thisEarthTube.totNodes);
619 2 : thisEarthTube.tLast.resize(thisEarthTube.totNodes);
620 2 : thisEarthTube.depthNode.resize(thisEarthTube.totNodes);
621 2 : thisEarthTube.tUndist.resize(thisEarthTube.totNodes);
622 2 : Real64 thickBase = (thisEarthTube.z - 3.0 * thisEarthTube.r1);
623 2 : Real64 thickTop = thickBase * thisEarthTubeParams.dimBoundAbove / float(thisEarthTubeParams.numNodesAbove);
624 2 : Real64 thickBottom = thickBase * thisEarthTubeParams.dimBoundBelow / float(thisEarthTubeParams.numNodesBelow);
625 2 : Real64 thickEarthTube = 4.0 * thisEarthTube.r1;
626 2 : Real64 deltat = state.dataGlobal->TimeStepZone;
627 2 : Real64 thermDiff = thisEarthTube.SoilThermDiff / Constant::rHoursInDay; // convert to "per hour" from "per day"
628 :
629 : // Node equations determine the _Coeff terms--see Engineering Referenve for details on these equation types
630 : // Note that node numbers are shifted for c++ arrays that go from 0 to numNodes-1.
631 : // Node Type 1 (Top Node)
632 2 : Real64 commonTerm = thermDiff * deltat / (thickTop * thickTop);
633 2 : thisEarthTube.aCoeff[0] = 0.0; // no a0 value
634 2 : thisEarthTube.bCoeff[0] = 1.0 + 3.0 * commonTerm;
635 2 : thisEarthTube.cCoeff[0] = -1.0 * commonTerm;
636 2 : thisEarthTube.dMult0 = 2.0 * commonTerm; // does not include temperatures (upper boundary or previous time step)--added later
637 : // Node Type 2 (Generic Top Section Node)
638 13 : for (int nodeNum = 1; nodeNum <= thisEarthTubeParams.numNodesAbove - 2; ++nodeNum) {
639 11 : thisEarthTube.aCoeff[nodeNum] = -1.0 * commonTerm;
640 11 : thisEarthTube.bCoeff[nodeNum] = 1.0 + 2.0 * commonTerm;
641 11 : thisEarthTube.cCoeff[nodeNum] = -1.0 * commonTerm;
642 : }
643 : // Node Type 3 (Last Top Section Node)
644 2 : int thisNode = thisEarthTubeParams.numNodesAbove - 1;
645 2 : Real64 commonTerm2 = 2.0 * thermDiff * deltat / (thickTop + thickEarthTube) / thickTop;
646 2 : thisEarthTube.aCoeff[thisNode] = -1.0 * commonTerm;
647 2 : thisEarthTube.bCoeff[thisNode] = 1.0 + commonTerm + commonTerm2;
648 2 : thisEarthTube.cCoeff[thisNode] = -1.0 * commonTerm2;
649 : // Node Type 4 (Earth Tube Node)
650 2 : thisNode = thisEarthTubeParams.numNodesAbove;
651 2 : commonTerm = 2.0 * thermDiff * deltat / (thickTop + thickEarthTube) / thickEarthTube;
652 2 : commonTerm2 = 2.0 * thermDiff * deltat / (thickBottom + thickEarthTube) / thickEarthTube;
653 2 : thisEarthTube.aCoeff[thisNode] = -1.0 * commonTerm;
654 2 : thisEarthTube.bCoeff[thisNode] = 1.0 + commonTerm + commonTerm2; // does not include earth tube air flow term--added later
655 2 : thisEarthTube.cCoeff[thisNode] = -1.0 * commonTerm2;
656 : // Node Type 5 (First Bottom Section Node)
657 2 : thisNode = thisEarthTubeParams.numNodesAbove + 1;
658 2 : commonTerm = thermDiff * deltat / (thickBottom * thickBottom);
659 2 : commonTerm2 = 2.0 * thermDiff * deltat / (thickBottom + thickEarthTube) / thickBottom;
660 2 : thisEarthTube.aCoeff[thisNode] = -1.0 * commonTerm2;
661 2 : thisEarthTube.bCoeff[thisNode] = 1.0 + commonTerm + commonTerm2;
662 2 : thisEarthTube.cCoeff[thisNode] = -1.0 * commonTerm;
663 : // Node Type 6 (Generic Bottom Section Node)
664 11 : for (int nodeNum = thisNode + 1; nodeNum <= thisEarthTube.totNodes - 2; ++nodeNum) {
665 9 : thisEarthTube.aCoeff[nodeNum] = -1.0 * commonTerm;
666 9 : thisEarthTube.bCoeff[nodeNum] = 1.0 + 2.0 * commonTerm;
667 9 : thisEarthTube.cCoeff[nodeNum] = -1.0 * commonTerm;
668 : }
669 : // Node Type 7 (Last Bottom Section Node, i.e. Last Node)
670 2 : thisNode = thisEarthTube.totNodes - 1; // shifted due to c++ arrays that go from 0 to numNodes-1
671 2 : thisEarthTube.aCoeff[thisNode] = -1.0 * commonTerm;
672 2 : thisEarthTube.bCoeff[thisNode] = 1.0 + 3.0 * commonTerm;
673 2 : thisEarthTube.cCoeff[thisNode] = 0.0; // no cN value
674 2 : thisEarthTube.dMultN = 2.0 * commonTerm; // does not include previous temperature and earth tube air flow terms--added later
675 :
676 : // Initialize node temperatures using undisturbed temperature equation and node depths
677 : // First, nodes above the earth tube
678 2 : thisEarthTube.depthNode[thisEarthTubeParams.numNodesAbove - 1] = thisEarthTube.z - 0.5 * (thickEarthTube + thickTop);
679 15 : for (int nodeNum = thisEarthTubeParams.numNodesAbove - 2; nodeNum >= 0; --nodeNum) {
680 13 : thisEarthTube.depthNode[nodeNum] = thisEarthTube.depthNode[nodeNum + 1] - thickTop;
681 : }
682 : // Now, the earth tube node
683 2 : thisEarthTube.depthNode[thisEarthTubeParams.numNodesAbove] = thisEarthTube.z;
684 : // Finally the nodes below the earth tube
685 2 : thisEarthTube.depthNode[thisEarthTubeParams.numNodesAbove + 1] = thisEarthTube.z + 0.5 * (thickEarthTube + thickBottom);
686 13 : for (int nodeNumBelow = 2; nodeNumBelow <= thisEarthTubeParams.numNodesBelow; ++nodeNumBelow) {
687 11 : int nodeNum = thisEarthTubeParams.numNodesAbove + nodeNumBelow;
688 11 : thisEarthTube.depthNode[nodeNum] = thisEarthTube.depthNode[nodeNum - 1] + thickBottom;
689 : }
690 2 : thisEarthTube.depthUpperBound = thisEarthTube.depthNode[0] - 0.5 * thickTop;
691 2 : thisEarthTube.depthLowerBound = thisEarthTube.depthNode[thisEarthTube.totNodes - 1] + 0.5 * thickBottom;
692 :
693 : // Calculate constant part of air flow term at earth tube node. Note that diffusiity/conductivity = 1/(density*specific_heat)
694 2 : thisEarthTube.airFlowCoeff = state.dataGlobal->TimeStepZone * thermDiff / thisEarthTube.SoilThermCond / thickEarthTube /
695 2 : thisEarthTubeParams.width / thisEarthTube.PipeLength;
696 :
697 : // Calculate some initial values in the Thomas algorithm. This includes c' when effectiveness is zero (entire c').
698 : // For any other effectiveness, c' will be the same as c' when effectiveness for is zero for the nodes above the earth
699 : // tube. So, the c' for effectiveness of zero (cPrime0) can be reused as needed.
700 32 : for (int nodeNum = 0; nodeNum <= thisEarthTube.totNodes - 1; ++nodeNum) {
701 30 : thisEarthTube.cCoeff0[nodeNum] = thisEarthTube.cCoeff[nodeNum];
702 : }
703 2 : thisEarthTube.initCPrime0();
704 :
705 2 : auto &zone = state.dataHeatBal->Zone(thisEarthTube.ZonePtr);
706 32 : for (int nodeNum = 1; nodeNum <= thisEarthTube.totNodes; ++nodeNum) {
707 90 : SetupOutputVariable(state,
708 60 : format("Earth Tube Node Temperature {}", nodeNum),
709 : Constant::Units::C,
710 30 : thisEarthTube.tCurrent[nodeNum - 1],
711 : OutputProcessor::TimeStepType::Zone,
712 : OutputProcessor::StoreType::Average,
713 30 : zone.Name);
714 90 : SetupOutputVariable(state,
715 60 : format("Earth Tube Undisturbed Ground Temperature {}", nodeNum),
716 : Constant::Units::C,
717 30 : thisEarthTube.tUndist[nodeNum - 1],
718 : OutputProcessor::TimeStepType::Zone,
719 : OutputProcessor::StoreType::Average,
720 30 : zone.Name);
721 : }
722 4 : SetupOutputVariable(state,
723 : "Earth Tube Upper Boundary Ground Temperature",
724 : Constant::Units::C,
725 2 : thisEarthTube.tUpperBound,
726 : OutputProcessor::TimeStepType::Zone,
727 : OutputProcessor::StoreType::Average,
728 2 : zone.Name);
729 4 : SetupOutputVariable(state,
730 : "Earth Tube Lower Boundary Ground Temperature",
731 : Constant::Units::C,
732 2 : thisEarthTube.tLowerBound,
733 : OutputProcessor::TimeStepType::Zone,
734 : OutputProcessor::StoreType::Average,
735 2 : zone.Name);
736 : }
737 : } // ...end of firstTimeInits block
738 :
739 : Real64 timeElapsedLoc =
740 4296 : state.dataGlobal->HourOfDay + state.dataGlobal->TimeStep * state.dataGlobal->TimeStepZone + state.dataHVACGlobal->SysTimeElapsed;
741 4296 : if (state.dataEarthTube->timeElapsed !=
742 : timeElapsedLoc) { // time changed, update last with "current", avoids duplicate initializations and improper updates
743 4171 : if (state.dataGlobal->BeginDayFlag || state.dataGlobal->BeginEnvrnFlag) {
744 : // update all of the undisturbed temperatures (only need to do this once per day because the equation only changes as the day changes
745 228 : for (int etNum = 1; etNum <= totEarthTube; ++etNum) {
746 171 : auto &thisEarthTube = state.dataEarthTube->EarthTubeSys(etNum);
747 171 : if (thisEarthTube.ModelType != EarthTubeModelType::Vertical) {
748 113 : continue; // Skip earth tubes that do not use vertical solution
749 : }
750 58 : thisEarthTube.tUpperBound = thisEarthTube.calcUndisturbedGroundTemperature(state, thisEarthTube.depthUpperBound);
751 58 : thisEarthTube.tLowerBound = thisEarthTube.calcUndisturbedGroundTemperature(state, thisEarthTube.depthLowerBound);
752 928 : for (int nodeNum = 0; nodeNum <= thisEarthTube.totNodes - 1; ++nodeNum) {
753 870 : thisEarthTube.tUndist[nodeNum] = thisEarthTube.calcUndisturbedGroundTemperature(state, thisEarthTube.depthNode[nodeNum]);
754 : }
755 : }
756 : } // ...end of BeginDayFlag block
757 :
758 8321 : if (state.dataGlobal->BeginEnvrnFlag ||
759 4150 : (!state.dataGlobal->WarmupFlag && state.dataGlobal->BeginDayFlag && state.dataGlobal->DayOfSim == 1)) {
760 108 : for (int etNum = 1; etNum <= totEarthTube; ++etNum) {
761 81 : auto &thisEarthTube = state.dataEarthTube->EarthTubeSys(etNum);
762 81 : if (thisEarthTube.ModelType != EarthTubeModelType::Vertical) {
763 53 : continue; // Skip earth tubes that do not use vertical solution
764 : }
765 448 : for (int nodeNum = 0; nodeNum <= thisEarthTube.totNodes - 1; ++nodeNum) {
766 420 : thisEarthTube.tLast[nodeNum] = thisEarthTube.tUndist[nodeNum];
767 420 : thisEarthTube.tCurrent[nodeNum] = thisEarthTube.tLast[nodeNum];
768 : }
769 : }
770 : }
771 :
772 16684 : for (int etNum = 1; etNum <= totEarthTube; ++etNum) {
773 12513 : auto &thisEarthTube = state.dataEarthTube->EarthTubeSys(etNum);
774 12513 : if (thisEarthTube.ModelType != EarthTubeModelType::Vertical) {
775 8299 : continue; // Skip earth tubes that do not use vertical solution
776 : }
777 67424 : for (int nodeNum = 0; nodeNum <= thisEarthTube.totNodes - 1; ++nodeNum) {
778 63210 : thisEarthTube.tLast[nodeNum] = thisEarthTube.tCurrent[nodeNum];
779 : }
780 : }
781 : }
782 4296 : state.dataEarthTube->timeElapsed = timeElapsedLoc;
783 4296 : }
784 :
785 2 : void EarthTubeData::initCPrime0()
786 : {
787 : // Calculate c' for when effectiveness is zero. Will use these values when there is no air flow through the earth tube
788 : // and also use the values in the top portion of the solution (before the earth tube node) since these will not change.
789 2 : this->cPrime0[0] = this->cCoeff0[0] / this->bCoeff[0];
790 28 : for (int i = 1; i <= this->totNodes - 2; ++i) {
791 26 : this->cPrime0[i] = this->cCoeff0[i] / (this->bCoeff[i] - this->aCoeff[i] * this->cPrime0[i - 1]);
792 : }
793 2 : cPrime0[this->totNodes - 1] = 0.0;
794 2 : }
795 :
796 4296 : void CalcEarthTube(EnergyPlusData &state)
797 : {
798 :
799 : // SUBROUTINE INFORMATION:
800 : // AUTHOR Kwang Ho Lee
801 : // DATE WRITTEN November 2005
802 :
803 : // PURPOSE OF THIS SUBROUTINE:
804 : // This subroutine simulates the components making up the EarthTube unit.
805 :
806 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
807 : Real64 Process1; // Variable Used in the Middle of the Calculation
808 : Real64 GroundTempt; // Ground Temperature between Depth z at time t
809 :
810 : Real64 AirThermCond; // Thermal Conductivity of Air (W/mC)
811 : Real64 AirKinemVisco; // Kinematic Viscosity of Air (m2/s)
812 : Real64 AirThermDiffus; // Thermal Diffusivity of Air (m2/s)
813 : Real64 Re; // Reynolds Number for Flow Inside Pipe
814 : Real64 Pr; // Prandtl Number for Flow Inside Pipe
815 : Real64 Nu; // Nusselt Number for Flow Inside Pipe
816 : Real64 fa; // Friction Factor of Pipe
817 : Real64 PipeHeatTransCoef; // Convective Heat Transfer Coefficient at Inner Pipe Surface
818 : Real64 Rc; // Thermal Resistance due to Convection between Air and Pipe Inner Surface
819 : Real64 Rp; // Thermal Resistance due to Conduction between Pipe Inner and Outer Surface
820 : Real64 Rs; // Thermal Resistance due to Conduction between Pipe Outer Surface and Soil
821 : Real64 Rt; // Total Thermal Resistance between Pipe Air and Soil
822 : Real64 OverallHeatTransCoef; // Overall Heat Transfer Coefficient of Earth Tube
823 : Real64 AverPipeAirVel; // Average Pipe Air Velocity (m/s)
824 : Real64 AirMassFlowRate; // Actual Mass Flow Rate of Air inside Pipe
825 : Real64 AirSpecHeat; // Specific Heat of Air
826 : Real64 AirDensity; // Density of Air
827 : Real64 EVF;
828 :
829 4296 : int numEarthTubes = (int)state.dataEarthTube->EarthTubeSys.size();
830 4296 : Real64 outTdb = state.dataEnvrn->OutDryBulbTemp;
831 17184 : for (int Loop = 1; Loop <= numEarthTubes; ++Loop) {
832 12888 : auto &thisEarthTube = state.dataEarthTube->EarthTubeSys(Loop);
833 12888 : int NZ = thisEarthTube.ZonePtr;
834 12888 : auto &thisZoneHB = state.dataZoneTempPredictorCorrector->zoneHeatBalance(NZ);
835 12888 : thisZoneHB.MCPTE = 0.0;
836 12888 : thisZoneHB.MCPE = 0.0;
837 12888 : thisZoneHB.EAMFL = 0.0;
838 12888 : thisZoneHB.EAMFLxHumRat = 0.0;
839 12888 : thisEarthTube.FanPower = 0.0;
840 :
841 : // Don't simulate for Basic Solution if the zone is below the minimum temperature limit, above the maximum temperature limit
842 : // or below the temperature difference limit
843 25776 : bool tempShutDown = thisZoneHB.MAT < thisEarthTube.MinTemperature || thisZoneHB.MAT > thisEarthTube.MaxTemperature ||
844 12888 : std::abs(thisZoneHB.MAT - outTdb) < thisEarthTube.DelTemperature;
845 : // check for Basic model and some temperature limit preventing the earth tube from running
846 12888 : if ((thisEarthTube.ModelType == EarthTubeModelType::Basic) && (tempShutDown)) {
847 2012 : continue;
848 : }
849 :
850 10876 : AirDensity = Psychrometrics::PsyRhoAirFnPbTdbW(state, state.dataEnvrn->OutBaroPress, outTdb, state.dataEnvrn->OutHumRat);
851 10876 : AirSpecHeat = Psychrometrics::PsyCpAirFnW(state.dataEnvrn->OutHumRat);
852 10876 : if (tempShutDown) {
853 281 : EVF = 0.0;
854 : } else {
855 10595 : EVF = thisEarthTube.DesignLevel * thisEarthTube.availSched->getCurrentVal();
856 : }
857 10876 : thisZoneHB.MCPE =
858 21752 : EVF * AirDensity * AirSpecHeat *
859 21752 : (thisEarthTube.ConstantTermCoef +
860 10876 : std::abs(outTdb - state.dataZoneTempPredictorCorrector->zoneHeatBalance(NZ).MAT) * thisEarthTube.TemperatureTermCoef +
861 10876 : state.dataEnvrn->WindSpeed * (thisEarthTube.VelocityTermCoef + state.dataEnvrn->WindSpeed * thisEarthTube.VelocitySQTermCoef));
862 :
863 10876 : thisZoneHB.EAMFL = thisZoneHB.MCPE / AirSpecHeat;
864 10876 : if (thisEarthTube.FanEfficiency > 0.0) {
865 10876 : thisEarthTube.FanPower = thisZoneHB.EAMFL * thisEarthTube.FanPressure / (thisEarthTube.FanEfficiency * AirDensity);
866 : }
867 :
868 10876 : AverPipeAirVel = EVF / Constant::Pi / pow_2(thisEarthTube.r1);
869 10876 : AirMassFlowRate = EVF * AirDensity;
870 :
871 10876 : if (thisEarthTube.ModelType == EarthTubeModelType::Basic) {
872 : // Calculation of Ground Temperature at Depth z at time t for Basic model
873 6496 : GroundTempt = thisEarthTube.calcUndisturbedGroundTemperature(state, thisEarthTube.z);
874 6496 : thisEarthTube.GroundTempt = GroundTempt;
875 : }
876 :
877 : // Calculation of Convective Heat Transfer Coefficient at Inner Pipe Surface
878 10876 : AirThermCond = 0.02442 + 0.6992 * outTdb / 10000.0;
879 10876 : AirKinemVisco = (0.1335 + 0.000925 * outTdb) / 10000.0;
880 10876 : AirThermDiffus = (0.0014 * outTdb + 0.1872) / 10000.0;
881 10876 : Re = 2.0 * thisEarthTube.r1 * AverPipeAirVel / AirKinemVisco;
882 10876 : Pr = AirKinemVisco / AirThermDiffus;
883 10876 : if (Re <= 2300.0) {
884 281 : Nu = 3.66;
885 10595 : } else if (Re <= 4000.0) {
886 0 : fa = std::pow(1.58 * std::log(Re) - 3.28, -2);
887 0 : Process1 = (fa / 2.0) * (Re - 1000.0) * Pr / (1.0 + 12.7 * std::sqrt(fa / 2.0) * (std::pow(Pr, 2.0 / 3.0) - 1.0));
888 0 : Nu = (Process1 - 3.66) / (1700.0) * Re + (4000.0 * 3.66 - 2300.0 * Process1) / 1700.0;
889 : } else {
890 10595 : fa = std::pow(1.58 * std::log(Re) - 3.28, -2);
891 10595 : Nu = (fa / 2.0) * (Re - 1000.0) * Pr / (1.0 + 12.7 * std::sqrt(fa / 2.0) * (std::pow(Pr, 2.0 / 3.0) - 1.0));
892 : }
893 10876 : PipeHeatTransCoef = Nu * AirThermCond / 2.0 / thisEarthTube.r1;
894 :
895 : // Calculation of Thermal Resistance and Overall Heat Transfer Coefficient
896 10876 : Rc = 1.0 / 2.0 / Constant::Pi / thisEarthTube.r1 / PipeHeatTransCoef;
897 10876 : Rp = std::log((thisEarthTube.r1 + thisEarthTube.r2) / thisEarthTube.r1) / 2.0 / Constant::Pi / thisEarthTube.PipeThermCond;
898 10876 : if (thisEarthTube.r3 > 0.0) {
899 6496 : Rs = std::log((thisEarthTube.r1 + thisEarthTube.r2 + thisEarthTube.r3) / (thisEarthTube.r1 + thisEarthTube.r2)) / 2.0 / Constant::Pi /
900 6496 : thisEarthTube.SoilThermCond;
901 : } else { // for the Vertical solution .r3 was reset to zero for this
902 4380 : Rs = 0.0;
903 : }
904 10876 : Rt = Rc + Rp + Rs;
905 10876 : OverallHeatTransCoef = 1.0 / Rt;
906 :
907 10876 : switch (thisEarthTube.ModelType) {
908 4380 : case EarthTubeModelType::Vertical: {
909 : // First calculate term that will need to be added at the diagonal for flow and then solve the matrix for new temperatures
910 : Real64 eff; // effectiveness
911 4380 : if (AirMassFlowRate > 0.0) {
912 : // Calculate the NTU parameter: NTU = UA/[(Mdot*Cp)min] where Mdot*Cp is for the air side
913 : // where: U = OverallHeatTransCoef
914 : // A = 2*Pi*r1*TubeLength
915 4099 : Real64 NTU =
916 4099 : OverallHeatTransCoef * 2.0 * Constant::Pi * thisEarthTube.r1 * thisEarthTube.PipeLength / (AirMassFlowRate * AirSpecHeat);
917 :
918 : // Effectiveness is 1 - e(-NTU)
919 4099 : Real64 constexpr maxExpPower(50.0); // Maximum power after which EXP argument would be zero for DP variables
920 4099 : if (NTU > maxExpPower) {
921 0 : eff = 1.0;
922 : } else {
923 4099 : eff = 1.0 - std::exp(-NTU);
924 : }
925 : } else { // if no flow, then eff is zero
926 281 : eff = 0.0;
927 : }
928 :
929 4380 : Real64 airFlowTerm = AirMassFlowRate * AirSpecHeat * eff * thisEarthTube.airFlowCoeff;
930 4380 : thisEarthTube.calcVerticalEarthTube(state, airFlowTerm);
931 :
932 4380 : int nodeET = state.dataEarthTube->EarthTubePars(thisEarthTube.vertParametersPtr).numNodesAbove;
933 4380 : if (eff <= 0.0) { // no flow--air temperature leaving earth tube is the same as what went in
934 281 : thisEarthTube.InsideAirTemp = outTdb;
935 4099 : } else if (eff >= 1.0) { // effectiveness is one so leaving temperature is the same as the ground node temperatre
936 0 : thisEarthTube.InsideAirTemp = thisEarthTube.tCurrent[nodeET];
937 : } else { // the temperature is between the inlet and ground temperatures
938 4099 : thisEarthTube.InsideAirTemp = outTdb - eff * (outTdb - thisEarthTube.tCurrent[nodeET]);
939 : }
940 :
941 4380 : } break;
942 6496 : case EarthTubeModelType::Basic: { // Basic model
943 6496 : if (AirMassFlowRate * AirSpecHeat == 0.0) {
944 0 : thisEarthTube.InsideAirTemp = GroundTempt;
945 :
946 : } else {
947 :
948 : // Calculation of Pipe Outlet Air Temperature
949 6496 : if (outTdb > GroundTempt) {
950 6484 : Process1 =
951 6484 : (std::log(std::abs(outTdb - GroundTempt)) * AirMassFlowRate * AirSpecHeat - OverallHeatTransCoef * thisEarthTube.PipeLength) /
952 6484 : (AirMassFlowRate * AirSpecHeat);
953 6484 : thisEarthTube.InsideAirTemp = std::exp(Process1) + GroundTempt;
954 12 : } else if (outTdb == GroundTempt) {
955 0 : thisEarthTube.InsideAirTemp = GroundTempt;
956 : } else {
957 12 : Process1 =
958 12 : (std::log(std::abs(outTdb - GroundTempt)) * AirMassFlowRate * AirSpecHeat - OverallHeatTransCoef * thisEarthTube.PipeLength) /
959 12 : (AirMassFlowRate * AirSpecHeat);
960 12 : thisEarthTube.InsideAirTemp = GroundTempt - std::exp(Process1);
961 : }
962 : }
963 6496 : } break;
964 0 : default: { // should never get here
965 0 : assert(false);
966 : } break;
967 : }
968 :
969 10876 : thisEarthTube.CalcEarthTubeHumRat(state, NZ);
970 : }
971 4296 : }
972 :
973 7482 : Real64 EarthTubeData::calcUndisturbedGroundTemperature(EnergyPlusData &state, Real64 depth)
974 : {
975 7482 : return this->AverSoilSurTemp -
976 14964 : this->ApmlSoilSurTemp * std::exp(-depth * std::sqrt(Constant::Pi / 365.0 / this->SoilThermDiff)) *
977 7482 : std::cos(2.0 * Constant::Pi / 365.0 *
978 7482 : (state.dataEnvrn->DayOfYear - this->SoilSurPhaseConst - depth / 2.0 * std::sqrt(365.0 / Constant::Pi / this->SoilThermDiff)));
979 : }
980 :
981 4380 : void EarthTubeData::calcVerticalEarthTube(EnergyPlusData &state, Real64 airFlowTerm)
982 : {
983 : // Perform matrix calculations to model the earth tube using the vertical solution.
984 : // At this point, temperatures have already been shifted so tLast is correct and
985 : // undisturbed ground temperature have also been calculated. We need to assign/update
986 : // vectors of coefficients and then perform the Thomas algorithm.
987 : // Note that airFlowTerm is mdot_a*cp_a*eff*deltat/rho_soil/cp_soil/nodethickness_et/width/length
988 :
989 4380 : int nodeET = state.dataEarthTube->EarthTubePars(this->vertParametersPtr).numNodesAbove;
990 4380 : int nodeLast = this->totNodes - 1; // minus one because c++ arrays start at 0
991 :
992 : // First, calculate cPrime in the forward sweep.
993 : // If airFlowTerm is zero, there is no flow so we can use can use cPrime0 for cPrime.
994 4380 : if (airFlowTerm <= 0.0) {
995 4442 : for (int nodeNum = 0; nodeNum <= nodeLast; ++nodeNum) {
996 4161 : this->cPrime[nodeNum] = this->cPrime0[nodeNum];
997 : }
998 : } else { // there is positive flow so calculate cPrime
999 4099 : this->cPrime[0] = this->cCoeff[0] / this->bCoeff[0];
1000 61539 : for (int nodeNum = 1; nodeNum <= nodeLast; ++nodeNum) {
1001 57440 : Real64 addTerm = 0.0;
1002 57440 : if (nodeNum == nodeET) {
1003 4099 : addTerm = airFlowTerm;
1004 : }
1005 57440 : this->cPrime[nodeNum] = this->cCoeff[nodeNum] / (this->bCoeff[nodeNum] + addTerm - this->aCoeff[nodeNum] * this->cPrime[nodeNum - 1]);
1006 : }
1007 : }
1008 :
1009 : // Second, set-up dCoeff
1010 4380 : this->dCoeff[0] = this->tLast[0] + this->dMult0 * this->tUpperBound;
1011 61320 : for (int nodeNum = 1; nodeNum <= nodeLast - 1; ++nodeNum) {
1012 56940 : if (nodeNum != nodeET) {
1013 52560 : this->dCoeff[nodeNum] = this->tLast[nodeNum];
1014 : } else {
1015 4380 : this->dCoeff[nodeNum] = this->tLast[nodeNum] + airFlowTerm * state.dataEnvrn->OutDryBulbTemp;
1016 : }
1017 : }
1018 4380 : this->dCoeff[nodeLast] = this->tLast[nodeLast] + this->dMultN * this->tLowerBound;
1019 :
1020 : // Third, calculate dPrime in the forward sweep.
1021 4380 : this->dPrime[0] = this->dCoeff[0] / this->bCoeff[0];
1022 65700 : for (int nodeNum = 1; nodeNum <= nodeLast; ++nodeNum) {
1023 61320 : Real64 addTerm = 0.0;
1024 61320 : if (nodeNum == nodeET) {
1025 4380 : addTerm = airFlowTerm;
1026 : }
1027 122640 : this->dPrime[nodeNum] = (this->dCoeff[nodeNum] - this->aCoeff[nodeNum] * this->dPrime[nodeNum - 1]) /
1028 61320 : (this->bCoeff[nodeNum] + addTerm - this->aCoeff[nodeNum] * this->cPrime[nodeNum - 1]);
1029 : }
1030 :
1031 : // Finally, obtain the solution (tCurrent) by back substitution.
1032 4380 : this->tCurrent[nodeLast] = this->dPrime[nodeLast];
1033 65700 : for (int nodeNum = nodeLast - 1; nodeNum >= 0; --nodeNum) {
1034 61320 : this->tCurrent[nodeNum] = this->dPrime[nodeNum] - this->cPrime[nodeNum] * this->tCurrent[nodeNum + 1];
1035 : }
1036 4380 : }
1037 :
1038 10876 : void EarthTubeData::CalcEarthTubeHumRat(EnergyPlusData &state, int const NZ)
1039 : { // Zone number (index)
1040 :
1041 : // SUBROUTINE INFORMATION:
1042 : // AUTHOR Kwang Ho Lee
1043 : // DATE WRITTEN November 2005
1044 : // MODIFIED Rick Strand, June 2017 (made this a separate subroutine)
1045 :
1046 : // PURPOSE OF THIS SUBROUTINE:
1047 : // This subroutine determines the leaving humidity ratio for the EarthTube
1048 : // and calculates parameters associated with humidity ratio.
1049 :
1050 10876 : Real64 InsideDewPointTemp = Psychrometrics::PsyTdpFnWPb(state, state.dataEnvrn->OutHumRat, state.dataEnvrn->OutBaroPress);
1051 10876 : Real64 InsideHumRat = 0.0;
1052 10876 : auto &thisZoneHB = state.dataZoneTempPredictorCorrector->zoneHeatBalance(NZ);
1053 :
1054 10876 : if (this->InsideAirTemp >= InsideDewPointTemp) {
1055 1099 : InsideHumRat = state.dataEnvrn->OutHumRat;
1056 1099 : Real64 const InsideEnthalpy = Psychrometrics::PsyHFnTdbW(this->InsideAirTemp, state.dataEnvrn->OutHumRat);
1057 : // Intake fans will add some heat to the air, raising the temperature for an intake fan...
1058 1099 : if (this->FanType == Ventilation::Intake) {
1059 : Real64 OutletAirEnthalpy;
1060 337 : if (thisZoneHB.EAMFL == 0.0) {
1061 136 : OutletAirEnthalpy = InsideEnthalpy;
1062 : } else {
1063 201 : OutletAirEnthalpy = InsideEnthalpy + this->FanPower / thisZoneHB.EAMFL;
1064 : }
1065 337 : this->AirTemp = Psychrometrics::PsyTdbFnHW(OutletAirEnthalpy, state.dataEnvrn->OutHumRat);
1066 : } else {
1067 762 : this->AirTemp = this->InsideAirTemp;
1068 : }
1069 1099 : thisZoneHB.MCPTE = thisZoneHB.MCPE * this->AirTemp;
1070 :
1071 : } else {
1072 9777 : InsideHumRat = Psychrometrics::PsyWFnTdpPb(state, this->InsideAirTemp, state.dataEnvrn->OutBaroPress);
1073 9777 : Real64 const InsideEnthalpy = Psychrometrics::PsyHFnTdbW(this->InsideAirTemp, InsideHumRat);
1074 : // Intake fans will add some heat to the air, raising the temperature for an intake fan...
1075 9777 : if (this->FanType == Ventilation::Intake) {
1076 : Real64 OutletAirEnthalpy;
1077 3465 : if (thisZoneHB.EAMFL == 0.0) {
1078 0 : OutletAirEnthalpy = InsideEnthalpy;
1079 : } else {
1080 3465 : OutletAirEnthalpy = InsideEnthalpy + this->FanPower / thisZoneHB.EAMFL;
1081 : }
1082 3465 : this->AirTemp = Psychrometrics::PsyTdbFnHW(OutletAirEnthalpy, InsideHumRat);
1083 : } else {
1084 6312 : this->AirTemp = this->InsideAirTemp;
1085 : }
1086 9777 : thisZoneHB.MCPTE = thisZoneHB.MCPE * this->AirTemp;
1087 : }
1088 :
1089 10876 : this->HumRat = InsideHumRat;
1090 10876 : this->WetBulbTemp = Psychrometrics::PsyTwbFnTdbWPb(state, this->InsideAirTemp, InsideHumRat, state.dataEnvrn->OutBaroPress);
1091 10876 : thisZoneHB.EAMFLxHumRat = thisZoneHB.EAMFL * InsideHumRat;
1092 10876 : }
1093 :
1094 4296 : void ReportEarthTube(EnergyPlusData &state)
1095 : {
1096 :
1097 : // SUBROUTINE INFORMATION:
1098 : // AUTHOR Kwang Ho Lee
1099 : // DATE WRITTEN November 2005
1100 : // MODIFIED B. Griffith April 2010 added output reports
1101 :
1102 : // PURPOSE OF THIS SUBROUTINE: This subroutine fills remaining report variables.
1103 :
1104 4296 : Real64 const ReportingConstant = state.dataHVACGlobal->TimeStepSysSec;
1105 :
1106 17184 : for (int ZoneLoop = 1; ZoneLoop <= state.dataGlobal->NumOfZones; ++ZoneLoop) { // Start of zone loads report variable update loop ...
1107 12888 : auto &thisZone = state.dataEarthTube->ZnRptET(ZoneLoop);
1108 12888 : auto const &thisZoneHB = state.dataZoneTempPredictorCorrector->zoneHeatBalance(ZoneLoop);
1109 :
1110 : // Break the infiltration load into heat gain and loss components.
1111 : Real64 const AirDensity =
1112 12888 : Psychrometrics::PsyRhoAirFnPbTdbW(state, state.dataEnvrn->OutBaroPress, state.dataEnvrn->OutDryBulbTemp, state.dataEnvrn->OutHumRat);
1113 12888 : Real64 const CpAir = Psychrometrics::PsyCpAirFnW(state.dataEnvrn->OutHumRat);
1114 12888 : thisZone.EarthTubeVolume = (thisZoneHB.MCPE / CpAir / AirDensity) * ReportingConstant;
1115 12888 : thisZone.EarthTubeMass = (thisZoneHB.MCPE / CpAir) * ReportingConstant;
1116 12888 : thisZone.EarthTubeVolFlowRate = thisZoneHB.MCPE / CpAir / AirDensity;
1117 12888 : thisZone.EarthTubeVolFlowRateStd = thisZoneHB.MCPE / CpAir / state.dataEnvrn->StdRhoAir;
1118 12888 : thisZone.EarthTubeMassFlowRate = thisZoneHB.MCPE / CpAir;
1119 12888 : thisZone.EarthTubeWaterMassFlowRate = thisZoneHB.EAMFLxHumRat;
1120 :
1121 12888 : thisZone.EarthTubeFanElec = 0.0;
1122 12888 : thisZone.EarthTubeAirTemp = 0.0;
1123 25776 : for (auto const &thisEarthTube : state.dataEarthTube->EarthTubeSys) {
1124 25776 : if (thisEarthTube.ZonePtr == ZoneLoop) {
1125 12888 : thisZone.EarthTubeFanElec = thisEarthTube.FanPower * ReportingConstant;
1126 12888 : thisZone.EarthTubeFanElecPower = thisEarthTube.FanPower;
1127 :
1128 : // Break the EarthTube load into heat gain and loss components.
1129 12888 : if (thisZoneHB.ZT > thisEarthTube.AirTemp) {
1130 12766 : thisZone.EarthTubeHeatLoss = thisZoneHB.MCPE * (thisZoneHB.ZT - thisEarthTube.AirTemp) * ReportingConstant;
1131 12766 : thisZone.EarthTubeHeatLossRate = thisZoneHB.MCPE * (thisZoneHB.ZT - thisEarthTube.AirTemp);
1132 12766 : thisZone.EarthTubeHeatGain = 0.0;
1133 12766 : thisZone.EarthTubeHeatGainRate = 0.0;
1134 : } else {
1135 122 : thisZone.EarthTubeHeatGain = thisZoneHB.MCPE * (thisEarthTube.AirTemp - thisZoneHB.ZT) * ReportingConstant;
1136 122 : thisZone.EarthTubeHeatGainRate = thisZoneHB.MCPE * (thisEarthTube.AirTemp - thisZoneHB.ZT);
1137 122 : thisZone.EarthTubeHeatLoss = 0.0;
1138 122 : thisZone.EarthTubeHeatLossRate = 0.0;
1139 : }
1140 :
1141 12888 : thisZone.EarthTubeAirTemp = thisEarthTube.AirTemp;
1142 12888 : thisZone.EarthTubeWetBulbTemp = thisEarthTube.WetBulbTemp;
1143 12888 : thisZone.EarthTubeHumRat = thisEarthTube.HumRat;
1144 12888 : thisZone.EarthTubeOATreatmentPower = thisZoneHB.MCPE * (thisEarthTube.AirTemp - state.dataEnvrn->OutDryBulbTemp);
1145 12888 : break; // DO loop
1146 : }
1147 12888 : }
1148 :
1149 : } // ... end of zone loads report variable update loop.
1150 4296 : }
1151 :
1152 : } // namespace EnergyPlus::EarthTube
|
