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