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