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