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