Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2023, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : // C++ Headers
49 : #include <cmath>
50 : #include <vector>
51 :
52 : // ObjexxFCL Headers
53 : #include <ObjexxFCL/Array.functions.hh>
54 : #include <ObjexxFCL/Fmath.hh>
55 :
56 : // cpgfunction Headers
57 : #include <cpgfunction/boreholes.h>
58 : #include <cpgfunction/gfunction.h>
59 : #include <cpgfunction/segments.h>
60 :
61 : // JSON Headers
62 : #include <nlohmann/json.hpp>
63 :
64 : // EnergyPlus Headers
65 : #include <EnergyPlus/BranchNodeConnections.hh>
66 : #include <EnergyPlus/Data/EnergyPlusData.hh>
67 : #include <EnergyPlus/DataEnvironment.hh>
68 : #include <EnergyPlus/DataHVACGlobals.hh>
69 : #include <EnergyPlus/DataLoopNode.hh>
70 : #include <EnergyPlus/DataStringGlobals.hh>
71 : #include <EnergyPlus/DataSystemVariables.hh>
72 : #include <EnergyPlus/DisplayRoutines.hh>
73 : #include <EnergyPlus/FileSystem.hh>
74 : #include <EnergyPlus/FluidProperties.hh>
75 : #include <EnergyPlus/General.hh>
76 : #include <EnergyPlus/GroundHeatExchangers.hh>
77 : #include <EnergyPlus/GroundTemperatureModeling/GroundTemperatureModelManager.hh>
78 : #include <EnergyPlus/InputProcessing/InputProcessor.hh>
79 : #include <EnergyPlus/NodeInputManager.hh>
80 : #include <EnergyPlus/OutputProcessor.hh>
81 : #include <EnergyPlus/Plant/DataPlant.hh>
82 : #include <EnergyPlus/PlantUtilities.hh>
83 : #include <EnergyPlus/UtilityRoutines.hh>
84 :
85 : namespace EnergyPlus::GroundHeatExchangers {
86 : // MODULE INFORMATION:
87 : // AUTHOR Arun Murugappan, Dan Fisher
88 : // DATE WRITTEN September 2000
89 : // MODIFIED B. Griffith, Sept 2010,plant upgrades
90 : // Matt Mitchell, February 2015. Added Slinky GHX.
91 : // Moved models to object-oriented design.
92 : // RE-ENGINEERED na
93 :
94 : // PURPOSE OF THIS MODULE:
95 : // The module contains the data structures and routines to simulate the
96 : // operation of vertical closed-loop ground heat exchangers (GLHE) typically
97 : // used in low temperature geothermal heat pump systems.
98 :
99 : // METHODOLOGY EMPLOYED:
100 : // The borehole and fluid temperatures are calculated from the response to
101 : // the current heat transfer rate and the response to the history of past
102 : // applied heat pulses. The response to each pulse is calculated from a non-
103 : // dimensionalized response function, or G-function, that is specific to the
104 : // given borehole field arrangement, depth and spacing. The data defining
105 : // this function is read from input.
106 : // The heat pulse histories need to be recorded over an extended period (months).
107 : // To aid computational efficiency past pulses are continuously aggregated into
108 : // equivalent heat pulses of longer duration, as each pulse becomes less recent.
109 :
110 : // REFERENCES:
111 : // Eskilson, P. 'Thermal Analysis of Heat Extraction Boreholes' Ph.D. Thesis:
112 : // Dept. of Mathematical Physics, University of Lund, Sweden, June 1987.
113 : // Yavuzturk, C., J.D. Spitler. 1999. 'A Short Time Step Response Factor Model
114 : // for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2): 475-485.
115 : // Xiong, Z., D.E. Fisher, J.D. Spitler. 2015. 'Development and Validation of a Slinky
116 : // Ground Heat Exchanger.' Applied Energy. Vol 114, 57-69.
117 :
118 : // Using/Aliasing
119 : using namespace DataLoopNode;
120 :
121 : using namespace GroundTemperatureManager;
122 :
123 : // MODULE PARAMETER DEFINITIONS
124 : constexpr Real64 hrsPerMonth(730.0); // Number of hours in month
125 : constexpr Real64 maxTSinHr(60); // Max number of time step in a hour
126 : static constexpr std::array<std::string_view, 2> GFuncCalcMethodsStrs = {"UHFCALC", "UBHWTCALC"};
127 :
128 : //******************************************************************************
129 :
130 1 : GLHESlinky::GLHESlinky(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j)
131 : {
132 : // Check for duplicates
133 1 : for (auto &existingObj : state.dataGroundHeatExchanger->singleBoreholesVector) {
134 0 : if (objName == existingObj->name) {
135 0 : ShowFatalError(state, "Invalid input for " + this->moduleName + " object: Duplicate name found: " + existingObj->name);
136 : }
137 : }
138 :
139 1 : bool errorsFound = false;
140 :
141 1 : this->name = objName;
142 :
143 2 : std::string inletNodeName = UtilityRoutines::MakeUPPERCase(j["inlet_node_name"].get<std::string>());
144 2 : std::string outletNodeName = UtilityRoutines::MakeUPPERCase(j["outlet_node_name"].get<std::string>());
145 :
146 : // get inlet node num
147 1 : this->inletNodeNum = NodeInputManager::GetOnlySingleNode(state,
148 : inletNodeName,
149 : errorsFound,
150 : DataLoopNode::ConnectionObjectType::GroundHeatExchangerSlinky,
151 : this->name,
152 : DataLoopNode::NodeFluidType::Water,
153 : DataLoopNode::ConnectionType::Inlet,
154 : NodeInputManager::CompFluidStream::Primary,
155 1 : ObjectIsNotParent);
156 :
157 : // get outlet node num
158 1 : this->outletNodeNum = NodeInputManager::GetOnlySingleNode(state,
159 : outletNodeName,
160 : errorsFound,
161 : DataLoopNode::ConnectionObjectType::GroundHeatExchangerSlinky,
162 : this->name,
163 : DataLoopNode::NodeFluidType::Water,
164 : DataLoopNode::ConnectionType::Outlet,
165 : NodeInputManager::CompFluidStream::Primary,
166 1 : ObjectIsNotParent);
167 :
168 1 : this->available = true;
169 1 : this->on = true;
170 :
171 1 : BranchNodeConnections::TestCompSet(state, this->moduleName, this->name, inletNodeName, outletNodeName, "Condenser Water Nodes");
172 :
173 : // load data
174 1 : this->designFlow = j["design_flow_rate"].get<Real64>();
175 1 : PlantUtilities::RegisterPlantCompDesignFlow(state, this->inletNodeNum, this->designFlow);
176 :
177 1 : this->soil.k = j["soil_thermal_conductivity"].get<Real64>();
178 1 : this->soil.rho = j["soil_density"].get<Real64>();
179 1 : this->soil.cp = j["soil_specific_heat"].get<Real64>();
180 1 : this->soil.rhoCp = this->soil.rho * this->soil.cp;
181 1 : this->pipe.k = j["pipe_thermal_conductivity"].get<Real64>();
182 1 : this->pipe.rho = j["pipe_density"].get<Real64>();
183 1 : this->pipe.cp = j["pipe_specific_heat"].get<Real64>();
184 1 : this->pipe.outDia = j["pipe_outer_diameter"].get<Real64>();
185 1 : this->pipe.outRadius = this->pipe.outDia / 2.0;
186 1 : this->pipe.thickness = j["pipe_thickness"].get<Real64>();
187 :
188 2 : std::string const hxConfig = UtilityRoutines::MakeUPPERCase(j["heat_exchanger_configuration"].get<std::string>());
189 1 : if (UtilityRoutines::SameString(hxConfig, "VERTICAL")) {
190 1 : this->verticalConfig = true;
191 0 : } else if (UtilityRoutines::SameString(hxConfig, "HORIZONTAL")) {
192 0 : this->verticalConfig = false;
193 : }
194 :
195 1 : this->coilDiameter = j["coil_diameter"].get<Real64>();
196 1 : this->coilPitch = j["coil_pitch"].get<Real64>();
197 1 : this->trenchDepth = j["trench_depth"].get<Real64>();
198 1 : this->trenchLength = j["trench_length"].get<Real64>();
199 1 : this->numTrenches = j["number_of_trenches"].get<int>();
200 1 : this->trenchSpacing = j["horizontal_spacing_between_pipes"].get<Real64>();
201 1 : this->maxSimYears = j["maximum_length_of_simulation"].get<Real64>();
202 :
203 : // Need to add a response factor object for the slinky model
204 2 : std::shared_ptr<GLHEResponseFactors> thisRF(new GLHEResponseFactors);
205 1 : thisRF->name = "Response Factor Object Auto Generated No: " + fmt::to_string(state.dataGroundHeatExchanger->numAutoGeneratedResponseFactors + 1);
206 1 : this->myRespFactors = thisRF;
207 1 : state.dataGroundHeatExchanger->responseFactorsVector.push_back(thisRF);
208 :
209 : // Number of coils
210 1 : this->numCoils = static_cast<int>(this->trenchLength / this->coilPitch);
211 :
212 : // Total tube length
213 1 : this->totalTubeLength = DataGlobalConstants::Pi * this->coilDiameter * this->trenchLength * this->numTrenches / this->coilPitch;
214 :
215 : // Get g function data
216 1 : this->SubAGG = 15;
217 1 : this->AGG = 192;
218 :
219 : // Average coil depth
220 1 : if (this->verticalConfig) {
221 : // Vertical configuration
222 1 : if (this->trenchDepth - this->coilDiameter < 0.0) {
223 : // Error: part of the coil is above ground
224 0 : ShowSevereError(state, this->moduleName + "=\"" + this->name + "\", invalid value in field.");
225 0 : ShowContinueError(state, format("...{}=[{:.3R}].", "Trench Depth", this->trenchDepth));
226 0 : ShowContinueError(state, format("...{}=[{:.3R}].", "Coil Depth", this->coilDepth));
227 0 : ShowContinueError(state, "...Part of coil will be above ground.");
228 0 : errorsFound = true;
229 :
230 : } else {
231 : // Entire coil is below ground
232 1 : this->coilDepth = this->trenchDepth - (this->coilDiameter / 2.0);
233 : }
234 :
235 : } else {
236 : // Horizontal configuration
237 0 : this->coilDepth = this->trenchDepth;
238 : }
239 :
240 : // Thermal diffusivity of the ground
241 1 : this->soil.diffusivity = this->soil.k / this->soil.rhoCp;
242 :
243 1 : state.dataGroundHeatExchanger->prevTimeSteps.allocate(static_cast<int>((this->SubAGG + 1) * maxTSinHr + 1));
244 1 : state.dataGroundHeatExchanger->prevTimeSteps = 0.0;
245 :
246 1 : if (this->pipe.thickness >= this->pipe.outDia / 2.0) {
247 0 : ShowSevereError(state, this->moduleName + "=\"" + this->name + "\", invalid value in field.");
248 0 : ShowContinueError(state, format("...{}=[{:.3R}].", "Pipe Thickness", this->pipe.thickness));
249 0 : ShowContinueError(state, format("...{}=[{:.3R}].", "Pipe Outer Diameter", this->pipe.outDia));
250 0 : ShowContinueError(state, "...Radius will be <=0.");
251 0 : errorsFound = true;
252 : }
253 :
254 : // Initialize ground temperature model and get pointer reference
255 2 : std::string const gtmType = UtilityRoutines::MakeUPPERCase(j["undisturbed_ground_temperature_model_type"].get<std::string>());
256 2 : std::string const gtmName = UtilityRoutines::MakeUPPERCase(j["undisturbed_ground_temperature_model_name"].get<std::string>());
257 1 : this->groundTempModel = GetGroundTempModelAndInit(state, gtmType, gtmName);
258 :
259 : // Check for Errors
260 1 : if (errorsFound) {
261 0 : ShowFatalError(state, "Errors found in processing input for " + this->moduleName);
262 : }
263 1 : }
264 :
265 : //******************************************************************************
266 :
267 22 : GLHEVert::GLHEVert(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j)
268 : {
269 : // Check for duplicates
270 22 : for (auto &existingObj : state.dataGroundHeatExchanger->singleBoreholesVector) {
271 0 : if (objName == existingObj->name) {
272 0 : ShowFatalError(state, "Invalid input for " + this->moduleName + " object: Duplicate name found: " + existingObj->name);
273 : }
274 : }
275 :
276 22 : bool errorsFound = false;
277 :
278 22 : this->name = objName;
279 :
280 : // get inlet node num
281 44 : std::string const inletNodeName = UtilityRoutines::MakeUPPERCase(j["inlet_node_name"].get<std::string>());
282 22 : this->inletNodeNum = NodeInputManager::GetOnlySingleNode(state,
283 : inletNodeName,
284 : errorsFound,
285 : DataLoopNode::ConnectionObjectType::GroundHeatExchangerSystem,
286 : objName,
287 : DataLoopNode::NodeFluidType::Water,
288 : DataLoopNode::ConnectionType::Inlet,
289 : NodeInputManager::CompFluidStream::Primary,
290 22 : ObjectIsNotParent);
291 :
292 : // get outlet node num
293 44 : std::string const outletNodeName = UtilityRoutines::MakeUPPERCase(j["outlet_node_name"].get<std::string>());
294 22 : this->outletNodeNum = NodeInputManager::GetOnlySingleNode(state,
295 : outletNodeName,
296 : errorsFound,
297 : DataLoopNode::ConnectionObjectType::GroundHeatExchangerSystem,
298 : objName,
299 : DataLoopNode::NodeFluidType::Water,
300 : DataLoopNode::ConnectionType::Outlet,
301 : NodeInputManager::CompFluidStream::Primary,
302 22 : ObjectIsNotParent);
303 22 : this->available = true;
304 22 : this->on = true;
305 :
306 22 : BranchNodeConnections::TestCompSet(state, this->moduleName, objName, inletNodeName, outletNodeName, "Condenser Water Nodes");
307 :
308 22 : this->designFlow = j["design_flow_rate"].get<Real64>();
309 22 : PlantUtilities::RegisterPlantCompDesignFlow(state, this->inletNodeNum, this->designFlow);
310 :
311 22 : this->soil.k = j["ground_thermal_conductivity"].get<Real64>();
312 22 : this->soil.rhoCp = j["ground_thermal_heat_capacity"].get<Real64>();
313 :
314 22 : if (j.find("ghe_vertical_responsefactors_object_name") != j.end()) {
315 : // Response factors come from IDF object
316 20 : this->myRespFactors =
317 40 : GetResponseFactor(state, UtilityRoutines::MakeUPPERCase(j["ghe_vertical_responsefactors_object_name"].get<std::string>()));
318 20 : this->gFunctionsExist = true;
319 : }
320 :
321 : // no g-functions in the input file, so they need to be calculated
322 22 : if (!this->gFunctionsExist) {
323 :
324 : // g-function calculation method
325 2 : if (j.find("g_function_calculation_method") != j.end()) {
326 4 : std::string gFunctionMethodStr = UtilityRoutines::MakeUPPERCase(j["g_function_calculation_method"].get<std::string>());
327 2 : if (gFunctionMethodStr == "UHFCALC") {
328 1 : this->gFuncCalcMethod = GFuncCalcMethod::UniformHeatFlux;
329 1 : } else if (gFunctionMethodStr == "UBHWTCALC") {
330 1 : this->gFuncCalcMethod = GFuncCalcMethod::UniformBoreholeWallTemp;
331 : } else {
332 0 : errorsFound = true;
333 0 : ShowSevereError(state, fmt::format("g-Function Calculation Method: \"{}\" is invalid", gFunctionMethodStr));
334 : }
335 : }
336 :
337 : // get borehole data from array or individual borehole instance objects
338 2 : if (j.find("ghe_vertical_array_object_name") != j.end()) {
339 : // Response factors come from array object
340 4 : this->myRespFactors = BuildAndGetResponseFactorObjectFromArray(
341 6 : state, GetVertArray(state, UtilityRoutines::MakeUPPERCase(j["ghe_vertical_array_object_name"].get<std::string>())));
342 : } else {
343 0 : if (j.find("vertical_well_locations") == j.end()) {
344 : // No ResponseFactors, GHEArray, or SingleBH object are referenced
345 0 : ShowSevereError(state, "No GHE:ResponseFactors, GHE:Vertical:Array, or GHE:Vertical:Single objects found");
346 0 : ShowFatalError(state, "Check references to these objects for GHE:System object: " + this->name);
347 : }
348 :
349 0 : auto const vars = j.at("vertical_well_locations");
350 :
351 : // Calculate response factors from individual boreholes
352 0 : std::vector<std::shared_ptr<GLHEVertSingle>> tempVectOfBHObjects;
353 :
354 0 : for (auto const &var : vars) {
355 0 : if (!var.at("ghe_vertical_single_object_name").empty()) {
356 : std::shared_ptr<GLHEVertSingle> tempBHptr =
357 0 : GetSingleBH(state, UtilityRoutines::MakeUPPERCase(var.at("ghe_vertical_single_object_name").get<std::string>()));
358 0 : tempVectOfBHObjects.push_back(tempBHptr);
359 : } else {
360 0 : break;
361 : }
362 : }
363 :
364 0 : this->myRespFactors = BuildAndGetResponseFactorsObjectFromSingleBHs(state, tempVectOfBHObjects);
365 :
366 0 : if (!this->myRespFactors) {
367 0 : errorsFound = true;
368 0 : ShowSevereError(state, "GroundHeatExchanger:Vertical:Single objects not found.");
369 : }
370 : }
371 : }
372 :
373 22 : this->bhDiameter = this->myRespFactors->props->bhDiameter;
374 22 : this->bhRadius = this->bhDiameter / 2.0;
375 22 : this->bhLength = this->myRespFactors->props->bhLength;
376 22 : this->bhUTubeDist = this->myRespFactors->props->bhUTubeDist;
377 :
378 : // pull pipe and grout data up from response factor struct for simplicity
379 22 : this->pipe.outDia = this->myRespFactors->props->pipe.outDia;
380 22 : this->pipe.innerDia = this->myRespFactors->props->pipe.innerDia;
381 22 : this->pipe.outRadius = this->pipe.outDia / 2;
382 22 : this->pipe.innerRadius = this->pipe.innerDia / 2;
383 22 : this->pipe.thickness = this->myRespFactors->props->pipe.thickness;
384 22 : this->pipe.k = this->myRespFactors->props->pipe.k;
385 22 : this->pipe.rhoCp = this->myRespFactors->props->pipe.rhoCp;
386 :
387 22 : this->grout.k = this->myRespFactors->props->grout.k;
388 22 : this->grout.rhoCp = this->myRespFactors->props->grout.rhoCp;
389 :
390 22 : this->myRespFactors->gRefRatio = this->bhRadius / this->bhLength;
391 :
392 : // Number of simulation years from RunPeriod
393 22 : this->myRespFactors->maxSimYears = state.dataEnvrn->MaxNumberSimYears;
394 :
395 : // total tube length
396 22 : this->totalTubeLength = this->myRespFactors->numBoreholes * this->myRespFactors->props->bhLength;
397 :
398 : // ground thermal diffusivity
399 22 : this->soil.diffusivity = this->soil.k / this->soil.rhoCp;
400 :
401 : // multipole method constants
402 22 : this->theta_1 = this->bhUTubeDist / (2 * this->bhRadius);
403 22 : this->theta_2 = this->bhRadius / this->pipe.outRadius;
404 22 : this->theta_3 = 1 / (2 * this->theta_1 * this->theta_2);
405 22 : this->sigma = (this->grout.k - this->soil.k) / (this->grout.k + this->soil.k);
406 :
407 22 : this->SubAGG = 15;
408 22 : this->AGG = 192;
409 :
410 : // Allocation of all the dynamic arrays
411 22 : this->QnMonthlyAgg.dimension(static_cast<int>(this->myRespFactors->maxSimYears * 12), 0.0);
412 22 : this->QnHr.dimension(730 + this->AGG + this->SubAGG, 0.0);
413 22 : this->QnSubHr.dimension(static_cast<int>((this->SubAGG + 1) * maxTSinHr + 1), 0.0);
414 22 : this->LastHourN.dimension(this->SubAGG + 1, 0);
415 :
416 22 : state.dataGroundHeatExchanger->prevTimeSteps.allocate(static_cast<int>((this->SubAGG + 1) * maxTSinHr + 1));
417 22 : state.dataGroundHeatExchanger->prevTimeSteps = 0.0;
418 :
419 : // Initialize ground temperature model and get pointer reference
420 22 : this->groundTempModel =
421 44 : GetGroundTempModelAndInit(state,
422 44 : UtilityRoutines::MakeUPPERCase(j["undisturbed_ground_temperature_model_type"].get<std::string>()),
423 66 : UtilityRoutines::MakeUPPERCase(j["undisturbed_ground_temperature_model_name"].get<std::string>()));
424 :
425 : // Check for Errors
426 22 : if (errorsFound) {
427 0 : ShowFatalError(state, "Errors found in processing input for " + this->moduleName);
428 : }
429 22 : }
430 :
431 : //******************************************************************************
432 :
433 0 : GLHEVertSingle::GLHEVertSingle(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j)
434 : {
435 : // Check for duplicates
436 0 : for (auto &existingObj : state.dataGroundHeatExchanger->singleBoreholesVector) {
437 0 : if (objName == existingObj->name) {
438 0 : ShowFatalError(state, "Invalid input for " + this->moduleName + " object: Duplicate name found: " + existingObj->name);
439 : }
440 : }
441 :
442 0 : this->name = objName;
443 0 : this->props = GetVertProps(state, UtilityRoutines::MakeUPPERCase(j["ghe_vertical_properties_object_name"].get<std::string>()));
444 0 : this->xLoc = j["x_location"].get<Real64>();
445 0 : this->yLoc = j["y_location"].get<Real64>();
446 0 : this->dl_i = 0.0;
447 0 : this->dl_ii = 0.0;
448 0 : this->dl_j = 0.0;
449 0 : }
450 :
451 : //******************************************************************************
452 :
453 2 : GLHEVertArray::GLHEVertArray(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j)
454 : {
455 : // Check for duplicates
456 2 : for (auto &existingObj : state.dataGroundHeatExchanger->vertArraysVector) {
457 0 : if (objName == existingObj->name) {
458 0 : ShowFatalError(state, "Invalid input for " + this->moduleName + " object: Duplicate name found: " + existingObj->name);
459 : }
460 : }
461 :
462 2 : this->name = objName;
463 2 : this->props = GetVertProps(state, UtilityRoutines::MakeUPPERCase(j["ghe_vertical_properties_object_name"].get<std::string>()));
464 2 : this->numBHinXDirection = j["number_of_boreholes_in_x_direction"].get<int>();
465 2 : this->numBHinYDirection = j["number_of_boreholes_in_y_direction"].get<int>();
466 2 : this->bhSpacing = j["borehole_spacing"].get<Real64>();
467 2 : }
468 :
469 : //******************************************************************************
470 :
471 20 : GLHEResponseFactors::GLHEResponseFactors(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j)
472 : {
473 :
474 : // Check for duplicates
475 40 : for (auto &existingObj : state.dataGroundHeatExchanger->vertPropsVector) {
476 20 : if (objName == existingObj->name) {
477 0 : ShowFatalError(state, "Invalid input for " + this->moduleName + " object: Duplicate name found: " + existingObj->name);
478 : }
479 : }
480 :
481 20 : this->name = objName;
482 20 : this->props = GetVertProps(state, UtilityRoutines::MakeUPPERCase(j["ghe_vertical_properties_object_name"].get<std::string>()));
483 20 : this->numBoreholes = j["number_of_boreholes"].get<int>();
484 20 : this->gRefRatio = j["g_function_reference_ratio"].get<Real64>();
485 20 : this->maxSimYears = state.dataEnvrn->MaxNumberSimYears;
486 :
487 40 : auto const vars = j.at("g_functions");
488 40 : std::vector<Real64> tmpLntts;
489 40 : std::vector<Real64> tmpGvals;
490 720 : for (auto const &var : vars) {
491 700 : tmpLntts.push_back(var.at("g_function_ln_t_ts_value").get<Real64>());
492 700 : tmpGvals.push_back(var.at("g_function_g_value").get<Real64>());
493 : }
494 :
495 20 : bool errorsFound = false;
496 :
497 20 : if (tmpLntts.size() == tmpLntts.size()) {
498 20 : this->numGFuncPairs = static_cast<int>(tmpLntts.size());
499 : } else {
500 0 : errorsFound = true;
501 0 : ShowSevereError(state, "Errors found processing response factor input for Response Factor= " + this->name);
502 0 : ShowSevereError(state, "Uneven number of g-function pairs");
503 : }
504 :
505 20 : this->LNTTS.dimension(this->numGFuncPairs, 0.0);
506 20 : this->GFNC.dimension(this->numGFuncPairs, 0.0);
507 :
508 720 : for (size_t i = 1; i <= tmpLntts.size(); ++i) {
509 700 : this->LNTTS(i) = tmpLntts[i - 1];
510 700 : this->GFNC(i) = tmpGvals[i - 1];
511 : }
512 :
513 20 : if (errorsFound) {
514 0 : ShowFatalError(state, "Errors found in processing input for " + this->moduleName);
515 : }
516 20 : }
517 :
518 : //******************************************************************************
519 :
520 22 : GLHEVertProps::GLHEVertProps(EnergyPlusData &state, std::string const &objName, nlohmann::json const &j)
521 : {
522 :
523 : // Check for duplicates
524 22 : for (auto &existingObj : state.dataGroundHeatExchanger->vertPropsVector) {
525 0 : if (objName == existingObj->name) {
526 0 : ShowFatalError(state, "Invalid input for " + this->moduleName + " object: Duplicate name found: " + existingObj->name);
527 : }
528 : }
529 :
530 : // Load data from JSON
531 22 : this->name = objName;
532 22 : this->bhTopDepth = j["depth_of_top_of_borehole"].get<Real64>();
533 22 : this->bhLength = j["borehole_length"].get<Real64>();
534 22 : this->bhDiameter = j["borehole_diameter"].get<Real64>();
535 22 : this->grout.k = j["grout_thermal_conductivity"].get<Real64>();
536 22 : this->grout.rhoCp = j["grout_thermal_heat_capacity"].get<Real64>();
537 22 : this->pipe.k = j["pipe_thermal_conductivity"].get<Real64>();
538 22 : this->pipe.rhoCp = j["pipe_thermal_heat_capacity"].get<Real64>();
539 22 : this->pipe.outDia = j["pipe_outer_diameter"].get<Real64>();
540 22 : this->pipe.thickness = j["pipe_thickness"].get<Real64>();
541 22 : this->bhUTubeDist = j["u_tube_distance"].get<Real64>();
542 :
543 : // Verify u-tube spacing is valid
544 22 : if (this->bhUTubeDist < this->pipe.outDia) {
545 1 : ShowWarningError(state, "Borehole shank spacing is less than the pipe diameter. U-tube spacing is reference from the u-tube pipe center.");
546 1 : ShowWarningError(state, "Shank spacing is set to the outer pipe diameter.");
547 1 : this->bhUTubeDist = this->pipe.outDia;
548 : }
549 :
550 : // Set remaining data derived from previous inputs
551 22 : this->pipe.innerDia = this->pipe.outDia - 2 * this->pipe.thickness;
552 22 : this->pipe.outRadius = this->pipe.outDia / 2;
553 22 : this->pipe.innerRadius = this->pipe.innerDia / 2;
554 22 : }
555 :
556 : //******************************************************************************
557 :
558 30 : std::shared_ptr<GLHEVertProps> GetVertProps(EnergyPlusData &state, std::string const &objectName)
559 : {
560 : // Check if this instance of this model has already been retrieved
561 30 : for (auto &thisProp : state.dataGroundHeatExchanger->vertPropsVector) {
562 : // Check if the type and name match
563 30 : if (objectName == thisProp->name) {
564 30 : return thisProp;
565 : }
566 : }
567 :
568 0 : ShowSevereError(state, fmt::format("Object=GroundHeatExchanger:Vertical:Properties, Name={} - not found.", objectName));
569 0 : ShowFatalError(state, "Preceding errors cause program termination");
570 :
571 : // needed to silence compiler, but should never get here
572 0 : return nullptr;
573 : }
574 :
575 : //******************************************************************************
576 :
577 0 : std::shared_ptr<GLHEVertSingle> GetSingleBH(EnergyPlusData &state, std::string const &objectName)
578 : {
579 : // Check if this instance of this model has already been retrieved
580 0 : for (auto &thisBH : state.dataGroundHeatExchanger->singleBoreholesVector) {
581 : // Check if the type and name match
582 0 : if (objectName == thisBH->name) {
583 0 : return thisBH;
584 : }
585 : }
586 :
587 0 : ShowSevereError(state, fmt::format("Object=GroundHeatExchanger:Vertical:Single, Name={} - not found.", objectName));
588 0 : ShowFatalError(state, "Preceding errors cause program termination");
589 :
590 : // needed to silence compiler, but should never get here
591 0 : return nullptr;
592 : }
593 :
594 : //******************************************************************************
595 :
596 2 : std::shared_ptr<GLHEVertArray> GetVertArray(EnergyPlusData &state, std::string const &objectName)
597 : {
598 : // Check if this instance of this model has already been retrieved
599 2 : for (auto &thisProp : state.dataGroundHeatExchanger->vertArraysVector) {
600 : // Check if the type and name match
601 2 : if (objectName == thisProp->name) {
602 2 : return thisProp;
603 : }
604 : }
605 :
606 0 : ShowSevereError(state, fmt::format("Object=GroundHeatExchanger:Vertical:Array, Name={} - not found.", objectName));
607 0 : ShowFatalError(state, "Preceding errors cause program termination");
608 :
609 : // needed to silence compiler, but should never get here
610 0 : return nullptr;
611 : }
612 :
613 : //******************************************************************************
614 :
615 20 : std::shared_ptr<GLHEResponseFactors> GetResponseFactor(EnergyPlusData &state, std::string const &objectName)
616 : {
617 : // Check if this instance of this model has already been retrieved
618 20 : for (auto &thisRF : state.dataGroundHeatExchanger->responseFactorsVector) {
619 : // Check if the type and name match
620 20 : if (objectName == thisRF->name) {
621 20 : return thisRF;
622 : }
623 : }
624 :
625 0 : ShowSevereError(state, fmt::format("Object=GroundHeatExchanger:ResponseFactors, Name={} - not found.", objectName));
626 0 : ShowFatalError(state, "Preceding errors cause program termination");
627 :
628 : // needed to silence compiler, but should never get here
629 0 : return nullptr;
630 : }
631 :
632 : //******************************************************************************
633 :
634 2 : std::shared_ptr<GLHEResponseFactors> BuildAndGetResponseFactorObjectFromArray(EnergyPlusData &state,
635 : std::shared_ptr<GLHEVertArray> const &arrayObjectPtr)
636 : {
637 : // Make new response factor object and store it for later use
638 2 : std::shared_ptr<GLHEResponseFactors> thisRF(new GLHEResponseFactors);
639 2 : thisRF->name = arrayObjectPtr->name;
640 2 : thisRF->props = arrayObjectPtr->props;
641 :
642 : // Build out new instances of the vertical BH objects which correspond to this object
643 2 : int xLoc = 0;
644 2 : int bhCounter = 0;
645 6 : for (int xBH = 1; xBH <= arrayObjectPtr->numBHinXDirection; ++xBH) {
646 4 : int yLoc = 0;
647 12 : for (int yBH = 1; yBH <= arrayObjectPtr->numBHinYDirection; ++yBH) {
648 8 : bhCounter += 1;
649 16 : std::shared_ptr<GLHEVertSingle> thisBH(new GLHEVertSingle);
650 8 : thisBH->name = format("{} BH {} loc: ({}, {})", thisRF->name, bhCounter, xLoc, yLoc);
651 8 : thisBH->props = GetVertProps(state, arrayObjectPtr->props->name);
652 8 : thisBH->xLoc = xLoc;
653 8 : thisBH->yLoc = yLoc;
654 8 : thisRF->myBorholes.push_back(thisBH);
655 8 : state.dataGroundHeatExchanger->singleBoreholesVector.push_back(thisBH);
656 8 : yLoc += arrayObjectPtr->bhSpacing;
657 8 : thisRF->numBoreholes += 1;
658 : }
659 4 : xLoc += arrayObjectPtr->bhSpacing;
660 : }
661 :
662 2 : SetupBHPointsForResponseFactorsObject(thisRF);
663 2 : state.dataGroundHeatExchanger->responseFactorsVector.push_back(thisRF);
664 2 : return thisRF;
665 : }
666 :
667 : //******************************************************************************
668 :
669 : std::shared_ptr<GLHEResponseFactors>
670 0 : BuildAndGetResponseFactorsObjectFromSingleBHs(EnergyPlusData &state, std::vector<std::shared_ptr<GLHEVertSingle>> const &singleBHsForRFVect)
671 : {
672 : // Make new response factor object and store it for later use
673 0 : std::shared_ptr<GLHEResponseFactors> thisRF(new GLHEResponseFactors);
674 0 : thisRF->name = format("Response Factor Object Auto Generated No: {}", state.dataGroundHeatExchanger->numAutoGeneratedResponseFactors + 1);
675 :
676 : // Make new props object which has the mean values of the other props objects referenced by the individual BH objects
677 0 : std::shared_ptr<GLHEVertProps> thisProps(new GLHEVertProps);
678 0 : thisProps->name = format("Response Factor Auto Generated Mean Props No: {}", state.dataGroundHeatExchanger->numAutoGeneratedResponseFactors + 1);
679 0 : int numBH = singleBHsForRFVect.size();
680 0 : for (auto &thisBH : state.dataGroundHeatExchanger->singleBoreholesVector) {
681 0 : thisProps->bhDiameter += thisBH->props->bhDiameter;
682 0 : thisProps->bhLength += thisBH->props->bhLength;
683 0 : thisProps->bhTopDepth += thisBH->props->bhTopDepth;
684 0 : thisProps->bhUTubeDist += thisBH->props->bhUTubeDist;
685 :
686 0 : thisProps->grout.cp += thisBH->props->grout.cp;
687 0 : thisProps->grout.diffusivity += thisBH->props->grout.diffusivity;
688 0 : thisProps->grout.k += thisBH->props->grout.k;
689 0 : thisProps->grout.rho += thisBH->props->grout.rho;
690 0 : thisProps->grout.rhoCp += thisBH->props->grout.rhoCp;
691 :
692 0 : thisProps->pipe.cp += thisBH->props->pipe.cp;
693 0 : thisProps->pipe.diffusivity += thisBH->props->pipe.diffusivity;
694 0 : thisProps->pipe.k += thisBH->props->pipe.k;
695 0 : thisProps->pipe.rho += thisBH->props->pipe.rho;
696 0 : thisProps->pipe.rhoCp += thisBH->props->pipe.rhoCp;
697 :
698 0 : thisProps->pipe.outDia += thisBH->props->pipe.outDia;
699 0 : thisProps->pipe.thickness += thisBH->props->pipe.thickness;
700 :
701 0 : thisProps->pipe.innerDia += (thisBH->props->pipe.outDia - 2 * thisBH->props->pipe.thickness);
702 :
703 0 : thisRF->myBorholes.push_back(thisBH);
704 : }
705 :
706 : // normalize by number of bh
707 0 : thisProps->bhDiameter /= numBH;
708 0 : thisProps->bhLength /= numBH;
709 0 : thisProps->bhTopDepth /= numBH;
710 0 : thisProps->bhUTubeDist /= numBH;
711 :
712 0 : thisProps->grout.cp /= numBH;
713 0 : thisProps->grout.diffusivity /= numBH;
714 0 : thisProps->grout.k /= numBH;
715 0 : thisProps->grout.rho /= numBH;
716 0 : thisProps->grout.rhoCp /= numBH;
717 :
718 0 : thisProps->pipe.cp /= numBH;
719 0 : thisProps->pipe.diffusivity /= numBH;
720 0 : thisProps->pipe.k /= numBH;
721 0 : thisProps->pipe.rho /= numBH;
722 0 : thisProps->pipe.rhoCp /= numBH;
723 :
724 0 : thisProps->pipe.outDia /= numBH;
725 0 : thisProps->pipe.thickness /= numBH;
726 :
727 0 : thisProps->pipe.innerDia /= numBH;
728 :
729 0 : thisRF->props = thisProps;
730 0 : thisRF->numBoreholes = thisRF->myBorholes.size();
731 0 : state.dataGroundHeatExchanger->vertPropsVector.push_back(thisProps);
732 :
733 0 : SetupBHPointsForResponseFactorsObject(thisRF);
734 :
735 0 : state.dataGroundHeatExchanger->responseFactorsVector.push_back(thisRF);
736 :
737 0 : state.dataGroundHeatExchanger->numAutoGeneratedResponseFactors += 1;
738 :
739 0 : return thisRF;
740 : }
741 :
742 : //******************************************************************************
743 :
744 2 : void SetupBHPointsForResponseFactorsObject(std::shared_ptr<GLHEResponseFactors> &thisRF)
745 : {
746 10 : for (auto &thisBH : thisRF->myBorholes) {
747 :
748 : // Using Simpson's rule the number of points (n+1) must be odd, therefore an even number of panels is required
749 : // Starting from i = 0 to i <= NumPanels produces an odd number of points
750 8 : constexpr int numPanels_i = 50;
751 8 : constexpr int numPanels_ii = 50;
752 8 : constexpr int numPanels_j = 560;
753 :
754 8 : thisBH->dl_i = thisBH->props->bhLength / numPanels_i;
755 416 : for (int i = 0; i <= numPanels_i; ++i) {
756 408 : MyCartesian newPoint;
757 408 : newPoint.x = thisBH->xLoc;
758 408 : newPoint.y = thisBH->yLoc;
759 408 : newPoint.z = thisBH->props->bhTopDepth + (i * thisBH->dl_i);
760 408 : thisBH->pointLocations_i.push_back(newPoint);
761 : }
762 :
763 8 : thisBH->dl_ii = thisBH->props->bhLength / numPanels_ii;
764 416 : for (int i = 0; i <= numPanels_ii; ++i) {
765 408 : MyCartesian newPoint;
766 : // For case when bh is being compared to itself, shift points by 1 radius in the horizontal plane
767 408 : newPoint.x = thisBH->xLoc + (thisBH->props->bhDiameter / 2.0) / sqrt(2.0);
768 408 : newPoint.y = thisBH->yLoc + (thisBH->props->bhDiameter / 2.0) / (-sqrt(2.0));
769 408 : newPoint.z = thisBH->props->bhTopDepth + (i * thisBH->dl_ii);
770 408 : thisBH->pointLocations_ii.push_back(newPoint);
771 : }
772 :
773 8 : thisBH->dl_j = thisBH->props->bhLength / numPanels_j;
774 4496 : for (int i = 0; i <= numPanels_j; ++i) {
775 4488 : MyCartesian newPoint;
776 4488 : newPoint.x = thisBH->xLoc;
777 4488 : newPoint.y = thisBH->yLoc;
778 4488 : newPoint.z = thisBH->props->bhTopDepth + (i * thisBH->dl_j);
779 4488 : thisBH->pointLocations_j.push_back(newPoint);
780 : }
781 : }
782 2 : }
783 :
784 : //******************************************************************************
785 :
786 115 : void GLHEBase::onInitLoopEquip(EnergyPlusData &state, [[maybe_unused]] const PlantLocation &calledFromLocation)
787 : {
788 115 : this->initGLHESimVars(state);
789 115 : }
790 :
791 : //******************************************************************************
792 :
793 398751 : void GLHEBase::simulate(EnergyPlusData &state,
794 : [[maybe_unused]] const PlantLocation &calledFromLocation,
795 : [[maybe_unused]] bool const FirstHVACIteration,
796 : [[maybe_unused]] Real64 &CurLoad,
797 : [[maybe_unused]] bool const RunFlag)
798 : {
799 :
800 398751 : if (this->needToSetupOutputVars) {
801 23 : this->setupOutput(state);
802 23 : this->needToSetupOutputVars = false;
803 : }
804 :
805 398751 : if (state.dataGlobal->KickOffSimulation) {
806 2023 : this->initGLHESimVars(state);
807 : } else {
808 396728 : this->initGLHESimVars(state);
809 396728 : this->calcGroundHeatExchanger(state);
810 396728 : this->updateGHX(state);
811 : }
812 398751 : }
813 :
814 : //******************************************************************************
815 :
816 23 : PlantComponent *GLHEBase::factory(EnergyPlusData &state, DataPlant::PlantEquipmentType objectType, std::string const &objectName)
817 : {
818 23 : if (state.dataGroundHeatExchanger->GetInput) {
819 23 : GetGroundHeatExchangerInput(state);
820 23 : state.dataGroundHeatExchanger->GetInput = false;
821 : }
822 23 : if (objectType == DataPlant::PlantEquipmentType::GrndHtExchgSystem) {
823 22 : for (auto &ghx : state.dataGroundHeatExchanger->verticalGLHE) {
824 22 : if (ghx.name == objectName) {
825 22 : return &ghx;
826 : }
827 : }
828 1 : } else if (objectType == DataPlant::PlantEquipmentType::GrndHtExchgSlinky) {
829 1 : for (auto &ghx : state.dataGroundHeatExchanger->slinkyGLHE) {
830 1 : if (ghx.name == objectName) {
831 1 : return &ghx;
832 : }
833 : }
834 : }
835 :
836 : // If we didn't find it, fatal
837 0 : ShowFatalError(state, "Ground Heat Exchanger Factory: Error getting inputs for GHX named: " + objectName);
838 : // Shut up the compiler
839 0 : return nullptr;
840 : }
841 :
842 : //******************************************************************************
843 :
844 5035536 : std::vector<Real64> GLHEVert::distances(MyCartesian const &point_i, MyCartesian const &point_j)
845 : {
846 10071072 : std::vector<Real64> sumVals;
847 :
848 : // Calculate the distance between points
849 5035536 : sumVals.push_back(pow_2(point_i.x - point_j.x));
850 5035536 : sumVals.push_back(pow_2(point_i.y - point_j.y));
851 5035536 : sumVals.push_back(pow_2(point_i.z - point_j.z));
852 :
853 5035536 : Real64 sumTot = 0;
854 5035536 : std::vector<Real64> retVals;
855 20142144 : std::for_each(sumVals.begin(), sumVals.end(), [&](Real64 n) { sumTot += n; });
856 5035536 : retVals.push_back(std::sqrt(sumTot));
857 :
858 : // Calculate distance to mirror point
859 5035536 : sumVals.pop_back();
860 5035536 : sumVals.push_back(pow_2(point_i.z - (-point_j.z)));
861 :
862 5035536 : sumTot = 0;
863 20142144 : std::for_each(sumVals.begin(), sumVals.end(), [&](Real64 n) { sumTot += n; });
864 5035536 : retVals.push_back(std::sqrt(sumTot));
865 :
866 10071072 : return retVals;
867 : }
868 :
869 : //******************************************************************************
870 :
871 5035536 : Real64 GLHEVert::calcResponse(std::vector<Real64> const &dists, Real64 const currTime)
872 : {
873 5035536 : Real64 pointToPointResponse = erfc(dists[0] / (2 * sqrt(this->soil.diffusivity * currTime))) / dists[0];
874 5035536 : Real64 pointToReflectedResponse = erfc(dists[1] / (2 * sqrt(this->soil.diffusivity * currTime))) / dists[1];
875 :
876 5035536 : return pointToPointResponse - pointToReflectedResponse;
877 : }
878 :
879 : //******************************************************************************
880 :
881 8976 : Real64 GLHEVert::integral(MyCartesian const &point_i, std::shared_ptr<GLHEVertSingle> const &bh_j, Real64 const currTime)
882 : {
883 :
884 : // This code could be optimized in a number of ways.
885 : // The first, most simple way would be to precompute the distances from point i to point j, then store them for reuse.
886 : // The second, more intensive method would be to break the calcResponse calls out into four different parts.
887 : // The first point, last point, odd points, and even points. Then multiply the odd/even points by their respective coefficient for the
888 : // Simpson's method. After that, all points are summed together and divided by 3.
889 :
890 8976 : Real64 sum_f = 0;
891 8976 : int i = 0;
892 8976 : int const lastIndex_j = static_cast<int>(bh_j->pointLocations_j.size() - 1u);
893 5044512 : for (auto &point_j : bh_j->pointLocations_j) {
894 10071072 : std::vector<Real64> dists = distances(point_i, point_j);
895 5035536 : Real64 const f = calcResponse(dists, currTime);
896 :
897 : // Integrate using Simpson's
898 5035536 : if (i == 0 || i == lastIndex_j) {
899 17952 : sum_f += f;
900 5017584 : } else if (isEven(i)) {
901 2504304 : sum_f += 2 * f;
902 : } else {
903 2513280 : sum_f += 4 * f;
904 : }
905 :
906 5035536 : ++i;
907 : }
908 :
909 8976 : return (bh_j->dl_j / 3.0) * sum_f;
910 : }
911 :
912 : //******************************************************************************
913 :
914 176 : Real64 GLHEVert::doubleIntegral(std::shared_ptr<GLHEVertSingle> const &bh_i, std::shared_ptr<GLHEVertSingle> const &bh_j, Real64 const currTime)
915 : {
916 :
917 : // Similar optimizations as discussed above could happen here
918 :
919 176 : if (bh_i == bh_j) {
920 :
921 44 : Real64 sum_f = 0;
922 44 : int i = 0;
923 44 : int const lastIndex = static_cast<int>(bh_i->pointLocations_ii.size() - 1u);
924 2288 : for (auto &thisPoint : bh_i->pointLocations_ii) {
925 :
926 2244 : Real64 f = integral(thisPoint, bh_j, currTime);
927 :
928 : // Integrate using Simpson's
929 2244 : if (i == 0 || i == lastIndex) {
930 88 : sum_f += f;
931 2156 : } else if (isEven(i)) {
932 1056 : sum_f += 2 * f;
933 : } else {
934 1100 : sum_f += 4 * f;
935 : }
936 :
937 2244 : ++i;
938 : }
939 :
940 44 : return (bh_i->dl_ii / 3.0) * sum_f;
941 :
942 : } else {
943 :
944 132 : Real64 sum_f = 0;
945 132 : int i = 0;
946 132 : int const lastIndex = static_cast<int>(bh_i->pointLocations_i.size() - 1u);
947 6864 : for (auto &thisPoint : bh_i->pointLocations_i) {
948 :
949 6732 : Real64 f = integral(thisPoint, bh_j, currTime);
950 :
951 : // Integrate using Simpson's
952 6732 : if (i == 0 || i == lastIndex) {
953 264 : sum_f += f;
954 6468 : } else if (isEven(i)) {
955 3168 : sum_f += 2 * f;
956 : } else {
957 3300 : sum_f += 4 * f;
958 : }
959 :
960 6732 : ++i;
961 : }
962 :
963 132 : return (bh_i->dl_i / 3.0) * sum_f;
964 : }
965 : }
966 :
967 : //******************************************************************************
968 :
969 2 : void GLHEVert::calcLongTimestepGFunctions(EnergyPlusData &state)
970 : {
971 2 : switch (this->gFuncCalcMethod) {
972 1 : case GFuncCalcMethod::UniformHeatFlux:
973 1 : this->calcUniformHeatFluxGFunctions(state);
974 1 : break;
975 1 : case GFuncCalcMethod::UniformBoreholeWallTemp:
976 1 : this->calcUniformBHWallTempGFunctions(state);
977 1 : break;
978 0 : default:
979 0 : assert(false);
980 : }
981 2 : }
982 :
983 : //******************************************************************************
984 :
985 1 : void GLHEVert::calcUniformBHWallTempGFunctions(EnergyPlusData &state)
986 : {
987 : // construct boreholes vector
988 2 : std::vector<gt::boreholes::Borehole> boreholes;
989 5 : for (auto &bh : this->myRespFactors->myBorholes) {
990 4 : boreholes.emplace_back(bh->props->bhLength, bh->props->bhTopDepth, bh->props->bhDiameter / 2.0, bh->xLoc, bh->yLoc);
991 : }
992 :
993 : // convert time to a std::vector from an Array1D
994 2 : std::vector<double> time;
995 12 : for (auto &v : this->myRespFactors->time) {
996 11 : time.push_back(v);
997 : }
998 :
999 : // Obtain number of segments by adaptive discretization
1000 2 : gt::segments::adaptive adptDisc;
1001 1 : int nSegments = adptDisc.discretize(this->bhLength, this->totalTubeLength);
1002 :
1003 2 : this->myRespFactors->GFNC =
1004 3 : gt::gfunction::uniform_borehole_wall_temperature(boreholes, time, this->soil.diffusivity, nSegments, true, state.dataGlobal->numThread);
1005 1 : }
1006 :
1007 : //******************************************************************************
1008 :
1009 2 : void GLHEVert::calcGFunctions(EnergyPlusData &state)
1010 : {
1011 :
1012 : // No other choice than to calculate the g-functions here
1013 2 : this->setupTimeVectors();
1014 2 : this->calcShortTimestepGFunctions(state);
1015 2 : this->calcLongTimestepGFunctions(state);
1016 2 : this->combineShortAndLongTimestepGFunctions();
1017 :
1018 : // save data for later
1019 2 : if (state.files.outputControl.glhe && !state.dataSysVars->DisableGLHECaching) {
1020 2 : myCacheData["Response Factors"]["time"] = std::vector<Real64>(this->myRespFactors->time.begin(), this->myRespFactors->time.end());
1021 2 : myCacheData["Response Factors"]["LNTTS"] = std::vector<Real64>(this->myRespFactors->LNTTS.begin(), this->myRespFactors->LNTTS.end());
1022 2 : myCacheData["Response Factors"]["GFNC"] = std::vector<Real64>(this->myRespFactors->GFNC.begin(), this->myRespFactors->GFNC.end());
1023 2 : writeGLHECacheToFile(state);
1024 : }
1025 2 : }
1026 :
1027 : //******************************************************************************
1028 :
1029 2 : void GLHEVert::setupTimeVectors()
1030 : {
1031 :
1032 2 : constexpr int numDaysInYear(365);
1033 2 : constexpr Real64 lnttsStepSize = 0.5;
1034 :
1035 : // Minimum simulation time for which finite line source method is applicable
1036 2 : constexpr Real64 lntts_min_for_long_timestep = -8.5;
1037 :
1038 : // Time scale constant
1039 2 : Real64 const t_s = pow_2(this->bhLength) / (9 * this->soil.diffusivity);
1040 :
1041 : // Temporary vector for holding the LNTTS vals
1042 4 : std::vector<Real64> tempLNTTS;
1043 :
1044 2 : tempLNTTS.push_back(lntts_min_for_long_timestep);
1045 :
1046 : // Determine how many g-function pairs to generate based on user defined maximum simulation time
1047 : while (true) {
1048 22 : Real64 maxPossibleSimTime = exp(tempLNTTS.back()) * t_s;
1049 22 : if (maxPossibleSimTime <
1050 22 : this->myRespFactors->maxSimYears * numDaysInYear * DataGlobalConstants::HoursInDay * DataGlobalConstants::SecInHour) {
1051 20 : tempLNTTS.push_back(tempLNTTS.back() + lnttsStepSize);
1052 : } else {
1053 2 : break;
1054 : }
1055 20 : }
1056 :
1057 : // Setup the arrays
1058 2 : this->myRespFactors->time.dimension(tempLNTTS.size(), 0.0);
1059 2 : this->myRespFactors->LNTTS.dimension(tempLNTTS.size(), 0.0);
1060 2 : this->myRespFactors->GFNC.dimension(tempLNTTS.size(), 0.0);
1061 :
1062 2 : int index = 1;
1063 24 : for (auto &thisLNTTS : tempLNTTS) {
1064 22 : this->myRespFactors->time(index) = exp(thisLNTTS) * t_s;
1065 22 : this->myRespFactors->LNTTS(index) = thisLNTTS;
1066 22 : ++index;
1067 : }
1068 2 : }
1069 :
1070 : //******************************************************************************
1071 :
1072 1 : void GLHEVert::calcUniformHeatFluxGFunctions(EnergyPlusData &state)
1073 : {
1074 1 : DisplayString(state, "Initializing GroundHeatExchanger:System: " + this->name);
1075 :
1076 : // Calculate the g-functions
1077 12 : for (size_t lntts_index = 1; lntts_index <= this->myRespFactors->LNTTS.size(); ++lntts_index) {
1078 55 : for (auto &bh_i : this->myRespFactors->myBorholes) {
1079 44 : Real64 sum_T_ji = 0;
1080 220 : for (auto &bh_j : this->myRespFactors->myBorholes) {
1081 176 : sum_T_ji += doubleIntegral(bh_i, bh_j, this->myRespFactors->time(lntts_index));
1082 : }
1083 44 : this->myRespFactors->GFNC(lntts_index) += sum_T_ji;
1084 : }
1085 11 : this->myRespFactors->GFNC(lntts_index) /= (2 * this->totalTubeLength);
1086 :
1087 22 : std::stringstream ss;
1088 11 : ss << std::fixed << std::setprecision(1) << float(lntts_index) / this->myRespFactors->LNTTS.size() * 100;
1089 :
1090 11 : DisplayString(state, "...progress: " + ss.str() + "%");
1091 : }
1092 1 : }
1093 :
1094 : //******************************************************************************
1095 :
1096 2 : void GLHEVert::calcShortTimestepGFunctions(EnergyPlusData &state)
1097 : {
1098 : using FluidProperties::GetDensityGlycol;
1099 : using FluidProperties::GetSpecificHeatGlycol;
1100 :
1101 : // SUBROUTINE PARAMETER DEFINITIONS:
1102 2 : constexpr const char *RoutineName("calcShortTimestepGFunctions");
1103 :
1104 : enum class CellType
1105 : {
1106 : Invalid = -1,
1107 : FLUID,
1108 : CONVECTION,
1109 : PIPE,
1110 : GROUT,
1111 : SOIL,
1112 : Num
1113 : };
1114 :
1115 : struct Cell
1116 : {
1117 :
1118 : ~Cell() = default;
1119 :
1120 : CellType type;
1121 : Real64 radius_center;
1122 : Real64 radius_outer;
1123 : Real64 radius_inner;
1124 : Real64 thickness;
1125 : Real64 vol;
1126 : Real64 conductivity;
1127 : Real64 rhoCp;
1128 : Real64 temperature;
1129 : Real64 temperature_prev_ts;
1130 :
1131 1070 : Cell()
1132 1070 : : type(), radius_center(0.0), radius_outer(0.0), radius_inner(0.0), thickness(0.0), vol(0.0), conductivity(0.0), rhoCp(0.0),
1133 1070 : temperature(0.0), temperature_prev_ts(0.0)
1134 : {
1135 1070 : }
1136 : };
1137 :
1138 : // vector to hold 1-D cells
1139 4 : std::vector<Cell> Cells;
1140 :
1141 : // setup pipe, convection, and fluid layer geometries
1142 2 : constexpr int num_pipe_cells = 4;
1143 2 : constexpr int num_conv_cells = 1;
1144 2 : constexpr int num_fluid_cells = 3;
1145 2 : Real64 const pipe_thickness = this->pipe.thickness;
1146 2 : Real64 const pcf_cell_thickness = pipe_thickness / num_pipe_cells;
1147 2 : Real64 const radius_pipe_out = std::sqrt(2) * this->pipe.outRadius;
1148 2 : Real64 const radius_pipe_in = radius_pipe_out - pipe_thickness;
1149 2 : Real64 const radius_conv = radius_pipe_in - num_conv_cells * pcf_cell_thickness;
1150 2 : Real64 const radius_fluid = radius_conv - (num_fluid_cells - 0.5) * pcf_cell_thickness; // accounts for half thickness of boundary cell
1151 :
1152 : // setup grout layer geometry
1153 2 : constexpr int num_grout_cells = 27;
1154 2 : Real64 const radius_grout = this->bhRadius;
1155 2 : Real64 const grout_cell_thickness = (radius_grout - radius_pipe_out) / num_grout_cells;
1156 :
1157 : // setup soil layer geometry
1158 2 : constexpr int num_soil_cells = 500;
1159 2 : constexpr Real64 radius_soil = 10;
1160 2 : Real64 const soil_cell_thickness = (radius_soil - radius_grout) / num_soil_cells;
1161 :
1162 : // use design flow rate
1163 2 : this->massFlowRate = this->designMassFlow;
1164 :
1165 : // calculate equivalent thermal resistance between borehole wall and fluid
1166 2 : Real64 bhResistance = calcBHAverageResistance(state);
1167 2 : Real64 bhConvectionResistance = calcPipeConvectionResistance(state);
1168 2 : Real64 bh_equivalent_resistance_tube_grout = bhResistance - bhConvectionResistance / 2.0;
1169 2 : Real64 bh_equivalent_resistance_convection = bhResistance - bh_equivalent_resistance_tube_grout;
1170 :
1171 2 : Real64 initial_temperature = this->inletTemp;
1172 6 : Real64 cpFluid_init = GetSpecificHeatGlycol(state,
1173 2 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
1174 : initial_temperature,
1175 2 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
1176 4 : RoutineName);
1177 6 : Real64 fluidDensity_init = GetDensityGlycol(state,
1178 2 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
1179 : initial_temperature,
1180 2 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
1181 4 : RoutineName);
1182 :
1183 : // initialize the fluid cells
1184 8 : for (int i = 0; i < num_fluid_cells; ++i) {
1185 6 : Cell thisCell;
1186 6 : thisCell.type = CellType::FLUID;
1187 6 : thisCell.thickness = pcf_cell_thickness;
1188 6 : thisCell.radius_center = radius_fluid + i * thisCell.thickness;
1189 :
1190 : // boundary cell is only half thickness
1191 6 : if (i == 0) {
1192 2 : thisCell.radius_inner = thisCell.radius_center;
1193 : } else {
1194 4 : thisCell.radius_inner = thisCell.radius_center - thisCell.thickness / 2.0;
1195 : }
1196 :
1197 6 : thisCell.radius_outer = thisCell.radius_center + thisCell.thickness / 2.0;
1198 6 : thisCell.conductivity = 200;
1199 6 : thisCell.rhoCp = 2.0 * cpFluid_init * fluidDensity_init * pow_2(this->pipe.innerRadius) / (pow_2(radius_conv) - pow_2(radius_fluid));
1200 6 : Cells.push_back(thisCell);
1201 : }
1202 :
1203 : // initialize the convection cells
1204 4 : for (int i = 0; i < num_conv_cells; ++i) {
1205 2 : Cell thisCell;
1206 2 : thisCell.thickness = pcf_cell_thickness;
1207 2 : thisCell.radius_inner = radius_conv + i * thisCell.thickness;
1208 2 : thisCell.radius_center = thisCell.radius_inner + thisCell.thickness / 2.0;
1209 2 : thisCell.radius_outer = thisCell.radius_inner + thisCell.thickness;
1210 2 : thisCell.conductivity = log(radius_pipe_in / radius_conv) / (2 * DataGlobalConstants::Pi * bh_equivalent_resistance_convection);
1211 2 : thisCell.rhoCp = 1;
1212 2 : Cells.push_back(thisCell);
1213 : }
1214 :
1215 : // initialize pipe cells
1216 10 : for (int i = 0; i < num_pipe_cells; ++i) {
1217 8 : Cell thisCell;
1218 8 : thisCell.type = CellType::PIPE;
1219 8 : thisCell.thickness = pcf_cell_thickness;
1220 8 : thisCell.radius_inner = radius_pipe_in + i * thisCell.thickness;
1221 8 : thisCell.radius_center = thisCell.radius_inner + thisCell.thickness / 2.0;
1222 8 : thisCell.radius_outer = thisCell.radius_inner + thisCell.thickness;
1223 8 : thisCell.conductivity = log(radius_grout / radius_pipe_in) / (2 * DataGlobalConstants::Pi * bh_equivalent_resistance_tube_grout);
1224 8 : thisCell.rhoCp = this->pipe.rhoCp;
1225 8 : Cells.push_back(thisCell);
1226 : }
1227 :
1228 : // initialize grout cells
1229 56 : for (int i = 0; i < num_grout_cells; ++i) {
1230 54 : Cell thisCell;
1231 54 : thisCell.type = CellType::GROUT;
1232 54 : thisCell.thickness = grout_cell_thickness;
1233 54 : thisCell.radius_inner = radius_pipe_out + i * thisCell.thickness;
1234 54 : thisCell.radius_center = thisCell.radius_inner + thisCell.thickness / 2.0;
1235 54 : thisCell.radius_outer = thisCell.radius_inner + thisCell.thickness;
1236 54 : thisCell.conductivity = log(radius_grout / radius_pipe_in) / (2 * DataGlobalConstants::Pi * bh_equivalent_resistance_tube_grout);
1237 54 : thisCell.rhoCp = grout.rhoCp;
1238 54 : Cells.push_back(thisCell);
1239 : }
1240 :
1241 : // initialize soil cells
1242 1002 : for (int i = 0; i < num_soil_cells; ++i) {
1243 1000 : Cell thisCell;
1244 1000 : thisCell.type = CellType::SOIL;
1245 1000 : thisCell.thickness = soil_cell_thickness;
1246 1000 : thisCell.radius_inner = radius_grout + i * thisCell.thickness;
1247 1000 : thisCell.radius_center = thisCell.radius_inner + thisCell.thickness / 2.0;
1248 1000 : thisCell.radius_outer = thisCell.radius_inner + thisCell.thickness;
1249 1000 : thisCell.conductivity = this->soil.k;
1250 1000 : thisCell.rhoCp = this->soil.rhoCp;
1251 1000 : Cells.push_back(thisCell);
1252 : }
1253 :
1254 : // other non-geometric specific setup
1255 1072 : for (auto &thisCell : Cells) {
1256 1070 : thisCell.vol = DataGlobalConstants::Pi * (pow_2(thisCell.radius_outer) - pow_2(thisCell.radius_inner));
1257 1070 : thisCell.temperature = initial_temperature;
1258 : }
1259 :
1260 : // set upper limit of time for the short time-step g-function calcs so there is some overlap
1261 2 : Real64 constexpr lntts_max_for_short_timestep = -9.0;
1262 2 : Real64 const t_s = pow_2(this->bhLength) / (9.0 * this->soil.diffusivity);
1263 :
1264 2 : Real64 constexpr time_step = 500;
1265 2 : Real64 const time_max_for_short_timestep = exp(lntts_max_for_short_timestep) * t_s;
1266 2 : Real64 total_time = 0;
1267 :
1268 2 : Real64 constexpr heat_flux = 40.4;
1269 :
1270 : // time step loop
1271 1102 : while (total_time < time_max_for_short_timestep) {
1272 :
1273 294800 : for (auto &thisCell : Cells) {
1274 294250 : thisCell.temperature_prev_ts = thisCell.temperature;
1275 : }
1276 :
1277 1100 : std::vector<Real64> a;
1278 1100 : std::vector<Real64> b;
1279 1100 : std::vector<Real64> c;
1280 1100 : std::vector<Real64> d;
1281 :
1282 : // setup tdma matrices
1283 550 : int num_cells = Cells.size();
1284 294800 : for (int cell_index = 0; cell_index < num_cells; ++cell_index) {
1285 294250 : if (cell_index == 0) {
1286 : // heat flux BC
1287 :
1288 550 : auto &thisCell = Cells[cell_index];
1289 550 : auto &eastCell = Cells[cell_index + 1];
1290 :
1291 550 : Real64 FE1 = log(thisCell.radius_outer / thisCell.radius_center) / (2 * DataGlobalConstants::Pi * thisCell.conductivity);
1292 550 : Real64 FE2 = log(eastCell.radius_center / eastCell.radius_inner) / (2 * DataGlobalConstants::Pi * eastCell.conductivity);
1293 550 : Real64 AE = 1 / (FE1 + FE2);
1294 :
1295 550 : Real64 AD = thisCell.rhoCp * thisCell.vol / time_step;
1296 :
1297 550 : a.push_back(0);
1298 550 : b.push_back(-AE / AD - 1);
1299 550 : c.push_back(AE / AD);
1300 550 : d.push_back(-thisCell.temperature_prev_ts - heat_flux / AD);
1301 :
1302 293700 : } else if (cell_index == num_cells - 1) {
1303 : // const ground temp bc
1304 :
1305 550 : auto &thisCell = Cells[cell_index];
1306 :
1307 550 : a.push_back(0);
1308 550 : b.push_back(1);
1309 550 : c.push_back(0);
1310 550 : d.push_back(thisCell.temperature_prev_ts);
1311 :
1312 : } else {
1313 : // all other cells
1314 :
1315 293150 : auto &westCell = Cells[cell_index - 1];
1316 293150 : auto &thisCell = Cells[cell_index];
1317 293150 : auto &eastCell = Cells[cell_index + 1];
1318 :
1319 293150 : Real64 FE1 = log(thisCell.radius_outer / thisCell.radius_center) / (2 * DataGlobalConstants::Pi * thisCell.conductivity);
1320 293150 : Real64 FE2 = log(eastCell.radius_center / eastCell.radius_inner) / (2 * DataGlobalConstants::Pi * eastCell.conductivity);
1321 293150 : Real64 AE = 1 / (FE1 + FE2);
1322 :
1323 293150 : Real64 FW1 = log(westCell.radius_outer / westCell.radius_center) / (2 * DataGlobalConstants::Pi * westCell.conductivity);
1324 293150 : Real64 FW2 = log(thisCell.radius_center / thisCell.radius_inner) / (2 * DataGlobalConstants::Pi * thisCell.conductivity);
1325 293150 : Real64 AW = -1 / (FW1 + FW2);
1326 :
1327 293150 : Real64 AD = thisCell.rhoCp * thisCell.vol / time_step;
1328 :
1329 293150 : a.push_back(-AW / AD);
1330 293150 : b.push_back(AW / AD - AE / AD - 1);
1331 293150 : c.push_back(AE / AD);
1332 293150 : d.push_back(-thisCell.temperature_prev_ts);
1333 : }
1334 : } // end tdma setup
1335 :
1336 : // solve for new temperatures
1337 1100 : std::vector<Real64> new_temps = TDMA(a, b, c, d);
1338 :
1339 294800 : for (int cell_index = 0; cell_index < num_cells; ++cell_index) {
1340 294250 : Cells[cell_index].temperature = new_temps[cell_index];
1341 : }
1342 :
1343 : // calculate bh wall temp
1344 550 : Real64 T_bhWall = 0.0;
1345 19250 : for (int cell_index = 0; cell_index < num_cells; ++cell_index) {
1346 19250 : auto &leftCell = Cells[cell_index];
1347 19250 : auto &rightCell = Cells[cell_index + 1];
1348 :
1349 19250 : if (leftCell.type == CellType::GROUT && rightCell.type == CellType::SOIL) {
1350 :
1351 550 : Real64 left_conductance = 2 * DataGlobalConstants::Pi * leftCell.conductivity / log(leftCell.radius_outer / leftCell.radius_inner);
1352 : Real64 right_conductance =
1353 550 : 2 * DataGlobalConstants::Pi * rightCell.conductivity / log(rightCell.radius_center / leftCell.radius_inner);
1354 :
1355 550 : T_bhWall =
1356 550 : (left_conductance * leftCell.temperature + right_conductance * rightCell.temperature) / (left_conductance + right_conductance);
1357 550 : break;
1358 : }
1359 : }
1360 :
1361 550 : total_time += time_step;
1362 :
1363 1100 : GFNC_shortTimestep.push_back(2 * DataGlobalConstants::Pi * this->soil.k *
1364 550 : ((Cells[0].temperature - initial_temperature) / heat_flux - bhResistance));
1365 550 : LNTTS_shortTimestep.push_back(log(total_time / t_s));
1366 :
1367 : } // end timestep loop
1368 2 : }
1369 :
1370 : //******************************************************************************
1371 :
1372 550 : std::vector<Real64> TDMA(std::vector<Real64> a, std::vector<Real64> b, std::vector<Real64> c, std::vector<Real64> d)
1373 : {
1374 : // from: https://en.wikibooks.org/wiki/Algorithm_Implementation/Linear_Algebra/Tridiagonal_matrix_algorithm#C.2B.2B
1375 :
1376 550 : int n = static_cast<int>(d.size() - 1u);
1377 :
1378 550 : c[0] /= b[0];
1379 550 : d[0] /= b[0];
1380 :
1381 293700 : for (int i = 1; i < n; ++i) {
1382 293150 : c[i] /= b[i] - a[i] * c[i - 1];
1383 293150 : d[i] = (d[i] - a[i] * d[i - 1]) / (b[i] - a[i] * c[i - 1]);
1384 : }
1385 :
1386 550 : d[n] = (d[n] - a[n] * d[n - 1]) / (b[n] - a[n] * c[n - 1]);
1387 :
1388 294250 : for (int i = n; i-- > 0;) {
1389 293700 : d[i] -= c[i] * d[i + 1];
1390 : }
1391 :
1392 550 : return d;
1393 : }
1394 :
1395 : //******************************************************************************
1396 :
1397 2 : void GLHEVert::combineShortAndLongTimestepGFunctions()
1398 : {
1399 4 : std::vector<Real64> GFNC_combined;
1400 4 : std::vector<Real64> LNTTS_combined;
1401 :
1402 2 : Real64 const t_s = pow_2(this->bhLength) / (9.0 * this->soil.diffusivity);
1403 :
1404 : // Nothing to do. Just put the short time step g-functions on the combined vector
1405 2 : int num_shortTimestepGFunctions = GFNC_shortTimestep.size();
1406 552 : for (int index_shortTS = 0; index_shortTS < num_shortTimestepGFunctions; ++index_shortTS) {
1407 550 : GFNC_combined.push_back(GFNC_shortTimestep[index_shortTS]);
1408 550 : LNTTS_combined.push_back(LNTTS_shortTimestep[index_shortTS]);
1409 : }
1410 :
1411 : // Add the rest of the long time-step g-functions to the combined curve
1412 24 : for (int index_longTS = this->myRespFactors->GFNC.l(); index_longTS <= this->myRespFactors->GFNC.u(); ++index_longTS) {
1413 22 : GFNC_combined.push_back(this->myRespFactors->GFNC(index_longTS));
1414 22 : LNTTS_combined.push_back(this->myRespFactors->LNTTS(index_longTS));
1415 : }
1416 :
1417 : // Move combined values into right data struct
1418 2 : this->myRespFactors->time.deallocate();
1419 2 : this->myRespFactors->LNTTS.deallocate();
1420 2 : this->myRespFactors->GFNC.deallocate();
1421 :
1422 2 : this->myRespFactors->time.dimension(GFNC_combined.size(), 0.0);
1423 2 : this->myRespFactors->LNTTS.dimension(GFNC_combined.size(), 0.0);
1424 2 : this->myRespFactors->GFNC.dimension(GFNC_combined.size(), 0.0);
1425 :
1426 574 : for (unsigned int index = 0; index < GFNC_combined.size(); ++index) {
1427 572 : this->myRespFactors->time[index] = exp(LNTTS_combined[index]) * t_s;
1428 572 : this->myRespFactors->LNTTS[index] = LNTTS_combined[index];
1429 572 : this->myRespFactors->GFNC[index] = GFNC_combined[index];
1430 : }
1431 2 : }
1432 :
1433 3 : void GLHEBase::makeThisGLHECacheAndCompareWithFileCache(EnergyPlusData &state)
1434 : {
1435 3 : if (state.files.outputControl.glhe && !state.dataSysVars->DisableGLHECaching) {
1436 3 : makeThisGLHECacheStruct();
1437 3 : readCacheFileAndCompareWithThisGLHECache(state);
1438 : }
1439 3 : }
1440 :
1441 : //******************************************************************************
1442 :
1443 2 : void GLHEVert::makeThisGLHECacheStruct()
1444 : {
1445 : // For convenience
1446 2 : auto &d = myCacheData["Phys Data"];
1447 :
1448 2 : d["Flow Rate"] = this->designFlow;
1449 2 : d["Soil k"] = this->soil.k;
1450 2 : d["Soil rhoCp"] = this->soil.rhoCp;
1451 2 : d["BH Top Depth"] = this->myRespFactors->props->bhTopDepth;
1452 2 : d["BH Length"] = this->myRespFactors->props->bhLength;
1453 2 : d["BH Diameter"] = this->myRespFactors->props->bhDiameter;
1454 2 : d["Grout k"] = this->myRespFactors->props->grout.k;
1455 2 : d["Grout rhoCp"] = this->myRespFactors->props->grout.rhoCp;
1456 2 : d["Pipe k"] = this->myRespFactors->props->pipe.k;
1457 2 : d["Pipe rhoCP"] = this->myRespFactors->props->pipe.rhoCp;
1458 2 : d["Pipe Diameter"] = this->myRespFactors->props->pipe.outDia;
1459 2 : d["Pipe Thickness"] = this->myRespFactors->props->pipe.thickness;
1460 2 : d["U-tube Dist"] = this->myRespFactors->props->bhUTubeDist;
1461 2 : d["Max Simulation Years"] = this->myRespFactors->maxSimYears;
1462 2 : d["g-Function Calc Method"] = GroundHeatExchangers::GFuncCalcMethodsStrs[int(this->gFuncCalcMethod)];
1463 :
1464 2 : auto &d_bh_data = d["BH Data"];
1465 2 : int i = 0;
1466 10 : for (auto &thisBH : this->myRespFactors->myBorholes) {
1467 8 : ++i;
1468 8 : auto &d_bh = d_bh_data[fmt::format("BH {}", i)];
1469 8 : d_bh["X-Location"] = thisBH->xLoc;
1470 8 : d_bh["Y-Location"] = thisBH->yLoc;
1471 : }
1472 2 : }
1473 :
1474 : //******************************************************************************
1475 :
1476 2 : void GLHEVert::readCacheFileAndCompareWithThisGLHECache(EnergyPlusData &state)
1477 : {
1478 :
1479 2 : if (!(state.files.outputControl.glhe && FileSystem::fileExists(state.dataStrGlobals->outputGLHEFilePath))) {
1480 : // if the file doesn't exist, there are no data to read
1481 2 : return;
1482 : }
1483 : // file exists -- read data and load if possible
1484 :
1485 0 : auto const cached_json = FileSystem::readJSON(state.dataStrGlobals->outputGLHEFilePath);
1486 :
1487 0 : for (auto const &existing_data : cached_json) {
1488 0 : if (myCacheData["Phys Data"] == existing_data["Phys Data"]) {
1489 0 : myCacheData["Response Factors"] = existing_data["Response Factors"];
1490 0 : gFunctionsExist = true;
1491 0 : break;
1492 : }
1493 : }
1494 :
1495 0 : if (gFunctionsExist) {
1496 : // Populate the time array
1497 0 : this->myRespFactors->time = Array1D<Real64>(myCacheData["Response Factors"]["time"].get<std::vector<Real64>>());
1498 :
1499 : // Populate the lntts array
1500 0 : this->myRespFactors->LNTTS = Array1D<Real64>(myCacheData["Response Factors"]["LNTTS"].get<std::vector<Real64>>());
1501 :
1502 : // Populate the g-function array
1503 0 : this->myRespFactors->GFNC = Array1D<Real64>(myCacheData["Response Factors"]["GFNC"].get<std::vector<Real64>>());
1504 : }
1505 : }
1506 :
1507 : //******************************************************************************
1508 :
1509 2 : void GLHEVert::writeGLHECacheToFile(EnergyPlusData &state) const
1510 : {
1511 :
1512 : // For convenience
1513 : using json = nlohmann::json;
1514 :
1515 4 : json cached_json;
1516 2 : if (FileSystem::fileExists(state.dataStrGlobals->outputGLHEFilePath)) {
1517 : // file exists -- add data
1518 : // open file
1519 0 : cached_json = FileSystem::readJSON(state.dataStrGlobals->outputGLHEFilePath);
1520 :
1521 : // add current data
1522 0 : cached_json.emplace(fmt::format("GHLE {}", cached_json.size() + 1), myCacheData);
1523 : } else {
1524 : // file doesn't exist -- add data
1525 : // add current data
1526 2 : cached_json.emplace("GHLE 1", myCacheData);
1527 : }
1528 2 : FileSystem::writeFile<FileSystem::FileTypes::GLHE>(state.dataStrGlobals->outputGLHEFilePath, cached_json, 2);
1529 2 : }
1530 :
1531 : //******************************************************************************
1532 :
1533 1 : void GLHESlinky::calcGFunctions(EnergyPlusData &state)
1534 : {
1535 : // SUBROUTINE INFORMATION:
1536 : // AUTHOR: Matt Mitchell
1537 : // DATE WRITTEN: February, 2015
1538 : // MODIFIED na
1539 : // RE-ENGINEERED na
1540 :
1541 : // PURPOSE OF THIS SUBROUTINE:
1542 : // calculates g-functions for the slinky ground heat exchanger model
1543 :
1544 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1545 1 : Real64 tLg_max(0.0);
1546 1 : constexpr Real64 tLg_min(-2);
1547 1 : constexpr Real64 tLg_grid(0.25);
1548 1 : constexpr Real64 ts(3600);
1549 : Real64 tLg;
1550 1 : constexpr Real64 convertYearsToSeconds(356 * 24 * 60 * 60);
1551 : Real64 fraction;
1552 2 : Array2D<Real64> valStored({0, this->numTrenches}, {0, this->numCoils}, -1.0);
1553 : int I0;
1554 : int J0;
1555 :
1556 1 : DisplayString(state, "Initializing GroundHeatExchanger:Slinky: " + this->name);
1557 :
1558 1 : this->X0.allocate(this->numCoils);
1559 1 : this->Y0.allocate(this->numTrenches);
1560 :
1561 : // Calculate the number of g-functions required
1562 1 : tLg_max = std::log10(this->maxSimYears * convertYearsToSeconds / ts);
1563 1 : int NPairs = static_cast<int>((tLg_max - tLg_min) / (tLg_grid) + 1);
1564 :
1565 : // Allocate and setup g-function arrays
1566 1 : this->myRespFactors->GFNC.allocate(NPairs);
1567 1 : this->myRespFactors->LNTTS.allocate(NPairs);
1568 1 : this->QnMonthlyAgg.allocate(static_cast<int>(this->maxSimYears * 12));
1569 1 : this->QnHr.allocate(730 + this->AGG + this->SubAGG);
1570 1 : this->QnSubHr.allocate(static_cast<int>((this->SubAGG + 1) * maxTSinHr + 1));
1571 1 : this->LastHourN.allocate(this->SubAGG + 1);
1572 :
1573 29 : for (int i = 1; i <= NPairs; ++i) {
1574 28 : this->myRespFactors->GFNC(i) = 0.0;
1575 28 : this->myRespFactors->LNTTS(i) = 0.0;
1576 : }
1577 :
1578 : // Calculate the number of loops (per trench) and number of trenches to be involved
1579 : // Due to the symmetry of a slinky GHX field, we need only calculate about
1580 : // on quarter of the rings' tube wall temperature perturbation to get the
1581 : // mean wall temperature perturbation of the entire slinky GHX field.
1582 1 : int numLC = std::ceil(this->numCoils / 2.0);
1583 1 : int numRC = std::ceil(this->numTrenches / 2.0);
1584 :
1585 : // Calculate coordinates (X0, Y0, Z0) of a ring's center
1586 201 : for (int coil = 1; coil <= this->numCoils; ++coil) {
1587 200 : this->X0(coil) = coilPitch * (coil - 1);
1588 : }
1589 16 : for (int trench = 1; trench <= this->numTrenches; ++trench) {
1590 15 : this->Y0(trench) = (trench - 1) * this->trenchSpacing;
1591 : }
1592 1 : this->Z0 = this->coilDepth;
1593 :
1594 : // If number of trenches is greater than 1, one quarter of the rings are involved.
1595 : // If number of trenches is 1, one half of the rings are involved.
1596 1 : if (this->numTrenches > 1) {
1597 1 : fraction = 0.25;
1598 : } else {
1599 0 : fraction = 0.5;
1600 : }
1601 :
1602 : // Calculate the corresponding time of each temperature response factor
1603 29 : for (int NT = 1; NT <= NPairs; ++NT) {
1604 28 : tLg = tLg_min + tLg_grid * (NT - 1);
1605 28 : Real64 t = std::pow(10, tLg) * ts;
1606 :
1607 : // Set the average temperature response of the whole field to zero
1608 28 : Real64 gFunc = 0;
1609 :
1610 28 : valStored = -1.0;
1611 :
1612 252 : for (int m1 = 1; m1 <= numRC; ++m1) {
1613 22624 : for (int n1 = 1; n1 <= numLC; ++n1) {
1614 358400 : for (int m = 1; m <= this->numTrenches; ++m) {
1615 67536000 : for (int n = 1; n <= this->numCoils; ++n) {
1616 :
1617 : // Zero out val after each iteration
1618 67200000 : Real64 doubleIntegralVal = 0.0;
1619 67200000 : Real64 midFieldVal = 0.0;
1620 :
1621 : // Calculate the distance between ring centers
1622 67200000 : Real64 disRing = distToCenter(m, n, m1, n1);
1623 :
1624 : // Save mm1 and nn1
1625 67200000 : int mm1 = std::abs(m - m1);
1626 67200000 : int nn1 = std::abs(n - n1);
1627 :
1628 : // If we're calculating a ring's temperature response to itself as a ring source,
1629 : // then we need some extra effort in calculating the double integral
1630 67200000 : if (m1 == m && n1 == n) {
1631 22400 : I0 = 33;
1632 22400 : J0 = 1089;
1633 : } else {
1634 67177600 : I0 = 33;
1635 67177600 : J0 = 561;
1636 : }
1637 :
1638 : Real64 gFuncin;
1639 :
1640 : // if the ring(n1, m1) is the near-field ring of the ring(n,m)
1641 67200000 : if (disRing <= 2.5 + this->coilDiameter) {
1642 : // if no calculated value has been stored
1643 1923628 : if (valStored(mm1, nn1) < 0) {
1644 924 : doubleIntegralVal = doubleIntegral(m, n, m1, n1, t, I0, J0);
1645 924 : valStored(mm1, nn1) = doubleIntegralVal;
1646 : // else: if a stored value is found for the combination of (m, n, m1, n1)
1647 : } else {
1648 1922704 : doubleIntegralVal = valStored(mm1, nn1);
1649 : }
1650 :
1651 : // due to symmetry, the temperature response of ring(n1, m1) should be 0.25, 0.5, or 1 times its calculated value
1652 1923628 : if (!isEven(this->numTrenches) && !isEven(this->numCoils) && m1 == numRC && n1 == numLC && this->numTrenches > 1.5) {
1653 0 : gFuncin = 0.25 * doubleIntegralVal;
1654 1923628 : } else if (!isEven(this->numTrenches) && m1 == numRC && this->numTrenches > 1.5) {
1655 250236 : gFuncin = 0.5 * doubleIntegralVal;
1656 1673392 : } else if (!isEven(this->numCoils) && n1 == numLC) {
1657 0 : gFuncin = 0.5 * doubleIntegralVal;
1658 : } else {
1659 1673392 : gFuncin = doubleIntegralVal;
1660 : }
1661 :
1662 : // if the ring(n1, m1) is in the far-field or the ring(n,m)
1663 65276372 : } else if (disRing > (10 + this->coilDiameter)) {
1664 51011604 : gFuncin = 0;
1665 :
1666 : // else the ring(n1, m1) is in the middle-field of the ring(n,m)
1667 : } else {
1668 : // if no calculated value have been stored
1669 14264768 : if (valStored(mm1, nn1) < 0.0) {
1670 6664 : midFieldVal = midFieldResponseFunction(m, n, m1, n1, t);
1671 6664 : valStored(mm1, nn1) = midFieldVal;
1672 : // if a stored value is found for the combination of (m, n, m1, n1), then
1673 : } else {
1674 14258104 : midFieldVal = valStored(mm1, nn1);
1675 : }
1676 :
1677 : // due to symmetry, the temperature response of ring(n1, m1) should be 0.25, 0.5, or 1 times its calculated value
1678 14264768 : if (!isEven(this->numTrenches) && !isEven(this->numCoils) && m1 == numRC && n1 == numLC && this->numTrenches > 1.5) {
1679 0 : gFuncin = 0.25 * midFieldVal;
1680 14264768 : } else if (!isEven(this->numTrenches) && m1 == numRC && this->numTrenches > 1.5) {
1681 2124864 : gFuncin = 0.5 * midFieldVal;
1682 12139904 : } else if (!isEven(this->numCoils) && n1 == numLC) {
1683 0 : gFuncin = 0.5 * midFieldVal;
1684 : } else {
1685 12139904 : gFuncin = midFieldVal;
1686 : }
1687 : }
1688 :
1689 67200000 : gFunc += gFuncin;
1690 :
1691 : } // n
1692 : } // m
1693 : } // n1
1694 : } // m1
1695 :
1696 28 : this->myRespFactors->GFNC(NT) =
1697 28 : (gFunc * (this->coilDiameter / 2.0)) / (4 * DataGlobalConstants::Pi * fraction * this->numTrenches * this->numCoils);
1698 28 : this->myRespFactors->LNTTS(NT) = tLg;
1699 :
1700 : } // NT time
1701 1 : }
1702 :
1703 : //******************************************************************************
1704 :
1705 1 : void GLHESlinky::makeThisGLHECacheStruct()
1706 : {
1707 1 : }
1708 :
1709 : //******************************************************************************
1710 :
1711 1 : void GLHESlinky::readCacheFileAndCompareWithThisGLHECache([[maybe_unused]] EnergyPlusData &state)
1712 : {
1713 1 : }
1714 :
1715 : //******************************************************************************
1716 :
1717 : Real64
1718 17593884 : GLHESlinky::nearFieldResponseFunction(int const m, int const n, int const m1, int const n1, Real64 const eta, Real64 const theta, Real64 const t)
1719 : {
1720 : // SUBROUTINE INFORMATION:
1721 : // AUTHOR: Matt Mitchell
1722 : // DATE WRITTEN: February, 2015
1723 : // MODIFIED na
1724 : // RE-ENGINEERED na
1725 :
1726 : // PURPOSE OF THIS SUBROUTINE:
1727 : // Calculates the temperature response of from one near-field point to another
1728 :
1729 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1730 17593884 : Real64 distance1 = distance(m, n, m1, n1, eta, theta);
1731 17593884 : Real64 sqrtAlphaT = std::sqrt(this->soil.diffusivity * t);
1732 :
1733 17593884 : if (!verticalConfig) {
1734 :
1735 0 : Real64 sqrtDistDepth = std::sqrt(pow_2(distance1) + 4 * pow_2(this->coilDepth));
1736 0 : Real64 errFunc1 = std::erfc(0.5 * distance1 / sqrtAlphaT);
1737 0 : Real64 errFunc2 = std::erfc(0.5 * sqrtDistDepth / sqrtAlphaT);
1738 :
1739 0 : return errFunc1 / distance1 - errFunc2 / sqrtDistDepth;
1740 :
1741 : } else {
1742 :
1743 17593884 : Real64 distance2 = distanceToFictRing(m, n, m1, n1, eta, theta);
1744 17593884 : Real64 errFunc1 = std::erfc(0.5 * distance1 / sqrtAlphaT);
1745 17593884 : Real64 errFunc2 = std::erfc(0.5 * distance2 / sqrtAlphaT);
1746 :
1747 17593884 : return errFunc1 / distance1 - errFunc2 / distance2;
1748 : }
1749 : }
1750 :
1751 : //******************************************************************************
1752 :
1753 6664 : Real64 GLHESlinky::midFieldResponseFunction(int const m, int const n, int const m1, int const n1, Real64 const t)
1754 : {
1755 : // SUBROUTINE INFORMATION:
1756 : // AUTHOR: Matt Mitchell
1757 : // DATE WRITTEN: February, 2015
1758 : // MODIFIED na
1759 : // RE-ENGINEERED na
1760 :
1761 : // PURPOSE OF THIS SUBROUTINE:
1762 : // Calculates the temperature response of from one mid-field point to another
1763 :
1764 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1765 6664 : Real64 sqrtAlphaT = std::sqrt(this->soil.diffusivity * t);
1766 :
1767 6664 : Real64 distance = distToCenter(m, n, m1, n1);
1768 6664 : Real64 sqrtDistDepth = std::sqrt(pow_2(distance) + 4 * pow_2(this->coilDepth));
1769 :
1770 6664 : Real64 errFunc1 = std::erfc(0.5 * distance / sqrtAlphaT);
1771 6664 : Real64 errFunc2 = std::erfc(0.5 * sqrtDistDepth / sqrtAlphaT);
1772 :
1773 6664 : return 4 * pow_2(DataGlobalConstants::Pi) * (errFunc1 / distance - errFunc2 / sqrtDistDepth);
1774 : }
1775 :
1776 : //******************************************************************************
1777 :
1778 17593884 : Real64 GLHESlinky::distance(int const m, int const n, int const m1, int const n1, Real64 const eta, Real64 const theta)
1779 : {
1780 : // SUBROUTINE INFORMATION:
1781 : // AUTHOR: Matt Mitchell
1782 : // DATE WRITTEN: February, 2015
1783 : // MODIFIED na
1784 : // RE-ENGINEERED na
1785 :
1786 : // PURPOSE OF THIS SUBROUTINE:
1787 : // Calculates the distance between any two points on any two loops
1788 :
1789 17593884 : Real64 const cos_theta = std::cos(theta);
1790 17593884 : Real64 const sin_theta = std::sin(theta);
1791 17593884 : Real64 const cos_eta = std::cos(eta);
1792 17593884 : Real64 const sin_eta = std::sin(eta);
1793 :
1794 17593884 : Real64 x = this->X0(n) + cos_theta * (this->coilDiameter / 2.0);
1795 17593884 : Real64 y = this->Y0(m) + sin_theta * (this->coilDiameter / 2.0);
1796 :
1797 17593884 : Real64 xIn = this->X0(n1) + cos_eta * (this->coilDiameter / 2.0 - this->pipe.outRadius);
1798 17593884 : Real64 yIn = this->Y0(m1) + sin_eta * (this->coilDiameter / 2.0 - this->pipe.outRadius);
1799 :
1800 17593884 : Real64 xOut = this->X0(n1) + cos_eta * (this->coilDiameter / 2.0 + this->pipe.outRadius);
1801 17593884 : Real64 yOut = this->Y0(m1) + sin_eta * (this->coilDiameter / 2.0 + this->pipe.outRadius);
1802 :
1803 17593884 : if (!verticalConfig) {
1804 :
1805 0 : return 0.5 * std::sqrt(pow_2(x - xIn) + pow_2(y - yIn)) + 0.5 * std::sqrt(pow_2(x - xOut) + pow_2(y - yOut));
1806 :
1807 : } else {
1808 :
1809 17593884 : Real64 z = this->Z0 + sin_theta * (this->coilDiameter / 2.0);
1810 :
1811 17593884 : Real64 zIn = this->Z0 + sin_eta * (this->coilDiameter / 2.0 - this->pipe.outRadius);
1812 17593884 : Real64 zOut = this->Z0 + sin_eta * (this->coilDiameter / 2.0 + this->pipe.outRadius);
1813 :
1814 17593884 : return 0.5 * std::sqrt(pow_2(x - xIn) + pow_2(this->Y0(m1) - this->Y0(m)) + pow_2(z - zIn)) +
1815 17593884 : 0.5 * std::sqrt(pow_2(x - xOut) + pow_2(this->Y0(m1) - this->Y0(m)) + pow_2(z - zOut));
1816 : }
1817 : }
1818 :
1819 : //******************************************************************************
1820 :
1821 17593884 : Real64 GLHESlinky::distanceToFictRing(int const m, int const n, int const m1, int const n1, Real64 const eta, Real64 const theta)
1822 : {
1823 : // SUBROUTINE INFORMATION:
1824 : // AUTHOR: Matt Mitchell
1825 : // DATE WRITTEN: February, 2015
1826 : // MODIFIED na
1827 : // RE-ENGINEERED na
1828 :
1829 : // PURPOSE OF THIS SUBROUTINE:
1830 : // Calculates the distance between any two points between real and fictitious rings
1831 :
1832 17593884 : Real64 const sin_theta = std::sin(theta);
1833 17593884 : Real64 const cos_theta = std::cos(theta);
1834 17593884 : Real64 const sin_eta = std::sin(eta);
1835 17593884 : Real64 const cos_eta = std::cos(eta);
1836 :
1837 17593884 : Real64 x = this->X0(n) + cos_theta * (this->coilDiameter / 2.0);
1838 : // Real64 y = Y0( m ) + sin_theta * (coilDiameter / 2.0);
1839 17593884 : Real64 z = this->Z0 + sin_theta * (this->coilDiameter / 2.0) + 2 * this->coilDepth;
1840 :
1841 17593884 : Real64 xIn = this->X0(n1) + cos_eta * (this->coilDiameter / 2.0 - this->pipe.outRadius);
1842 : // Real64 yIn = Y0( m1 ) + sin_eta * (coilDiameter / 2.0 - pipe.outRadius);
1843 17593884 : Real64 zIn = this->Z0 + sin_eta * (this->coilDiameter / 2.0 - this->pipe.outRadius);
1844 :
1845 17593884 : Real64 xOut = this->X0(n1) + cos_eta * (this->coilDiameter / 2.0 + this->pipe.outRadius);
1846 : // Real64 yOut = Y0( m1 ) + sin_eta * (coilDiameter / 2.0 + outRadius);
1847 17593884 : Real64 zOut = this->Z0 + sin_eta * (this->coilDiameter / 2.0 + this->pipe.outRadius);
1848 :
1849 17593884 : return 0.5 * std::sqrt(pow_2(x - xIn) + pow_2(this->Y0(m1) - this->Y0(m)) + pow_2(z - zIn)) +
1850 17593884 : 0.5 * std::sqrt(pow_2(x - xOut) + pow_2(this->Y0(m1) - this->Y0(m)) + pow_2(z - zOut));
1851 : }
1852 :
1853 : //******************************************************************************
1854 :
1855 67206664 : Real64 GLHESlinky::distToCenter(int const m, int const n, int const m1, int const n1)
1856 : {
1857 : // SUBROUTINE INFORMATION:
1858 : // AUTHOR: Matt Mitchell
1859 : // DATE WRITTEN: February, 2015
1860 : // MODIFIED na
1861 : // RE-ENGINEERED na
1862 :
1863 : // PURPOSE OF THIS SUBROUTINE:
1864 : // Calculates the center-to-center distance between rings
1865 :
1866 67206664 : return std::sqrt(pow_2(this->X0(n) - this->X0(n1)) + pow_2(this->Y0(m) - this->Y0(m1)));
1867 : }
1868 :
1869 : //******************************************************************************
1870 :
1871 84966236 : inline bool GLHEBase::isEven(int const val)
1872 : {
1873 : // SUBROUTINE INFORMATION:
1874 : // AUTHOR: Matt Mitchell
1875 : // DATE WRITTEN: February, 2015
1876 : // MODIFIED na
1877 : // RE-ENGINEERED na
1878 :
1879 : // PURPOSE OF THIS SUBROUTINE:
1880 : // Determines if an integer is even
1881 :
1882 84966236 : if (val % 2 == 0) {
1883 41306700 : return true;
1884 : } else {
1885 43659536 : return false;
1886 : }
1887 : }
1888 :
1889 : //******************************************************************************
1890 :
1891 30492 : Real64 GLHESlinky::integral(int const m, int const n, int const m1, int const n1, Real64 const t, Real64 const eta, Real64 const J0)
1892 : {
1893 : // SUBROUTINE INFORMATION:
1894 : // AUTHOR: Matt Mitchell
1895 : // DATE WRITTEN: February, 2015
1896 : // MODIFIED na
1897 : // RE-ENGINEERED na
1898 :
1899 : // PURPOSE OF THIS SUBROUTINE:
1900 : // Integrates the temperature response at one point based on
1901 : // input from other points
1902 :
1903 : // METHODOLOGY EMPLOYED:
1904 : // Simpson's 1/3 rule of integration
1905 :
1906 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1907 30492 : Real64 sumIntF(0.0);
1908 30492 : Real64 theta(0.0);
1909 30492 : constexpr Real64 theta1(0.0);
1910 30492 : constexpr Real64 theta2(2 * DataGlobalConstants::Pi);
1911 60984 : Array1D<Real64> f(J0, 0.0);
1912 :
1913 30492 : Real64 h = (theta2 - theta1) / (J0 - 1);
1914 :
1915 : // Calculate the function at various equally spaced x values
1916 17624376 : for (int j = 1; j <= J0; ++j) {
1917 :
1918 17593884 : theta = theta1 + (j - 1) * h;
1919 :
1920 17593884 : f(j) = nearFieldResponseFunction(m, n, m1, n1, eta, theta, t);
1921 :
1922 17593884 : if (j == 1 || j == J0) {
1923 60984 : f(j) = f(j);
1924 17532900 : } else if (isEven(j)) {
1925 8781696 : f(j) = 4 * f(j);
1926 : } else {
1927 8751204 : f(j) = 2 * f(j);
1928 : }
1929 :
1930 17593884 : sumIntF += f(j);
1931 : }
1932 :
1933 60984 : return (h / 3) * sumIntF;
1934 : }
1935 :
1936 : //******************************************************************************
1937 :
1938 924 : Real64 GLHESlinky::doubleIntegral(int const m, int const n, int const m1, int const n1, Real64 const t, int const I0, int const J0)
1939 : {
1940 : // SUBROUTINE INFORMATION:
1941 : // AUTHOR: Matt Mitchell
1942 : // DATE WRITTEN: February, 2015
1943 : // MODIFIED na
1944 : // RE-ENGINEERED na
1945 :
1946 : // PURPOSE OF THIS SUBROUTINE:
1947 : // Integrates the temperature response at one point based on
1948 : // input from other points
1949 :
1950 : // METHODOLOGY EMPLOYED:
1951 : // Simpson's 1/3 rule of integration
1952 :
1953 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
1954 924 : constexpr Real64 eta1(0.0);
1955 924 : constexpr Real64 eta2(2 * DataGlobalConstants::Pi);
1956 :
1957 924 : Real64 sumIntF(0.0);
1958 1848 : Array1D<Real64> g(I0, 0.0);
1959 :
1960 924 : Real64 h = (eta2 - eta1) / (I0 - 1);
1961 :
1962 : // Calculates the value of the function at various equally spaced values
1963 31416 : for (int i = 1; i <= I0; ++i) {
1964 :
1965 30492 : Real64 eta = eta1 + (i - 1) * h;
1966 30492 : g(i) = integral(m, n, m1, n1, t, eta, J0);
1967 :
1968 30492 : if (i == 1 || i == I0) {
1969 1848 : g(i) = g(i);
1970 28644 : } else if (isEven(i)) {
1971 14784 : g(i) = 4 * g(i);
1972 : } else {
1973 13860 : g(i) = 2 * g(i);
1974 : }
1975 :
1976 30492 : sumIntF += g(i);
1977 : }
1978 :
1979 1848 : return (h / 3) * sumIntF;
1980 : }
1981 :
1982 : //******************************************************************************
1983 :
1984 382208 : void GLHEVert::getAnnualTimeConstant()
1985 : {
1986 : // SUBROUTINE INFORMATION:
1987 : // AUTHOR: Matt Mitchell
1988 : // DATE WRITTEN: February, 2015
1989 : // MODIFIED na
1990 : // RE-ENGINEERED na
1991 :
1992 : // PURPOSE OF THIS SUBROUTINE:
1993 : // calculate annual time constant for ground conduction
1994 :
1995 382208 : constexpr Real64 hrInYear(8760);
1996 :
1997 382208 : this->timeSS = (pow_2(this->bhLength) / (9.0 * this->soil.diffusivity)) / DataGlobalConstants::SecInHour / hrInYear;
1998 382208 : this->timeSSFactor = this->timeSS * 8760.0;
1999 382208 : }
2000 :
2001 : //******************************************************************************
2002 :
2003 14520 : void GLHESlinky::getAnnualTimeConstant()
2004 : {
2005 : // SUBROUTINE INFORMATION:
2006 : // AUTHOR: Matt Mitchell
2007 : // DATE WRITTEN: February, 2015
2008 : // MODIFIED na
2009 : // RE-ENGINEERED na
2010 :
2011 : // PURPOSE OF THIS SUBROUTINE:
2012 : // calculate annual time constant for ground conduction
2013 :
2014 14520 : this->timeSSFactor = 1.0;
2015 14520 : }
2016 :
2017 : //******************************************************************************
2018 :
2019 396728 : void GLHEBase::calcGroundHeatExchanger(EnergyPlusData &state)
2020 : {
2021 : // SUBROUTINE INFORMATION:
2022 : // AUTHOR: Dan Fisher
2023 : // DATE WRITTEN: August, 2000
2024 : // MODIFIED Arun Murugappan
2025 : // RE-ENGINEERED na
2026 :
2027 : // PURPOSE OF THIS SUBROUTINE:
2028 : // This is the main routine to simulate the operation of vertical
2029 : // closed-loop ground heat exchangers (GLHE).
2030 :
2031 : // METHODOLOGY EMPLOYED:
2032 : // The borehole and fluid temperatures are calculated from the response to
2033 : // the current heat transfer rate and the response to the history of past
2034 : // applied heat pulses. The response to each pulse is calculated from a non-
2035 : // dimensionalized response function, or G-function, that is specific to the
2036 : // given borehole field arrangement, depth and spacing. The data defining
2037 : // this function is read from input.
2038 : // The heat pulse histories need to be recorded over an extended period (months).
2039 : // To aid computational efficiency past pulses are continuously aggregated into
2040 : // equivalent heat pulses of longer duration, as each pulse becomes less recent.
2041 :
2042 : // REFERENCES:
2043 : // Eskilson, P. 'Thermal Analysis of Heat Extraction Boreholes' Ph.D. Thesis:
2044 : // Dept. of Mathematical Physics, University of Lund, Sweden, June 1987.
2045 : // Yavuzturk, C., J.D. Spitler. 1999. 'A Short Time Step Response Factor Model
2046 : // for Vertical Ground Loop Heat Exchangers.' ASHRAE Transactions. 105(2): 475-485.
2047 :
2048 : // Using/Aliasing
2049 : using FluidProperties::GetDensityGlycol;
2050 : using FluidProperties::GetSpecificHeatGlycol;
2051 :
2052 : // SUBROUTINE ARGUMENT DEFINITIONS
2053 396728 : constexpr const char *RoutineName("CalcGroundHeatExchanger");
2054 :
2055 : // LOCAL PARAMETERS
2056 : Real64 fluidAveTemp;
2057 : Real64 tmpQnSubHourly; // current Qn sub-hourly value
2058 396728 : Real64 sumTotal(0.0); // sum of all the Qn (load) blocks
2059 :
2060 : // Calculate G-Functions
2061 396728 : if (this->firstTime) {
2062 23 : if (!gFunctionsExist) {
2063 3 : makeThisGLHECacheAndCompareWithFileCache(state);
2064 3 : if (!gFunctionsExist) {
2065 3 : calcGFunctions(state);
2066 3 : gFunctionsExist = true;
2067 : }
2068 : }
2069 23 : this->firstTime = false;
2070 : }
2071 :
2072 396728 : this->inletTemp = state.dataLoopNodes->Node(this->inletNodeNum).Temp;
2073 :
2074 1190184 : Real64 cpFluid = GetSpecificHeatGlycol(state,
2075 396728 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2076 : this->inletTemp,
2077 396728 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2078 793456 : RoutineName);
2079 :
2080 396728 : Real64 kGroundFactor = 2.0 * DataGlobalConstants::Pi * this->soil.k;
2081 :
2082 : // Get time constants
2083 396728 : getAnnualTimeConstant();
2084 :
2085 396728 : if (triggerDesignDayReset && state.dataGlobal->WarmupFlag) updateCurSimTime = true;
2086 396728 : if (state.dataGlobal->DayOfSim == 1 && updateCurSimTime) {
2087 90 : state.dataGroundHeatExchanger->currentSimTime = 0.0;
2088 90 : state.dataGroundHeatExchanger->prevTimeSteps = 0.0;
2089 90 : this->QnHr = 0.0;
2090 90 : this->QnMonthlyAgg = 0.0;
2091 90 : this->QnSubHr = 0.0;
2092 90 : this->LastHourN = 1;
2093 90 : state.dataGroundHeatExchanger->N = 1;
2094 90 : updateCurSimTime = false;
2095 90 : triggerDesignDayReset = false;
2096 : }
2097 :
2098 1190184 : state.dataGroundHeatExchanger->currentSimTime = (state.dataGlobal->DayOfSim - 1) * 24 + state.dataGlobal->HourOfDay - 1 +
2099 793456 : (state.dataGlobal->TimeStep - 1) * state.dataGlobal->TimeStepZone +
2100 396728 : state.dataHVACGlobal->SysTimeElapsed; //+ TimeStepsys
2101 396728 : state.dataGroundHeatExchanger->locHourOfDay =
2102 396728 : static_cast<int>(mod(state.dataGroundHeatExchanger->currentSimTime, DataGlobalConstants::HoursInDay) + 1);
2103 396728 : state.dataGroundHeatExchanger->locDayOfSim = static_cast<int>(state.dataGroundHeatExchanger->currentSimTime / 24 + 1);
2104 :
2105 396728 : if (state.dataGlobal->DayOfSim > 1) {
2106 288224 : updateCurSimTime = true;
2107 : }
2108 :
2109 396728 : if (!state.dataGlobal->WarmupFlag) {
2110 51096 : triggerDesignDayReset = true;
2111 : }
2112 :
2113 396728 : if (state.dataGroundHeatExchanger->currentSimTime <= 0.0) {
2114 912 : state.dataGroundHeatExchanger->prevTimeSteps = 0.0; // This resets history when rounding 24:00 hours during warmup avoids hard crash later
2115 912 : calcAggregateLoad(state); // Just allocates and initializes prevHour array
2116 912 : return;
2117 : }
2118 :
2119 : // Store currentSimTime in prevTimeSteps only if a time step occurs
2120 395816 : if (state.dataGroundHeatExchanger->prevTimeSteps(1) != state.dataGroundHeatExchanger->currentSimTime) {
2121 90600 : state.dataGroundHeatExchanger->prevTimeSteps =
2122 135900 : eoshift(state.dataGroundHeatExchanger->prevTimeSteps, -1, state.dataGroundHeatExchanger->currentSimTime);
2123 45300 : ++state.dataGroundHeatExchanger->N;
2124 : }
2125 :
2126 395816 : if (state.dataGroundHeatExchanger->N != PrevN) {
2127 45300 : PrevN = state.dataGroundHeatExchanger->N;
2128 45300 : this->QnSubHr = eoshift(this->QnSubHr, -1, this->lastQnSubHr);
2129 : }
2130 :
2131 395816 : calcAggregateLoad(state);
2132 :
2133 : // Update the heat exchanger resistance each time
2134 395816 : this->HXResistance = calcHXResistance(state);
2135 :
2136 395816 : if (state.dataGroundHeatExchanger->N == 1) {
2137 0 : if (this->massFlowRate <= 0.0) {
2138 0 : tmpQnSubHourly = 0.0;
2139 0 : fluidAveTemp = this->tempGround;
2140 0 : this->ToutNew = this->inletTemp;
2141 : } else {
2142 0 : Real64 gFuncVal = getGFunc(state.dataGroundHeatExchanger->currentSimTime / (this->timeSSFactor));
2143 :
2144 0 : Real64 C_1 = (this->totalTubeLength) / (2.0 * this->massFlowRate * cpFluid);
2145 0 : tmpQnSubHourly = (this->tempGround - this->inletTemp) / (gFuncVal / (kGroundFactor) + this->HXResistance + C_1);
2146 0 : fluidAveTemp = this->tempGround - tmpQnSubHourly * this->HXResistance;
2147 0 : this->ToutNew = this->tempGround - tmpQnSubHourly * (gFuncVal / (kGroundFactor) + this->HXResistance - C_1);
2148 : }
2149 : } else {
2150 : // no monthly super position
2151 395816 : if (state.dataGroundHeatExchanger->currentSimTime < (hrsPerMonth + this->AGG + this->SubAGG)) {
2152 :
2153 : // Calculate the Sub Hourly Superposition
2154 :
2155 : // same as above for sub-hourly( with no aggregation]
2156 395816 : Real64 sumQnSubHourly = 0.0;
2157 : int IndexN;
2158 395816 : if (int(state.dataGroundHeatExchanger->currentSimTime) < this->SubAGG) {
2159 66376 : IndexN = int(state.dataGroundHeatExchanger->currentSimTime) + 1;
2160 : } else {
2161 329440 : IndexN = this->SubAGG + 1;
2162 : }
2163 :
2164 395816 : int subHourlyLimit = state.dataGroundHeatExchanger->N - this->LastHourN(IndexN); // Check this when running simulation
2165 36024952 : for (int I = 1; I <= subHourlyLimit; ++I) {
2166 36024600 : if (I == subHourlyLimit) {
2167 395464 : if (int(state.dataGroundHeatExchanger->currentSimTime) >= this->SubAGG) {
2168 : Real64 gFuncVal =
2169 658880 : getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(I + 1)) /
2170 988320 : (this->timeSSFactor));
2171 329440 : Real64 RQSubHr = gFuncVal / (kGroundFactor);
2172 329440 : sumQnSubHourly += (this->QnSubHr(I) - this->QnHr(IndexN)) * RQSubHr;
2173 : } else {
2174 : Real64 gFuncVal =
2175 132048 : getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(I + 1)) /
2176 198072 : (this->timeSSFactor));
2177 66024 : Real64 RQSubHr = gFuncVal / (kGroundFactor);
2178 66024 : sumQnSubHourly += this->QnSubHr(I) * RQSubHr;
2179 : }
2180 395464 : break;
2181 : }
2182 : // prevTimeSteps(I+1) This is "I+1" because prevTimeSteps(1) = CurrentTimestep
2183 71258272 : Real64 gFuncVal = getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(I + 1)) /
2184 106887408 : (this->timeSSFactor));
2185 35629136 : Real64 RQSubHr = gFuncVal / (kGroundFactor);
2186 35629136 : sumQnSubHourly += (this->QnSubHr(I) - this->QnSubHr(I + 1)) * RQSubHr;
2187 : }
2188 :
2189 : // Calculate the Hourly Superposition
2190 : // same as above for hourly
2191 395816 : Real64 sumQnHourly = 0.0;
2192 :
2193 395816 : int hourlyLimit = int(state.dataGroundHeatExchanger->currentSimTime);
2194 28074888 : for (int I = this->SubAGG + 1; I <= hourlyLimit; ++I) {
2195 28005216 : if (I == hourlyLimit) {
2196 326144 : Real64 gFuncVal = getGFunc(state.dataGroundHeatExchanger->currentSimTime / (this->timeSSFactor));
2197 326144 : Real64 RQHour = gFuncVal / (kGroundFactor);
2198 326144 : sumQnHourly += this->QnHr(I) * RQHour;
2199 326144 : break;
2200 : }
2201 55358144 : Real64 gFuncVal = getGFunc((state.dataGroundHeatExchanger->currentSimTime - int(state.dataGroundHeatExchanger->currentSimTime) + I) /
2202 83037216 : (this->timeSSFactor));
2203 27679072 : Real64 RQHour = gFuncVal / (kGroundFactor);
2204 27679072 : sumQnHourly += (this->QnHr(I) - this->QnHr(I + 1)) * RQHour;
2205 : }
2206 :
2207 : // Find the total Sum of the Temperature difference due to all load blocks
2208 395816 : sumTotal = sumQnSubHourly + sumQnHourly;
2209 :
2210 : // Calculate the sub-hourly temperature due the Last Time steps Load
2211 : Real64 gFuncVal =
2212 395816 : getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(2)) / (this->timeSSFactor));
2213 395816 : Real64 RQSubHr = gFuncVal / kGroundFactor;
2214 :
2215 395816 : if (this->massFlowRate <= 0.0) {
2216 93130 : tmpQnSubHourly = 0.0;
2217 93130 : fluidAveTemp = this->tempGround - sumTotal; // Q(N)*RB = 0
2218 93130 : this->ToutNew = this->inletTemp;
2219 : } else {
2220 : // Dr.Sitler's Explicit set of equations to calculate the New Outlet Temperature of the U-Tube
2221 302686 : Real64 C0 = RQSubHr;
2222 302686 : Real64 C1 = this->tempGround - (sumTotal - this->QnSubHr(1) * RQSubHr);
2223 302686 : Real64 C2 = this->totalTubeLength / (2.0 * this->massFlowRate * cpFluid);
2224 302686 : Real64 C3 = this->massFlowRate * cpFluid / (this->totalTubeLength);
2225 302686 : tmpQnSubHourly = (C1 - this->inletTemp) / (this->HXResistance + C0 - C2 + (1 / C3));
2226 302686 : fluidAveTemp = C1 - (C0 + this->HXResistance) * tmpQnSubHourly;
2227 302686 : this->ToutNew = C1 + (C2 - C0 - this->HXResistance) * tmpQnSubHourly;
2228 : }
2229 :
2230 : } else { // Monthly Aggregation and super position
2231 :
2232 : // the number of months of simulation elapsed
2233 0 : int numOfMonths = static_cast<int>((state.dataGroundHeatExchanger->currentSimTime + 1) / hrsPerMonth);
2234 :
2235 : // The Month up to which the monthly blocks are superposed
2236 : int currentMonth;
2237 :
2238 0 : if (state.dataGroundHeatExchanger->currentSimTime < ((numOfMonths)*hrsPerMonth) + this->AGG + this->SubAGG) {
2239 0 : currentMonth = numOfMonths - 1;
2240 : } else {
2241 0 : currentMonth = numOfMonths;
2242 : }
2243 :
2244 : // Monthly superposition
2245 : // tmp variable which holds the sum of the Temperature difference due to Aggregated heat extraction/rejection step
2246 0 : Real64 sumQnMonthly = 0.0;
2247 :
2248 0 : for (int I = 1; I <= currentMonth; ++I) {
2249 0 : if (I == 1) {
2250 0 : Real64 gFuncVal = getGFunc(state.dataGroundHeatExchanger->currentSimTime / (this->timeSSFactor));
2251 0 : Real64 RQMonth = gFuncVal / (kGroundFactor);
2252 0 : sumQnMonthly += this->QnMonthlyAgg(I) * RQMonth;
2253 0 : continue;
2254 : }
2255 0 : Real64 gFuncVal = getGFunc((state.dataGroundHeatExchanger->currentSimTime - (I - 1) * hrsPerMonth) / (this->timeSSFactor));
2256 0 : Real64 RQMonth = gFuncVal / (kGroundFactor);
2257 0 : sumQnMonthly += (this->QnMonthlyAgg(I) - this->QnMonthlyAgg(I - 1)) * RQMonth;
2258 : }
2259 :
2260 : // Hourly Superposition
2261 0 : Real64 sumQnHourly = 0.0;
2262 0 : int hourlyLimit = int(state.dataGroundHeatExchanger->currentSimTime - currentMonth * hrsPerMonth);
2263 0 : for (int I = 1 + this->SubAGG; I <= hourlyLimit; ++I) {
2264 0 : if (I == hourlyLimit) {
2265 : Real64 gFuncVal =
2266 0 : getGFunc((state.dataGroundHeatExchanger->currentSimTime - int(state.dataGroundHeatExchanger->currentSimTime) + I) /
2267 0 : (this->timeSSFactor));
2268 0 : Real64 RQHour = gFuncVal / (kGroundFactor);
2269 0 : sumQnHourly += (this->QnHr(I) - this->QnMonthlyAgg(currentMonth)) * RQHour;
2270 0 : break;
2271 : }
2272 0 : Real64 gFuncVal = getGFunc((state.dataGroundHeatExchanger->currentSimTime - int(state.dataGroundHeatExchanger->currentSimTime) + I) /
2273 0 : (this->timeSSFactor));
2274 0 : Real64 RQHour = gFuncVal / (kGroundFactor);
2275 0 : sumQnHourly += (this->QnHr(I) - this->QnHr(I + 1)) * RQHour;
2276 : }
2277 :
2278 : // sub-hourly Superposition
2279 0 : int subHourlyLimit = state.dataGroundHeatExchanger->N - this->LastHourN(this->SubAGG + 1);
2280 0 : Real64 sumQnSubHourly = 0.0;
2281 0 : for (int I = 1; I <= subHourlyLimit; ++I) {
2282 0 : if (I == subHourlyLimit) {
2283 0 : Real64 gFuncVal = getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(I + 1)) /
2284 0 : (this->timeSSFactor));
2285 0 : Real64 RQSubHr = gFuncVal / (kGroundFactor);
2286 0 : sumQnSubHourly += (this->QnSubHr(I) - this->QnHr(this->SubAGG + 1)) * RQSubHr;
2287 0 : break;
2288 : }
2289 0 : Real64 gFuncVal = getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(I + 1)) /
2290 0 : (this->timeSSFactor));
2291 0 : Real64 RQSubHr = gFuncVal / (kGroundFactor);
2292 0 : sumQnSubHourly += (this->QnSubHr(I) - this->QnSubHr(I + 1)) * RQSubHr;
2293 : }
2294 :
2295 0 : sumTotal = sumQnMonthly + sumQnHourly + sumQnSubHourly;
2296 :
2297 : // Calculate the sub-hourly temperature due the Last Time steps Load
2298 : Real64 gFuncVal =
2299 0 : getGFunc((state.dataGroundHeatExchanger->currentSimTime - state.dataGroundHeatExchanger->prevTimeSteps(2)) / (this->timeSSFactor));
2300 0 : Real64 RQSubHr = gFuncVal / (kGroundFactor);
2301 :
2302 0 : if (this->massFlowRate <= 0.0) {
2303 0 : tmpQnSubHourly = 0.0;
2304 0 : fluidAveTemp = this->tempGround - sumTotal; // Q(N)*RB = 0
2305 0 : this->ToutNew = this->inletTemp;
2306 : } else {
2307 : // Explicit set of equations to calculate the New Outlet Temperature of the U-Tube
2308 0 : Real64 C0 = RQSubHr;
2309 0 : Real64 C1 = this->tempGround - (sumTotal - this->QnSubHr(1) * RQSubHr);
2310 0 : Real64 C2 = this->totalTubeLength / (2 * this->massFlowRate * cpFluid);
2311 0 : Real64 C3 = this->massFlowRate * cpFluid / (this->totalTubeLength);
2312 0 : tmpQnSubHourly = (C1 - this->inletTemp) / (this->HXResistance + C0 - C2 + (1 / C3));
2313 0 : fluidAveTemp = C1 - (C0 + this->HXResistance) * tmpQnSubHourly;
2314 0 : this->ToutNew = C1 + (C2 - C0 - this->HXResistance) * tmpQnSubHourly;
2315 : }
2316 : } // end of AGG OR NO AGG
2317 : } // end of N = 1 branch
2318 395816 : this->bhTemp = this->tempGround - sumTotal;
2319 :
2320 : // Load the QnSubHourly Array with a new value at end of every timestep
2321 395816 : this->lastQnSubHr = tmpQnSubHourly;
2322 395816 : this->outletTemp = this->ToutNew;
2323 395816 : this->QGLHE = tmpQnSubHourly * this->totalTubeLength;
2324 395816 : this->aveFluidTemp = fluidAveTemp;
2325 : }
2326 :
2327 : //******************************************************************************
2328 :
2329 396728 : void GLHEBase::updateGHX(EnergyPlusData &state)
2330 : {
2331 : // SUBROUTINE INFORMATION:
2332 : // AUTHOR: Matt Mitchell
2333 : // DATE WRITTEN: February, 2015
2334 : // MODIFIED: na
2335 : // RE-ENGINEERED: na
2336 :
2337 : // PURPOSE OF THIS SUBROUTINE:
2338 : // Updates the outlet node and check for out of bounds temperatures
2339 :
2340 : // Using/Aliasing
2341 : using FluidProperties::GetDensityGlycol;
2342 : using FluidProperties::GetSpecificHeatGlycol;
2343 : using PlantUtilities::SafeCopyPlantNode;
2344 :
2345 : // SUBROUTINE ARGUMENT DEFINITIONS
2346 396728 : constexpr const char *RoutineName("UpdateGroundHeatExchanger");
2347 :
2348 396728 : constexpr Real64 deltaTempLimit(100.0); // temp limit for warnings
2349 :
2350 396728 : SafeCopyPlantNode(state, this->inletNodeNum, this->outletNodeNum);
2351 :
2352 396728 : state.dataLoopNodes->Node(this->outletNodeNum).Temp = this->outletTemp;
2353 396728 : state.dataLoopNodes->Node(this->outletNodeNum).Enthalpy =
2354 1586912 : this->outletTemp * GetSpecificHeatGlycol(state,
2355 396728 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2356 : this->outletTemp,
2357 396728 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2358 396728 : RoutineName);
2359 :
2360 396728 : Real64 GLHEdeltaTemp = std::abs(this->outletTemp - this->inletTemp);
2361 :
2362 396728 : if (GLHEdeltaTemp > deltaTempLimit && this->numErrorCalls < state.dataGroundHeatExchanger->numVerticalGLHEs && !state.dataGlobal->WarmupFlag) {
2363 0 : Real64 fluidDensity = GetDensityGlycol(state,
2364 0 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2365 : this->inletTemp,
2366 0 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2367 0 : RoutineName);
2368 0 : this->designMassFlow = this->designFlow * fluidDensity;
2369 0 : ShowWarningError(state, "Check GLHE design inputs & g-functions for consistency");
2370 0 : ShowContinueError(state, "For GroundHeatExchanger: " + this->name + "GLHE delta Temp > 100C.");
2371 0 : ShowContinueError(state, "This can be encountered in cases where the GLHE mass flow rate is either significantly");
2372 0 : ShowContinueError(state, " lower than the design value, or cases where the mass flow rate rapidly changes.");
2373 0 : ShowContinueError(state, format("GLHE Current Flow Rate={:.3T}; GLHE Design Flow Rate={:.3T}", this->massFlowRate, this->designMassFlow));
2374 0 : ++this->numErrorCalls;
2375 : }
2376 396728 : }
2377 :
2378 : //******************************************************************************
2379 :
2380 396728 : void GLHEBase::calcAggregateLoad(EnergyPlusData &state)
2381 : {
2382 :
2383 : // SUBROUTINE INFORMATION:
2384 : // AUTHOR: Arun Murugappan
2385 : // DATE WRITTEN: August, 2000
2386 : // MODIFIED: na
2387 : // RE-ENGINEERED: na
2388 :
2389 : // PURPOSE OF THIS SUBROUTINE:
2390 : // Manages the heat transfer history.
2391 :
2392 : // METHODOLOGY EMPLOYED:
2393 : // The heat pulse histories need to be recorded over an extended period (months).
2394 : // To aid computational efficiency past pulses are continuously aggregated into
2395 : // equivalent heat pulses of longer duration, as each pulse becomes less recent.
2396 : // Past sub-hourly loads are re-aggregated into equivalent hourly and monthly loads.
2397 :
2398 : // REFERENCES:
2399 : // Eskilson, P. 'Thermal Analysis of Heat Extraction Boreholes' Ph.D. Thesis:
2400 : // Dept. of Mathematical Physics, University of Lund, Sweden, June 1987.
2401 : // Yavuzturk, C., J.D. Spitler. 1999. 'A Short Time Step Response Factor Model
2402 : // for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2): 475-485.
2403 :
2404 396728 : if (state.dataGroundHeatExchanger->currentSimTime <= 0.0) return;
2405 :
2406 : // FOR EVERY HOUR UPDATE THE HOURLY QN this->QnHr(J)
2407 : // THIS IS DONE BY AGGREGATING THE sub-hourly QN FROM THE PREVIOUS HOUR TO UNTIL THE CURRENT HOUR
2408 : // AND STORING IT IN verticalGLHE(GLHENum)%QnHr(J)
2409 :
2410 : // sub-hourly Qn IS NOT AGGREGATED . IT IS THE BASIC LOAD
2411 395816 : if (this->prevHour != state.dataGroundHeatExchanger->locHourOfDay) {
2412 7824 : Real64 SumQnHr = 0.0;
2413 : int J;
2414 52804 : for (J = 1; J <= (state.dataGroundHeatExchanger->N - this->LastHourN(1)); ++J) {
2415 44980 : SumQnHr +=
2416 44980 : this->QnSubHr(J) * std::abs(state.dataGroundHeatExchanger->prevTimeSteps(J) - state.dataGroundHeatExchanger->prevTimeSteps(J + 1));
2417 : }
2418 7824 : if (state.dataGroundHeatExchanger->prevTimeSteps(1) != state.dataGroundHeatExchanger->prevTimeSteps(J)) {
2419 7824 : SumQnHr /= std::abs(state.dataGroundHeatExchanger->prevTimeSteps(1) - state.dataGroundHeatExchanger->prevTimeSteps(J));
2420 : } else {
2421 0 : SumQnHr /= 0.05; // estimated small timestep
2422 : }
2423 7824 : this->QnHr = eoshift(this->QnHr, -1, SumQnHr);
2424 7824 : this->LastHourN = eoshift(this->LastHourN, -1, state.dataGroundHeatExchanger->N);
2425 : }
2426 :
2427 : // CHECK IF A MONTH PASSES...
2428 791632 : if (mod(((state.dataGroundHeatExchanger->locDayOfSim - 1) * DataGlobalConstants::HoursInDay + (state.dataGroundHeatExchanger->locHourOfDay)),
2429 395816 : hrsPerMonth) == 0 &&
2430 0 : this->prevHour != state.dataGroundHeatExchanger->locHourOfDay) {
2431 0 : Real64 MonthNum = static_cast<int>(
2432 0 : (state.dataGroundHeatExchanger->locDayOfSim * DataGlobalConstants::HoursInDay + state.dataGroundHeatExchanger->locHourOfDay) /
2433 0 : hrsPerMonth);
2434 0 : Real64 SumQnMonth = 0.0;
2435 0 : for (int J = 1; J <= int(hrsPerMonth); ++J) {
2436 0 : SumQnMonth += this->QnHr(J);
2437 : }
2438 0 : SumQnMonth /= hrsPerMonth;
2439 0 : this->QnMonthlyAgg(MonthNum) = SumQnMonth;
2440 : }
2441 395816 : if (this->prevHour != state.dataGroundHeatExchanger->locHourOfDay) {
2442 7824 : this->prevHour = state.dataGroundHeatExchanger->locHourOfDay;
2443 : }
2444 : }
2445 :
2446 : //******************************************************************************
2447 :
2448 23 : void GetGroundHeatExchangerInput(EnergyPlusData &state)
2449 : {
2450 : // SUBROUTINE INFORMATION:
2451 : // AUTHOR: Dan Fisher
2452 : // DATE WRITTEN: August, 2000
2453 : // MODIFIED Arun Murugappan
2454 : // RE-ENGINEERED na
2455 :
2456 23 : bool errorsFound = false;
2457 :
2458 : // GET NUMBER OF ALL EQUIPMENT TYPES
2459 23 : state.dataGroundHeatExchanger->numVerticalGLHEs =
2460 46 : state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "GroundHeatExchanger:System");
2461 23 : state.dataGroundHeatExchanger->numSlinkyGLHEs =
2462 46 : state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "GroundHeatExchanger:Slinky");
2463 23 : state.dataGroundHeatExchanger->numVertArray =
2464 46 : state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "GroundHeatExchanger:Vertical:Array");
2465 23 : state.dataGroundHeatExchanger->numVertProps =
2466 46 : state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "GroundHeatExchanger:Vertical:Properties");
2467 23 : state.dataGroundHeatExchanger->numResponseFactors =
2468 46 : state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "GroundHeatExchanger:ResponseFactors");
2469 23 : state.dataGroundHeatExchanger->numSingleBorehole =
2470 46 : state.dataInputProcessing->inputProcessor->getNumObjectsFound(state, "GroundHeatExchanger:Vertical:Single");
2471 :
2472 23 : if (state.dataGroundHeatExchanger->numVerticalGLHEs <= 0 && state.dataGroundHeatExchanger->numSlinkyGLHEs <= 0) {
2473 0 : ShowSevereError(state, "Error processing inputs for GLHE objects");
2474 0 : ShowContinueError(state, "Simulation indicated these objects were found, but input processor doesn't find any");
2475 0 : ShowContinueError(state, "Check inputs for GroundHeatExchanger:System and GroundHeatExchanger:Slinky");
2476 0 : ShowContinueError(state, "Also check plant/branch inputs for references to invalid/deleted objects");
2477 0 : errorsFound = true;
2478 : }
2479 :
2480 23 : if (state.dataGroundHeatExchanger->numVertProps > 0) {
2481 :
2482 44 : std::string currObj = "GroundHeatExchanger:Vertical:Properties";
2483 :
2484 44 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(currObj);
2485 22 : if (instances == state.dataInputProcessing->inputProcessor->epJSON.end()) {
2486 : ShowSevereError(state, // LCOV_EXCL_LINE
2487 : currObj + ": Somehow getNumObjectsFound was > 0 but epJSON.find found 0"); // LCOV_EXCL_LINE
2488 : }
2489 :
2490 22 : auto &instancesValue = instances.value();
2491 44 : for (auto it = instancesValue.begin(); it != instancesValue.end(); ++it) {
2492 22 : auto const &instance = it.value();
2493 22 : auto const &objName = it.key();
2494 44 : auto const &objNameUC = UtilityRoutines::MakeUPPERCase(objName);
2495 22 : state.dataInputProcessing->inputProcessor->markObjectAsUsed(currObj, objName);
2496 44 : std::shared_ptr<GLHEVertProps> thisObj(new GLHEVertProps(state, objNameUC, instance));
2497 22 : state.dataGroundHeatExchanger->vertPropsVector.push_back(thisObj);
2498 : }
2499 : }
2500 :
2501 23 : if (state.dataGroundHeatExchanger->numResponseFactors > 0) {
2502 :
2503 40 : std::string currObj = "GroundHeatExchanger:ResponseFactors";
2504 :
2505 40 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(currObj);
2506 20 : if (instances == state.dataInputProcessing->inputProcessor->epJSON.end()) {
2507 : ShowSevereError(state, // LCOV_EXCL_LINE
2508 : currObj + ": Somehow getNumObjectsFound was > 0 but epJSON.find found 0"); // LCOV_EXCL_LINE
2509 : }
2510 :
2511 20 : auto &instancesValue = instances.value();
2512 40 : for (auto it = instancesValue.begin(); it != instancesValue.end(); ++it) {
2513 20 : auto const &instance = it.value();
2514 20 : auto const &objName = it.key();
2515 40 : auto const &objNameUC = UtilityRoutines::MakeUPPERCase(objName);
2516 20 : state.dataInputProcessing->inputProcessor->markObjectAsUsed(currObj, objName);
2517 40 : std::shared_ptr<GLHEResponseFactors> thisObj(new GLHEResponseFactors(state, objNameUC, instance));
2518 20 : state.dataGroundHeatExchanger->responseFactorsVector.push_back(thisObj);
2519 : }
2520 : }
2521 :
2522 23 : if (state.dataGroundHeatExchanger->numVertArray > 0) {
2523 :
2524 4 : std::string currObj = "GroundHeatExchanger:Vertical:Array";
2525 :
2526 4 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(currObj);
2527 2 : if (instances == state.dataInputProcessing->inputProcessor->epJSON.end()) {
2528 : ShowSevereError(state, // LCOV_EXCL_LINE
2529 : currObj + ": Somehow getNumObjectsFound was > 0 but epJSON.find found 0"); // LCOV_EXCL_LINE
2530 : }
2531 :
2532 2 : auto &instancesValue = instances.value();
2533 4 : for (auto it = instancesValue.begin(); it != instancesValue.end(); ++it) {
2534 2 : auto const &instance = it.value();
2535 2 : auto const &objName = it.key();
2536 4 : auto const &objNameUC = UtilityRoutines::MakeUPPERCase(objName);
2537 2 : state.dataInputProcessing->inputProcessor->markObjectAsUsed(currObj, objName);
2538 4 : std::shared_ptr<GLHEVertArray> thisObj(new GLHEVertArray(state, objNameUC, instance));
2539 2 : state.dataGroundHeatExchanger->vertArraysVector.push_back(thisObj);
2540 : }
2541 : }
2542 :
2543 23 : if (state.dataGroundHeatExchanger->numSingleBorehole > 0) {
2544 :
2545 0 : std::string currObj = "GroundHeatExchanger:Vertical:Single";
2546 :
2547 0 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(currObj);
2548 0 : if (instances == state.dataInputProcessing->inputProcessor->epJSON.end()) {
2549 : ShowSevereError(state, // LCOV_EXCL_LINE
2550 : currObj + ": Somehow getNumObjectsFound was > 0 but epJSON.find found 0"); // LCOV_EXCL_LINE
2551 : }
2552 :
2553 0 : auto &instancesValue = instances.value();
2554 0 : for (auto it = instancesValue.begin(); it != instancesValue.end(); ++it) {
2555 0 : auto const &instance = it.value();
2556 0 : auto const &objName = it.key();
2557 0 : auto const &objNameUC = UtilityRoutines::MakeUPPERCase(objName);
2558 0 : state.dataInputProcessing->inputProcessor->markObjectAsUsed(currObj, objName);
2559 0 : std::shared_ptr<GLHEVertSingle> thisObj(new GLHEVertSingle(state, objNameUC, instance));
2560 0 : state.dataGroundHeatExchanger->singleBoreholesVector.push_back(thisObj);
2561 : }
2562 : }
2563 :
2564 23 : if (state.dataGroundHeatExchanger->numVerticalGLHEs > 0) {
2565 :
2566 44 : std::string currObj = "GroundHeatExchanger:System";
2567 :
2568 44 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(currObj);
2569 22 : if (instances == state.dataInputProcessing->inputProcessor->epJSON.end()) {
2570 : ShowSevereError(state, // LCOV_EXCL_LINE
2571 : currObj + ": Somehow getNumObjectsFound was > 0 but epJSON.find found 0"); // LCOV_EXCL_LINE
2572 : }
2573 :
2574 22 : auto &instancesValue = instances.value();
2575 44 : for (auto it = instancesValue.begin(); it != instancesValue.end(); ++it) {
2576 22 : auto const &instance = it.value();
2577 22 : auto const &objName = it.key();
2578 44 : auto const &objNameUC = UtilityRoutines::MakeUPPERCase(objName);
2579 22 : state.dataInputProcessing->inputProcessor->markObjectAsUsed(currObj, objName);
2580 22 : state.dataGroundHeatExchanger->verticalGLHE.emplace_back(state, objNameUC, instance);
2581 : }
2582 : }
2583 :
2584 : // SLINKY GLHE
2585 :
2586 23 : if (state.dataGroundHeatExchanger->numSlinkyGLHEs > 0) {
2587 :
2588 2 : std::string currObj = "GroundHeatExchanger:Slinky";
2589 :
2590 2 : auto const instances = state.dataInputProcessing->inputProcessor->epJSON.find(currObj);
2591 1 : if (instances == state.dataInputProcessing->inputProcessor->epJSON.end()) {
2592 : ShowSevereError(state, // LCOV_EXCL_LINE
2593 : currObj + ": Somehow getNumObjectsFound was > 0 but epJSON.find found 0"); // LCOV_EXCL_LINE
2594 : }
2595 :
2596 1 : auto &instancesValue = instances.value();
2597 2 : for (auto it = instancesValue.begin(); it != instancesValue.end(); ++it) {
2598 1 : auto const &instance = it.value();
2599 1 : auto const &objName = it.key();
2600 2 : auto const &objNameUC = UtilityRoutines::MakeUPPERCase(objName);
2601 1 : state.dataInputProcessing->inputProcessor->markObjectAsUsed(currObj, objName);
2602 1 : state.dataGroundHeatExchanger->slinkyGLHE.emplace_back(state, objNameUC, instance);
2603 : }
2604 : }
2605 23 : }
2606 :
2607 : //******************************************************************************
2608 :
2609 23 : void GLHEBase::setupOutput(EnergyPlusData &state)
2610 : {
2611 46 : SetupOutputVariable(state,
2612 : "Ground Heat Exchanger Average Borehole Temperature",
2613 : OutputProcessor::Unit::C,
2614 : this->bhTemp,
2615 : OutputProcessor::SOVTimeStepType::System,
2616 : OutputProcessor::SOVStoreType::Average,
2617 23 : this->name);
2618 46 : SetupOutputVariable(state,
2619 : "Ground Heat Exchanger Heat Transfer Rate",
2620 : OutputProcessor::Unit::W,
2621 : this->QGLHE,
2622 : OutputProcessor::SOVTimeStepType::System,
2623 : OutputProcessor::SOVStoreType::Average,
2624 23 : this->name);
2625 46 : SetupOutputVariable(state,
2626 : "Ground Heat Exchanger Inlet Temperature",
2627 : OutputProcessor::Unit::C,
2628 : this->inletTemp,
2629 : OutputProcessor::SOVTimeStepType::System,
2630 : OutputProcessor::SOVStoreType::Average,
2631 23 : this->name);
2632 46 : SetupOutputVariable(state,
2633 : "Ground Heat Exchanger Outlet Temperature",
2634 : OutputProcessor::Unit::C,
2635 : this->outletTemp,
2636 : OutputProcessor::SOVTimeStepType::System,
2637 : OutputProcessor::SOVStoreType::Average,
2638 23 : this->name);
2639 46 : SetupOutputVariable(state,
2640 : "Ground Heat Exchanger Mass Flow Rate",
2641 : OutputProcessor::Unit::kg_s,
2642 : this->massFlowRate,
2643 : OutputProcessor::SOVTimeStepType::System,
2644 : OutputProcessor::SOVStoreType::Average,
2645 23 : this->name);
2646 46 : SetupOutputVariable(state,
2647 : "Ground Heat Exchanger Average Fluid Temperature",
2648 : OutputProcessor::Unit::C,
2649 : this->aveFluidTemp,
2650 : OutputProcessor::SOVTimeStepType::System,
2651 : OutputProcessor::SOVStoreType::Average,
2652 23 : this->name);
2653 46 : SetupOutputVariable(state,
2654 : "Ground Heat Exchanger Farfield Ground Temperature",
2655 : OutputProcessor::Unit::C,
2656 : this->tempGround,
2657 : OutputProcessor::SOVTimeStepType::System,
2658 : OutputProcessor::SOVStoreType::Average,
2659 23 : this->name);
2660 23 : }
2661 :
2662 : //******************************************************************************
2663 :
2664 288216 : Real64 GLHEVert::calcBHAverageResistance(EnergyPlusData &state)
2665 : {
2666 : // Calculates the average thermal resistance of the borehole using the first-order multipole method.
2667 :
2668 : // Javed, S. & Spitler, J.D. 2016. 'Accuracy of Borehole Thermal Resistance Calculation Methods
2669 : // for Grouted Single U-tube Ground Heat Exchangers.' Applied Energy.187:790-806.
2670 :
2671 : // Equation 13
2672 :
2673 288216 : Real64 const beta = 2 * DataGlobalConstants::Pi * this->grout.k * calcPipeResistance(state);
2674 :
2675 288216 : Real64 const final_term_1 = log(this->theta_2 / (2 * this->theta_1 * pow(1 - pow_4(this->theta_1), this->sigma)));
2676 288216 : Real64 const num_final_term_2 = pow_2(this->theta_3) * pow_2(1 - (4 * this->sigma * pow_4(this->theta_1)) / (1 - pow_4(this->theta_1)));
2677 288216 : Real64 const den_final_term_2_pt_1 = (1 + beta) / (1 - beta);
2678 288216 : Real64 const den_final_term_2_pt_2 = pow_2(this->theta_3) * (1 + (16 * this->sigma * pow_4(this->theta_1)) / pow_2(1 - pow_4(this->theta_1)));
2679 288216 : Real64 const den_final_term_2 = den_final_term_2_pt_1 + den_final_term_2_pt_2;
2680 288216 : Real64 const final_term_2 = num_final_term_2 / den_final_term_2;
2681 :
2682 288216 : return (1 / (4 * DataGlobalConstants::Pi * this->grout.k)) * (beta + final_term_1 - final_term_2);
2683 : }
2684 :
2685 : //******************************************************************************
2686 :
2687 288214 : Real64 GLHEVert::calcBHTotalInternalResistance(EnergyPlusData &state)
2688 : {
2689 : // Calculates the total internal thermal resistance of the borehole using the first-order multipole method.
2690 :
2691 : // Javed, S. & Spitler, J.D. 2016. 'Accuracy of Borehole Thermal Resistance Calculation Methods
2692 : // for Grouted Single U-tube Ground Heat Exchangers.' Applied Energy. 187:790-806.
2693 :
2694 : // Equation 26
2695 :
2696 288214 : Real64 beta = 2 * DataGlobalConstants::Pi * this->grout.k * calcPipeResistance(state);
2697 :
2698 288214 : Real64 final_term_1 = log(pow(1 + pow_2(this->theta_1), this->sigma) / (this->theta_3 * pow(1 - pow_2(this->theta_1), this->sigma)));
2699 288214 : Real64 num_term_2 = pow_2(this->theta_3) * pow_2(1 - pow_4(this->theta_1) + 4 * this->sigma * pow_2(this->theta_1));
2700 288214 : Real64 den_term_2_pt_1 = (1 + beta) / (1 - beta) * pow_2(1 - pow_4(this->theta_1));
2701 288214 : Real64 den_term_2_pt_2 = pow_2(this->theta_3) * pow_2(1 - pow_4(this->theta_1));
2702 288214 : Real64 den_term_2_pt_3 = 8 * this->sigma * pow_2(this->theta_1) * pow_2(this->theta_3) * (1 + pow_4(this->theta_1));
2703 288214 : Real64 den_term_2 = den_term_2_pt_1 - den_term_2_pt_2 + den_term_2_pt_3;
2704 288214 : Real64 final_term_2 = num_term_2 / den_term_2;
2705 :
2706 288214 : return (1 / (DataGlobalConstants::Pi * this->grout.k)) * (beta + final_term_1 - final_term_2);
2707 : }
2708 :
2709 : //******************************************************************************
2710 :
2711 0 : Real64 GLHEVert::calcBHGroutResistance(EnergyPlusData &state)
2712 : {
2713 : // Calculates grout resistance. Use for validation.
2714 :
2715 : // Javed, S. & Spitler, J.D. 2016. 'Accuracy of Borehole Thermal Resistance Calculation Methods
2716 : // for Grouted Single U-tube Ground Heat Exchangers.' Applied Energy. 187:790-806.
2717 :
2718 : // Equation 3
2719 :
2720 0 : return calcBHAverageResistance(state) - calcPipeResistance(state) / 2.0;
2721 : }
2722 :
2723 : //******************************************************************************
2724 :
2725 381344 : Real64 GLHEVert::calcHXResistance(EnergyPlusData &state)
2726 : {
2727 : // Calculates the effective thermal resistance of the borehole assuming a uniform heat flux.
2728 :
2729 : // Javed, S. & Spitler, J.D. Calculation of Borehole Thermal Resistance. In 'Advances in
2730 : // Ground-Source Heat Pump Systems,' pp. 84. Rees, S.J. ed. Cambridge, MA. Elsevier Ltd. 2016.
2731 :
2732 : // Eq: 3-67
2733 :
2734 : using FluidProperties::GetSpecificHeatGlycol;
2735 :
2736 381344 : constexpr const char *RoutineName("calcBHResistance");
2737 :
2738 381344 : if (this->massFlowRate <= 0.0) {
2739 93130 : return 0;
2740 : } else {
2741 864642 : Real64 const cpFluid = GetSpecificHeatGlycol(state,
2742 288214 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2743 : this->inletTemp,
2744 288214 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2745 576428 : RoutineName);
2746 288214 : return calcBHAverageResistance(state) +
2747 288214 : 1 / (3 * calcBHTotalInternalResistance(state)) * pow_2(this->bhLength / (this->massFlowRate * cpFluid));
2748 : }
2749 : }
2750 :
2751 : //******************************************************************************
2752 :
2753 576430 : Real64 GLHEVert::calcPipeConductionResistance()
2754 : {
2755 : // Calculates the thermal resistance of a pipe, in [K/(W/m)].
2756 :
2757 : // Javed, S. & Spitler, J.D. 2016. 'Accuracy of Borehole Thermal Resistance Calculation Methods
2758 : // for Grouted Single U-tube Ground Heat Exchangers.' Applied Energy. 187:790-806.
2759 :
2760 576430 : return log(this->pipe.outDia / this->pipe.innerDia) / (2 * DataGlobalConstants::Pi * this->pipe.k);
2761 : }
2762 :
2763 : //******************************************************************************
2764 :
2765 576432 : Real64 GLHEVert::calcPipeConvectionResistance(EnergyPlusData &state)
2766 : {
2767 : // Calculates the convection resistance using Gnielinski and Petukov, in [K/(W/m)]
2768 :
2769 : // Gneilinski, V. 1976. 'New equations for heat and mass transfer in turbulent pipe and channel flow.'
2770 : // International Chemical Engineering 16(1976), pp. 359-368.
2771 :
2772 : using FluidProperties::GetConductivityGlycol;
2773 : using FluidProperties::GetSpecificHeatGlycol;
2774 : using FluidProperties::GetViscosityGlycol;
2775 :
2776 : // SUBROUTINE PARAMETER DEFINITIONS:
2777 576432 : constexpr const char *RoutineName("calcPipeConvectionResistance");
2778 :
2779 : // Get fluid props
2780 576432 : this->inletTemp = state.dataLoopNodes->Node(this->inletNodeNum).Temp;
2781 :
2782 1729296 : Real64 const cpFluid = GetSpecificHeatGlycol(state,
2783 576432 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2784 : this->inletTemp,
2785 576432 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2786 1152864 : RoutineName);
2787 1729296 : Real64 const kFluid = GetConductivityGlycol(state,
2788 576432 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2789 : this->inletTemp,
2790 576432 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2791 1152864 : RoutineName);
2792 1729296 : Real64 const fluidViscosity = GetViscosityGlycol(state,
2793 576432 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2794 : this->inletTemp,
2795 576432 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2796 1152864 : RoutineName);
2797 :
2798 : // Smoothing fit limits
2799 576432 : constexpr Real64 lower_limit = 2000;
2800 576432 : constexpr Real64 upper_limit = 4000;
2801 :
2802 576432 : Real64 const bhMassFlowRate = this->massFlowRate / this->myRespFactors->numBoreholes;
2803 576432 : Real64 const reynoldsNum = 4 * bhMassFlowRate / (fluidViscosity * DataGlobalConstants::Pi * this->pipe.innerDia);
2804 :
2805 576432 : Real64 nusseltNum = 0.0;
2806 576432 : if (reynoldsNum < lower_limit) {
2807 512976 : nusseltNum = 4.01; // laminar mean(4.36, 3.66)
2808 63456 : } else if (lower_limit <= reynoldsNum && reynoldsNum < upper_limit) {
2809 31660 : Real64 constexpr nu_low = 4.01; // laminar
2810 31660 : Real64 const f = frictionFactor(reynoldsNum); // turbulent
2811 31660 : Real64 const prandtlNum = (cpFluid * fluidViscosity) / (kFluid);
2812 31660 : Real64 const nu_high = (f / 8) * (reynoldsNum - 1000) * prandtlNum / (1 + 12.7 * std::sqrt(f / 8) * (pow(prandtlNum, 2.0 / 3.0) - 1));
2813 31660 : Real64 const sf = 1 / (1 + std::exp(-(reynoldsNum - 3000) / 150.0)); // smoothing function
2814 31660 : nusseltNum = (1 - sf) * nu_low + sf * nu_high;
2815 : } else {
2816 31796 : Real64 const f = frictionFactor(reynoldsNum);
2817 31796 : Real64 const prandtlNum = (cpFluid * fluidViscosity) / (kFluid);
2818 31796 : nusseltNum = (f / 8) * (reynoldsNum - 1000) * prandtlNum / (1 + 12.7 * std::sqrt(f / 8) * (pow(prandtlNum, 2.0 / 3.0) - 1));
2819 : }
2820 :
2821 576432 : Real64 h = nusseltNum * kFluid / this->pipe.innerDia;
2822 :
2823 576432 : return 1 / (h * DataGlobalConstants::Pi * this->pipe.innerDia);
2824 : }
2825 :
2826 : //******************************************************************************
2827 :
2828 63456 : Real64 GLHEVert::frictionFactor(Real64 const reynoldsNum)
2829 : {
2830 : // Calculates the friction factor in smooth tubes
2831 :
2832 : // Petukov, B.S. 1970. 'Heat transfer and friction in turbulent pipe flow with variable physical properties.'
2833 : // In Advances in Heat Transfer, ed. T.F. Irvine and J.P. Hartnett, Vol. 6. New York Academic Press.
2834 :
2835 : // limits picked be within about 1% of actual values
2836 63456 : constexpr Real64 lower_limit = 1500;
2837 63456 : constexpr Real64 upper_limit = 5000;
2838 :
2839 63456 : if (reynoldsNum < lower_limit) {
2840 0 : return 64.0 / reynoldsNum; // pure laminar flow
2841 63456 : } else if (lower_limit <= reynoldsNum && reynoldsNum < upper_limit) {
2842 : // pure laminar flow
2843 63456 : Real64 const f_low = 64.0 / reynoldsNum;
2844 : // pure turbulent flow
2845 63456 : Real64 const f_high = pow(0.79 * log(reynoldsNum) - 1.64, -2.0);
2846 63456 : Real64 const sf = 1 / (1 + exp(-(reynoldsNum - 3000.0) / 450.0)); // smoothing function
2847 63456 : return (1 - sf) * f_low + sf * f_high;
2848 : } else {
2849 0 : return pow(0.79 * log(reynoldsNum) - 1.64, -2.0); // pure turbulent flow
2850 : }
2851 : }
2852 :
2853 : //******************************************************************************
2854 :
2855 576430 : Real64 GLHEVert::calcPipeResistance(EnergyPlusData &state)
2856 : {
2857 : // Calculates the combined conduction and convection pipe resistance
2858 :
2859 : // Javed, S. & Spitler, J.D. 2016. 'Accuracy of Borehole Thermal Resistance Calculation Methods
2860 : // for Grouted Single U-tube Ground Heat Exchangers.' J. Energy Engineering. Draft in progress.
2861 :
2862 : // Equation 3
2863 :
2864 576430 : return calcPipeConductionResistance() + calcPipeConvectionResistance(state);
2865 : }
2866 :
2867 : //******************************************************************************
2868 :
2869 14472 : Real64 GLHESlinky::calcHXResistance(EnergyPlusData &state)
2870 : {
2871 :
2872 : // SUBROUTINE INFORMATION:
2873 : // AUTHOR Matt Mitchell
2874 : // DATE WRITTEN February, 2015
2875 : // MODIFIED
2876 : // RE-ENGINEERED
2877 :
2878 : // PURPOSE OF THIS SUBROUTINE:
2879 : // Calculates the resistance of the slinky HX from the fluid to the
2880 : // outer tube wall.
2881 :
2882 : using FluidProperties::GetConductivityGlycol;
2883 : using FluidProperties::GetDensityGlycol;
2884 : using FluidProperties::GetSpecificHeatGlycol;
2885 : using FluidProperties::GetViscosityGlycol;
2886 :
2887 : // SUBROUTINE PARAMETER DEFINITIONS:
2888 14472 : constexpr const char *RoutineName("CalcSlinkyGroundHeatExchanger");
2889 :
2890 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
2891 : Real64 nusseltNum;
2892 : Real64 Rconv;
2893 14472 : constexpr Real64 A(3150);
2894 14472 : constexpr Real64 B(350);
2895 14472 : constexpr Real64 laminarNusseltNo(4.364);
2896 :
2897 43416 : Real64 cpFluid = GetSpecificHeatGlycol(state,
2898 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2899 : this->inletTemp,
2900 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2901 28944 : RoutineName);
2902 43416 : Real64 kFluid = GetConductivityGlycol(state,
2903 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2904 : this->inletTemp,
2905 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2906 28944 : RoutineName);
2907 43416 : Real64 fluidDensity = GetDensityGlycol(state,
2908 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2909 : this->inletTemp,
2910 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2911 28944 : RoutineName);
2912 43416 : Real64 fluidViscosity = GetViscosityGlycol(state,
2913 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
2914 : this->inletTemp,
2915 14472 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
2916 28944 : RoutineName);
2917 :
2918 : // calculate mass flow rate
2919 14472 : Real64 singleSlinkyMassFlowRate = this->massFlowRate / this->numTrenches;
2920 :
2921 14472 : Real64 pipeInnerRad = this->pipe.outRadius - this->pipe.thickness;
2922 14472 : Real64 pipeInnerDia = 2.0 * pipeInnerRad;
2923 :
2924 14472 : if (singleSlinkyMassFlowRate == 0.0) {
2925 0 : Rconv = 0.0;
2926 : } else {
2927 : // Re=Rho*V*D/Mu
2928 28944 : Real64 reynoldsNum = fluidDensity * pipeInnerDia *
2929 28944 : (singleSlinkyMassFlowRate / fluidDensity / (DataGlobalConstants::Pi * pow_2(pipeInnerRad))) / fluidViscosity;
2930 14472 : Real64 prandtlNum = (cpFluid * fluidViscosity) / (kFluid);
2931 : // Convection Resistance
2932 14472 : if (reynoldsNum <= 2300) {
2933 0 : nusseltNum = laminarNusseltNo;
2934 14472 : } else if (reynoldsNum > 2300 && reynoldsNum <= 4000) {
2935 0 : Real64 sf = 0.5 + 0.5 * std::tanh((reynoldsNum - A) / B);
2936 0 : Real64 turbulentNusseltNo = 0.023 * std::pow(reynoldsNum, 0.8) * std::pow(prandtlNum, 0.35);
2937 0 : nusseltNum = laminarNusseltNo * (1 - sf) + turbulentNusseltNo * sf;
2938 : } else {
2939 14472 : nusseltNum = 0.023 * std::pow(reynoldsNum, 0.8) * std::pow(prandtlNum, 0.35);
2940 : }
2941 14472 : Real64 hci = nusseltNum * kFluid / pipeInnerDia;
2942 14472 : Rconv = 1.0 / (2.0 * DataGlobalConstants::Pi * pipeInnerDia * hci);
2943 : }
2944 :
2945 : // Conduction Resistance
2946 14472 : Real64 Rcond = std::log(this->pipe.outRadius / pipeInnerRad) / (2.0 * DataGlobalConstants::Pi * this->pipe.k) / 2.0; // pipe in parallel so /2
2947 :
2948 14472 : return Rcond + Rconv;
2949 : }
2950 :
2951 : //******************************************************************************
2952 :
2953 64425632 : Real64 GLHEBase::interpGFunc(Real64 const x_val) const
2954 : {
2955 : // Purpose: interpolate between g-function values, with linear extrapolation above and below range
2956 :
2957 64425632 : auto const &x = this->myRespFactors->LNTTS;
2958 64425632 : auto const &y = this->myRespFactors->GFNC;
2959 :
2960 64425632 : auto const &x_begin = x.begin();
2961 64425632 : auto const &x_end = x.end();
2962 64425632 : auto const &upper_it = std::upper_bound(x_begin, x_end, x_val);
2963 :
2964 64425632 : int l_idx = 0;
2965 64425632 : int u_idx = 0;
2966 :
2967 64425632 : if (upper_it == x_begin) {
2968 : // Linear extrapolation beyond the lower bound
2969 476096 : l_idx = 0;
2970 476096 : u_idx = 1;
2971 63949536 : } else if (upper_it == x_end) {
2972 : // Linear extrapolation beyond the upper bound
2973 0 : u_idx = x.size() - 1;
2974 0 : l_idx = u_idx - 1;
2975 : } else {
2976 : // In the middle of the range
2977 63949536 : u_idx = std::distance(x.begin(), upper_it);
2978 63949536 : l_idx = u_idx - 1;
2979 : }
2980 :
2981 64425632 : Real64 const x_low = x[l_idx];
2982 64425632 : Real64 const x_high = x[u_idx];
2983 64425632 : Real64 const y_low = y[l_idx];
2984 64425632 : Real64 const y_high = y[u_idx];
2985 :
2986 64425632 : return (x_val - x_low) / (x_high - x_low) * (y_high - y_low) + y_low;
2987 : }
2988 :
2989 : //******************************************************************************
2990 :
2991 1713112 : Real64 GLHESlinky::getGFunc(Real64 const time)
2992 : {
2993 : // SUBROUTINE INFORMATION:
2994 : // AUTHOR: Matt Mitchell
2995 : // DATE WRITTEN: February, 2015
2996 : // MODIFIED na
2997 : // RE-ENGINEERED na
2998 :
2999 : // PURPOSE OF THIS SUBROUTINE:
3000 : // Gets the g-function for slinky GHXs
3001 : // Note: Base 10 here.
3002 :
3003 1713112 : Real64 LNTTS = std::log10(time);
3004 :
3005 1713112 : return interpGFunc(LNTTS);
3006 : }
3007 :
3008 : //******************************************************************************
3009 :
3010 62712520 : Real64 GLHEVert::getGFunc(Real64 const time)
3011 : {
3012 : // SUBROUTINE INFORMATION:
3013 : // AUTHOR: Matt Mitchell
3014 : // DATE WRITTEN: February, 2015
3015 : // MODIFIED na
3016 : // RE-ENGINEERED na
3017 :
3018 : // PURPOSE OF THIS SUBROUTINE:
3019 : // Gets the g-function for vertical GHXs
3020 : // Note: Base e here.
3021 :
3022 62712520 : Real64 LNTTS = std::log(time);
3023 62712520 : Real64 gFuncVal = interpGFunc(LNTTS);
3024 62712520 : Real64 RATIO = this->bhRadius / this->bhLength;
3025 :
3026 62712520 : if (RATIO != this->myRespFactors->gRefRatio) {
3027 0 : gFuncVal -= std::log(this->bhRadius / (this->bhLength * this->myRespFactors->gRefRatio));
3028 : }
3029 :
3030 62712520 : return gFuncVal;
3031 : }
3032 :
3033 : //******************************************************************************
3034 :
3035 384240 : void GLHEVert::initGLHESimVars(EnergyPlusData &state)
3036 : {
3037 : // SUBROUTINE INFORMATION:
3038 : // AUTHOR: Dan Fisher
3039 : // DATE WRITTEN: August, 2000
3040 : // MODIFIED Arun Murugappan
3041 : // RE-ENGINEERED na
3042 :
3043 : // Using/Aliasing
3044 : using PlantUtilities::RegulateCondenserCompFlowReqOp;
3045 : using PlantUtilities::SetComponentFlowRate;
3046 :
3047 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
3048 768480 : Real64 currTime = ((state.dataGlobal->DayOfSim - 1) * 24 + (state.dataGlobal->HourOfDay - 1) +
3049 768480 : (state.dataGlobal->TimeStep - 1) * state.dataGlobal->TimeStepZone + state.dataHVACGlobal->SysTimeElapsed) *
3050 384240 : DataGlobalConstants::SecInHour;
3051 :
3052 : // Init more variables
3053 :
3054 384240 : if (this->myEnvrnFlag && state.dataGlobal->BeginEnvrnFlag) {
3055 119 : this->initEnvironment(state, currTime);
3056 : }
3057 :
3058 : // Calculate the average ground temperature over the depth of the borehole
3059 384240 : Real64 minDepth = this->myRespFactors->props->bhTopDepth;
3060 384240 : Real64 maxDepth = this->myRespFactors->props->bhLength + minDepth;
3061 384240 : Real64 oneQuarterDepth = minDepth + (maxDepth - minDepth) * 0.25;
3062 384240 : Real64 halfDepth = minDepth + (maxDepth - minDepth) * 0.5;
3063 384240 : Real64 threeQuarterDepth = minDepth + (maxDepth - minDepth) * 0.75;
3064 :
3065 384240 : this->tempGround = 0;
3066 :
3067 384240 : this->tempGround += this->groundTempModel->getGroundTempAtTimeInSeconds(state, minDepth, currTime);
3068 384240 : this->tempGround += this->groundTempModel->getGroundTempAtTimeInSeconds(state, maxDepth, currTime);
3069 384240 : this->tempGround += this->groundTempModel->getGroundTempAtTimeInSeconds(state, oneQuarterDepth, currTime);
3070 384240 : this->tempGround += this->groundTempModel->getGroundTempAtTimeInSeconds(state, halfDepth, currTime);
3071 384240 : this->tempGround += this->groundTempModel->getGroundTempAtTimeInSeconds(state, threeQuarterDepth, currTime);
3072 :
3073 384240 : this->tempGround /= 5;
3074 :
3075 384240 : this->massFlowRate = RegulateCondenserCompFlowReqOp(state, this->plantLoc, this->designMassFlow);
3076 :
3077 384240 : SetComponentFlowRate(state, this->massFlowRate, this->inletNodeNum, this->outletNodeNum, this->plantLoc);
3078 :
3079 : // Reset local environment init flag
3080 384240 : if (!state.dataGlobal->BeginEnvrnFlag) this->myEnvrnFlag = true;
3081 384240 : }
3082 :
3083 : //******************************************************************************
3084 :
3085 119 : void GLHEVert::initEnvironment(EnergyPlusData &state, [[maybe_unused]] Real64 const CurTime)
3086 : {
3087 :
3088 119 : constexpr const char *RoutineName("initEnvironment");
3089 :
3090 119 : this->myEnvrnFlag = false;
3091 :
3092 357 : Real64 fluidDensity = FluidProperties::GetDensityGlycol(state,
3093 119 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
3094 : 20.0,
3095 119 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
3096 238 : RoutineName);
3097 119 : this->designMassFlow = this->designFlow * fluidDensity;
3098 119 : PlantUtilities::InitComponentNodes(state, 0.0, this->designMassFlow, this->inletNodeNum, this->outletNodeNum);
3099 :
3100 119 : this->lastQnSubHr = 0.0;
3101 119 : state.dataLoopNodes->Node(this->inletNodeNum).Temp = this->tempGround;
3102 119 : state.dataLoopNodes->Node(this->outletNodeNum).Temp = this->tempGround;
3103 :
3104 : // zero out all history arrays
3105 119 : this->QnHr = 0.0;
3106 119 : this->QnMonthlyAgg = 0.0;
3107 119 : this->QnSubHr = 0.0;
3108 119 : this->LastHourN = 0;
3109 119 : state.dataGroundHeatExchanger->prevTimeSteps = 0.0;
3110 119 : state.dataGroundHeatExchanger->currentSimTime = 0.0;
3111 119 : this->QGLHE = 0.0;
3112 119 : this->prevHour = 1;
3113 119 : }
3114 :
3115 : //******************************************************************************
3116 :
3117 22 : void GLHEVert::oneTimeInit_new(EnergyPlusData &state)
3118 : {
3119 :
3120 : using PlantUtilities::ScanPlantLoopsForObject;
3121 :
3122 : // Locate the hx on the plant loops for later usage
3123 22 : bool errFlag = false;
3124 22 : ScanPlantLoopsForObject(state, this->name, DataPlant::PlantEquipmentType::GrndHtExchgSystem, this->plantLoc, errFlag, _, _, _, _, _);
3125 22 : if (errFlag) {
3126 0 : ShowFatalError(state, "initGLHESimVars: Program terminated due to previous condition(s).");
3127 : }
3128 22 : }
3129 0 : void GLHEVert::oneTimeInit([[maybe_unused]] EnergyPlusData &state)
3130 : {
3131 0 : }
3132 :
3133 : //******************************************************************************
3134 :
3135 14626 : void GLHESlinky::initGLHESimVars(EnergyPlusData &state)
3136 : {
3137 : // SUBROUTINE INFORMATION:
3138 : // AUTHOR: Dan Fisher
3139 : // DATE WRITTEN: August, 2000
3140 : // MODIFIED Arun Murugappan
3141 : // RE-ENGINEERED na
3142 :
3143 : // Using/Aliasing
3144 : using PlantUtilities::RegulateCondenserCompFlowReqOp;
3145 : using PlantUtilities::SetComponentFlowRate;
3146 : using namespace GroundTemperatureManager;
3147 :
3148 29252 : Real64 CurTime = ((state.dataGlobal->DayOfSim - 1) * 24 + (state.dataGlobal->HourOfDay - 1) +
3149 29252 : (state.dataGlobal->TimeStep - 1) * state.dataGlobal->TimeStepZone + state.dataHVACGlobal->SysTimeElapsed) *
3150 14626 : DataGlobalConstants::SecInHour;
3151 :
3152 : // Init more variables
3153 :
3154 14626 : if (this->myEnvrnFlag && state.dataGlobal->BeginEnvrnFlag) {
3155 6 : this->initEnvironment(state, CurTime);
3156 : }
3157 :
3158 14626 : this->tempGround = this->groundTempModel->getGroundTempAtTimeInSeconds(state, this->coilDepth, CurTime);
3159 :
3160 14626 : this->massFlowRate = RegulateCondenserCompFlowReqOp(state, this->plantLoc, this->designMassFlow);
3161 :
3162 14626 : SetComponentFlowRate(state, this->massFlowRate, this->inletNodeNum, this->outletNodeNum, this->plantLoc);
3163 :
3164 : // Reset local environment init flag
3165 14626 : if (!state.dataGlobal->BeginEnvrnFlag) this->myEnvrnFlag = true;
3166 14626 : }
3167 :
3168 : //******************************************************************************
3169 :
3170 6 : void GLHESlinky::initEnvironment(EnergyPlusData &state, Real64 const CurTime)
3171 : {
3172 :
3173 6 : constexpr const char *RoutineName("initEnvironment");
3174 :
3175 6 : this->myEnvrnFlag = false;
3176 :
3177 18 : Real64 fluidDensity = FluidProperties::GetDensityGlycol(state,
3178 6 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidName,
3179 : 20.0,
3180 6 : state.dataPlnt->PlantLoop(this->plantLoc.loopNum).FluidIndex,
3181 12 : RoutineName);
3182 6 : this->designMassFlow = this->designFlow * fluidDensity;
3183 6 : PlantUtilities::InitComponentNodes(state, 0.0, this->designMassFlow, this->inletNodeNum, this->outletNodeNum);
3184 :
3185 6 : this->lastQnSubHr = 0.0;
3186 6 : state.dataLoopNodes->Node(this->inletNodeNum).Temp = this->groundTempModel->getGroundTempAtTimeInSeconds(state, this->coilDepth, CurTime);
3187 6 : state.dataLoopNodes->Node(this->outletNodeNum).Temp = this->groundTempModel->getGroundTempAtTimeInSeconds(state, this->coilDepth, CurTime);
3188 :
3189 : // zero out all history arrays
3190 6 : this->QnHr = 0.0;
3191 6 : this->QnMonthlyAgg = 0.0;
3192 6 : this->QnSubHr = 0.0;
3193 6 : this->LastHourN = 0;
3194 6 : state.dataGroundHeatExchanger->prevTimeSteps = 0.0;
3195 6 : state.dataGroundHeatExchanger->currentSimTime = 0.0;
3196 6 : this->QGLHE = 0.0;
3197 6 : this->prevHour = 1;
3198 6 : }
3199 :
3200 : //******************************************************************************
3201 :
3202 1 : void GLHESlinky::oneTimeInit_new(EnergyPlusData &state)
3203 : {
3204 : using PlantUtilities::ScanPlantLoopsForObject;
3205 :
3206 : // Locate the hx on the plant loops for later usage
3207 1 : bool errFlag = false;
3208 1 : ScanPlantLoopsForObject(state, this->name, DataPlant::PlantEquipmentType::GrndHtExchgSlinky, this->plantLoc, errFlag, _, _, _, _, _);
3209 1 : if (errFlag) {
3210 0 : ShowFatalError(state, "initGLHESimVars: Program terminated due to previous condition(s).");
3211 : }
3212 1 : }
3213 0 : void GLHESlinky::oneTimeInit([[maybe_unused]] EnergyPlusData &state)
3214 : {
3215 0 : }
3216 :
3217 : //******************************************************************************
3218 :
3219 2313 : } // namespace EnergyPlus::GroundHeatExchangers
|
