GHE g-function calculation enhancements

.github/workflows/linux_build.yml

@@ -69,6 +69,7 @@ jobs:
                 -DUSE_PSYCH_STATS=ON \
                 -DUSE_PSYCH_ERRORS=OFF \
                 -DENABLE_PCH=OFF \
+                -DBUILD_CPGFUNCTION=ON \
                 ../
 
       - name: Build EnergyPlus

CMakeLists.txt

@@ -265,6 +265,9 @@ target_include_directories(project_options SYSTEM INTERFACE "${kiva_SOURCE_DIR}/
 add_subdirectory(scripts/dev/generate_embeddable_epJSON_schema)
 set_target_properties(generate_embeddedable_epJSON_schema PROPERTIES FOLDER "Internal")
 
+target_include_directories(project_options INTERFACE ${PROJECT_SOURCE_DIR}/third_party/cpgfunctionEP-0.3.2)
+target_include_directories(project_options INTERFACE ${PROJECT_SOURCE_DIR}/third_party/cpgfunctionEP-0.3.2/include)
+
 if(OPENGL_FOUND)
   set(BUILD_PENUMBRA_TESTING
       OFF

doc/engineering-reference/src/simulation-models-encyclopedic-reference-002/heat-exchangers.tex

@@ -929,7 +929,11 @@ \subsubsection{Long Time-Step Response Factors}\label{long-time-step-response-fa
 
 The \emph{g}-functions developed by Eskilson are given only down to values of $\ln t/t_s = -4.5$, which for typical system and soil types can be on the order of 100 days. The values can be extended by the applying a line source model down to a time of $t = 5r_b^2/\alpha$. This time varies from 3-6 hours for a typical borehole field. This is because the analytical line source model, based on which the Eskilson model was developed, does not give a prompt increase in borehole wall temperature at \(r = {r_b}\) . It gives acceptable results only after the non-dimensional times of \(\alpha t/r_b^2 > 5\). But to model short time responses of a borehole we need response factors which can give accurate results down to minutes.
 
-In order to generate the long time-step response factors on the fly, EnergyPlus uses the model developed by Marcotte \& Pasquier(2009) which uses a discretized line source model. The \emph{g}-functions are generated using the following equation. The boreholes are discretized into segments. The temperature response of each segment on all other segments is then used to determine response factor for that particular geometry. The model estimates surface effects by creating ``imaginary" boreholes which are mirrored about the ground surface.
+In order to generate the long time-step response functions on the fly, EnergyPlus has two different models that make two different approximations for the calculation of the response functions. The first model uses uniform heat flux boundary conditions (UHFcalc model) – that is, it assumes that all boreholes have the same uniform heat flux. The heat flux may vary over time, but it is always uniform throughout the field. The second model uses uniform borehole wall temperatures (UBHWTcalc model) – the heat flux will vary between boreholes and vary vertically for each borehole.
+
+Calculation of \emph{g}-functions with uniform heat flux boundary conditions is simpler, but accuracy can decrease as the borehole-to-borehole interference increases. Consider a case with a large rectangular borefield, say 9x16 boreholes – if there is a significant annual imbalance in the heat rejection/heat extraction, the interior of the field will become thermally saturated, and the heat flux in the interior of the field will decrease over time. See Spitler, et al. (2020) for an illustration. The accuracy of using this approximation is discussed by Malayappan and Spitler (2013) with application to sizing of ground heat exchangers. As shown there, use of the UHFcalc model tends to overpredict the long-term temperature change for fields with significant borehole-to-borehole interference and significant annual heat rejection/heat extraction. In general, the UBHWTcalc model can be recommended for all applications, even though for a smaller numbers of boreholes the UHFcalc model offers sufficient accuracy.
+
+For the UHFcalc model, EnergyPlus builds the model developed by Marcotte \& Pasquier(2009) which uses a discretized line source model. The \emph{g}-functions are generated using the following equation. The boreholes are discretized into segments. The temperature response of each segment on all other segments is then used to determine response factor for that particular geometry. The model estimates surface effects by creating ``imaginary" boreholes which are mirrored about the ground surface.
 
 \begin{equation}
 g = \frac{1}{2 H_T} \sum_{i=1}^N \sum_{j=1}^N \left( \int_{u_1^i}^{u_2^i} \int_{u_1^i}^{u_2^i} \frac{erfc\left( \frac{d(u_i, u_j)}{2\sqrt{\alpha t}}\right)}{d(u_i, u_j)} - \frac{erfc\left( \frac{d(u_i, u_j^{'})}{2\sqrt{\alpha t}}\right)}{d(u_i, u_j^{'})}\right)
@@ -953,6 +957,11 @@ \subsubsection{Long Time-Step Response Factors}\label{long-time-step-response-fa
 
 \(d(u_i, u_j^{'})\) is the distance between current point and point on other imaginary borehole.
 
+For the UBHWTcalc model, EnergyPlus builds on the model developed by Cimmino (2018a, 2018b, 2019), with improvements (Cook and Spitler 2021) to reduce memory consumption and increase speed. The general approach is similar to that of the UHFcalc model, but the heat flux is adjusted for each segment to enforce the uniform borehole wall temperature model. Additionally, the methodology of Cimmino (UBHWTcalc) is a semi-analytical solution, so the solution is dependent on the number of segments used. The number of segments are adaptively discretized to approximate a uniform inlet fluid temperature (UIFT) \emph{g}-function. 
+
+Cimmino (2015) introduced a UIFT \emph{g}-function calculation that most closely matches physical reality. As boreholes are usually plumbed in parallel, the entering fluid temperatures are approximately uniform, and the actual distribution of heat within the field then depends on both the heat transfer outside the borehole and the heat transfer within the borehole. Thus, the g-function calculated with the UIFT boundary conditions depends on the flow rate and borehole thermal resistance in addition to the geometric configuration of the boreholes. Like, the UBHWT \emph{g}-function calculation, the UIFT calculation also depends on the number of segments used – as the number of segments used increases, the solution converges. 
+
+So, though it’s possible to find a converged solution for the UIFT \emph{g}-function, it strictly speaking only applies when the borehole thermal resistance and flow rate remain fixed. Both commonly change during the system operation, as changing fluid temperatures affect the fluid properties. Therefore, the UIFT solution might be thought of as having a range of g-functions that change slightly with time. Therefore, in much the same way that radiation heat transfer within a zone is modeled with simplified methods because the occupants are likely to move the furniture, the UBHWT calculation can be used as a reasonable approximation to the UIFT calculation. As shown by Spitler, et al. (2020) the UIFT \emph{g}-function can closely match the UIFT \emph{g}-function, with considerably fewer segments. The number of segments required to closely approximate the UIFT \emph{g}-function varies with the number of boreholes and the borehole depth. Therefore, the UBHWT model uses an adaptive discretization algorithm that uses a smaller number of segments than the UIFT calculation. This provides excellent accuracy while requiring less computational time and memory. 
 
 \subsubsection{Short Time-Step Response Factors}\label{short-time-step-response-factors}
 
@@ -1129,15 +1138,57 @@ \subsubsection{Summary of Variable Short Time Step Response Factor Model}\label{
 
 \subsubsection{References}\label{references-2-006}
 
-Eskilson, P. 1987. Thermal Analysis of Heat Extraction Boreholes. Ph.D.~Thesis, Department of Mathematical Physics, University of Lund, Lund, Sweden.
+\hangindent=2em
+\hangafter=1
+\noindent Claesson, J., G. Helstr{\"o}m. 2011. Multipole method to calculate borehole thermal resistances in a borehole heat exchanger. HVAC\&R Research. 17(6), 895-911.
+
+\hangindent=2em
+\hangafter=1
+\noindent Cook, J. C. and J. D. Spitler. 2021. Faster computation of g-functions used for modeling of ground heat exchangers with reduced memory consumption. Building Simulation 2021. Bruges, Belgium, IBPSA.
+
+\hangindent=2em
+\hangafter=1
+\noindent Cimmino, M. (2015). "The effects of borehole thermal resistances and fluid flow rate on the g-functions of geo-thermal bore fields." International Journal of Heat and Mass Transfer 91: 1119-1127.
+
+\hangindent=2em
+\hangafter=1
+\noindent Cimmino, M. 2018a. "Fast calculation of the g-functions of geothermal borehole fields using similarities in the evaluation of the finite line source solution." Journal of Building Performance Simulation 11(6): 655-668.
+
+\hangindent=2em
+\hangafter=1
+\noindent Cimmino, M. 2018b. pygfunction: an open-source toolbox for the evaluation of thermal. eSim 2018, Montreál, IBPSA Canada.
+
+\hangindent=2em
+\hangafter=1
+\noindent Cimmino, M. 2019. "Semi-Analytical Method for g-Function Calculation of bore fields with series- and parallel-connected boreholes." Science and Technology for the Built Environment 25(8): 1007-1022.
+
+\hangindent=2em
+\hangafter=1
+\noindent Eskilson, P. 1987. Thermal Analysis of Heat Extraction Boreholes. Ph.D.~Thesis, Department of Mathematical Physics, University of Lund, Lund, Sweden.
+
+\hangindent=2em
+\hangafter=1
+\noindent Malayappan, V. and J. D. Spitler. 2013. Limitations of Using Uniform Heat Flux Assumptions in Sizing Vertical Borehole Heat Exchanger Fields. Clima 2013. Prague (Czech Republic).
+
+\hangindent=2em
+\hangafter=1
+\noindent Spitler, J. D., J. C. Cook and X. Liu. 2020. A Preliminary Investigation on the Cost Reduction Potential of Optimizing Bore Fields for Commercial Ground Source Heat Pump Systems. Proceedings, 45th Workshop on Geothermal Reservoir Engineering. Stanford, California, Stanford University. 
 
-Claesson, J., G. Helstr{\"o}m. 2011. Multipole method to calculate borehole thermal resistances in a borehole heat exchanger. HVAC\&R Research. 17(6), 895-911.
+\hangindent=2em
+\hangafter=1
+\noindent Spitler, J. D., J. C. Cook and X. Liu. 2020. FY20 Second Milestone Report for Advanced Techno-Economic Model-ing for Geothermal Heat Pump Applications in Residential, Commercial, and Industry Building, Oak Ridge Na-tional Laboratory.
 
-Yavuzturk, C. 1999. Modeling of Vertical Ground Loop Heat Exchangers for Ground Source Heat Pump Systems. Ph.D.~Thesis, Department of Mechanical and Aerospace Engineering, Oklahoma State University, Stillwater, Oklahoma.
+\hangindent=2em
+\hangafter=1
+\noindent Xu, X., J.D. Spitler. 2006. Modeling of Vertical Ground Loop Heat Exchangers with Variable Convective and Thermal Mass of Fluid. Proceedings of the 10th International Conference on Thermal Energy Storage-Ecostock 2006, Pomona, NJ.
 
-Yavuzturk, C., J.D. Spitler. 1999. A Short Time Step Response Factor Model for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2):475-485.
+\hangindent=2em
+\hangafter=1
+\noindent Yavuzturk, C. 1999. Modeling of Vertical Ground Loop Heat Exchangers for Ground Source Heat Pump Systems. Ph.D.~Thesis, Department of Mechanical and Aerospace Engineering, Oklahoma State University, Stillwater, Oklahoma.
 
-Xu, X., J.D. Spitler. 2006. Modeling of Vertical Ground Loop Heat Exchangers with Variable Convective and Thermal Mass of Fluid. Proceedings of the 10th International Conference on Thermal Energy Storage-Ecostock 2006, Pomona, NJ.
+\hangindent=2em
+\hangafter=1
+\noindent Yavuzturk, C., J.D. Spitler. 1999. A Short Time Step Response Factor Model for Vertical Ground Loop Heat Exchangers. ASHRAE Transactions. 105(2):475-485.
 
 \subsection{GroundHeatExchanger:Slinky}\label{groundheatexchangerslinky}
 

doc/input-output-reference/src/overview/group-condenser-equipment.tex

@@ -3026,6 +3026,10 @@ \subsubsection{Inputs}\label{inputs-10-002}
 
 The unique name of the \hyperref[groundheatexchangerresponsefactors]{GroundHeatExchanger:ResponseFactors} object used to define the third-party response factors. If present, the GHE:Vertical:Array and GHE:Vertical:Single objects defined will be ignored.
 
+\paragraph{Field: g-Function Calculation Method}
+
+If the GHE:ResponseFactors object is not provided, this input will be relevant to the calculation method used to generate the g-function values. There are two options for this input: ``UHFcalc" or ``UBHWTcalc". The default is UHFcalc, where a g-function is computed with the uniform heat flux boundary condition. The UBHWTcalc option computes the g-function with a uniform borehole wall temperature boundary condition. The g-function computed with a uniform borehole wall temperature boundary condition is expected to be a better approximation of the borehole physics and is also expected to be faster to compute.
+
 \paragraph{Field: GHE:Vertical:Array object name}
 
 The unique name of the \hyperref[groundheatexchangerverticalarray]{GroundHeatExchanger:Vertical:Array} object used to define a rectangular borehole field. If present, the GHE:Vertical:Single objects will be ignored.
@@ -3047,7 +3051,8 @@ \subsubsection{Inputs}\label{inputs-10-002}
     2.5,                !- Ground Thermal Conductivity {W/m-K}
     2.5E+06,            !- Ground Thermal Heat Capacity {J/m3-K}
     ,                   !- Response Factors Object Name
-    GHE-Array;
+    UBHWTcalc,          !- g-Function Calculation Method
+    GHE-Array;          !- GHE:Vertical:Array Object Name
 
   GroundHeatExchanger:Vertical:Properties,
     GHE-1 Props,        !- Name
@@ -3083,6 +3088,7 @@ \subsubsection{Inputs}\label{inputs-10-002}
     2.5,                !- Ground Thermal Conductivity {W/m-K}
     2.5E+06,            !- Ground Thermal Heat Capacity {J/m3-K}
     ,                   !- Response Factors Object Name
+    UBHWTcalc,          !- g-Function Calculation Method
     ,                   !- GHE Array Object Name
     GHE-1,              !- GHE Borehole Definition 1
     GHE-2,              !- GHE Borehole Definition 2