Control of Heat Transfer in a Vertical Ground Heat Exchanger for a Geothermal Heat Pump System

Abstract

This paper presents a mathematical model of heat transfer behavior between the liquid inside vertical underground geothermal pipes and the surrounding ground for heating (in the winter) and cooling (in the summer) modes in a ground heat exchanger (GHE) that can optimize its output temperature. The GHE’s output temperature reaches the appropriate value when the water velocity is lowered enough. Subsequently, the proposed model was applied to a case study of a 400-ton geothermal heat pump system (GHPS) at Oakland University, in both the heating and cooling modes, to assess its validity and improve the GHE’s performance. The model was implemented in MATLAB using an ordinary differential equation (ODE) solver. Four different water velocities were used to demonstrate the significant effect of velocity on the loop exit temperature. Model predictive control (MPC) was designed to optimize the GHE’s output temperature by controlling the water velocity, which could reduce the energy consumption used for heat and water circulating pumps. The results reveal that the acceptable range of the water velocity for Oakland University’s GHE was between 0.35 and 0.45 $\mathrm{m}/\mathrm{s}$, which ensured that the heat pump system delivered the proper temperature to provide the Human Health Building (HHB) with a comfortable temperature regardless of the season. The suggested water velocity ranges in vertical single U-tube pipes with diameters of De 25 $\mathrm{mm}$, De 32 $\mathrm{mm}$, and De 40 $\mathrm{mm}$ are between 0.33 and 0.43 $\mathrm{m}/\mathrm{s}$, 0.35 to 0.45 $\mathrm{m}/\mathrm{s}$, and 0.38 to 0.48 $\mathrm{m}/\mathrm{s}$, respectively.

1. Introduction

To take advantage of the existing heat contained within the earth, an efficient technology known as the geothermal heat pump system (GHPS) was proposed. GHPSs were developed to efficiently heat and cool buildings by harnessing this energy from the earth. Thus, the GHPS provides high-efficiency heating and cooling for buildings, as well as domestic hot water production that can be used in different applications (electricity production, greenhouse, etc.). The heating and cooling efficiencies of the GHPS outperform conventional heating and air conditioning systems, which are 30% to 70% higher than conventional heating, and 20% to 50% higher than air conditioning systems. Furthermore, the GHPS significantly reduces electricity use by approximately 75% in compassion to conventional heating or cooling systems [1]. It would be best if renewable energy sources were used to meet GHPSs’ power needs instead of electricity. Energy efficiency attained with this GHPS technology will have significant environmental benefits, such as decreased carbon dioxide (${\mathrm{CO}}_2$) emissions and a quieter atmosphere. Therefore, the GHPS is considered a clean, reliable, renewable, and sustainable source of energy that requires very low maintenance and is unaffected by climatic conditions. It is also environmentally friendly and greatly reduces ${\mathrm{CO}}_2$ emissions by up to 52% when compared to heating, ventilation, and air conditioning (HVAC) [1,2,3,4,5,6]. After incurring initial installation costs, GHPS is free to operate. Since there are few restrictions for GHPS applications, it can be installed virtually anywhere and operated continuously for both heating and cooling purposes to provide the building with a comfortable temperature regardless of the weather. Additionally, as the GHPS benefits from more heat energy, it would be best to install its ground heat exchanger (GHE) in an area with high thermal conductivity, preferably moist soil. The GHPS was first used in North America less than one hundred years ago and is now gaining popularity globally, especially in colder regions, because of its high energy efficiency and thus lower operating cost compared to conventional systems [7]. Additionally, countries utilizing GHPS aim to create environments with a low-carbon economy by reducing ${\mathrm{CO}}_2$ emissions. It should be highlighted that regular monitoring and control are necessary to keep the system operating at optimum efficiency. On the other hand, despite all the benefits of the GHPS, its installation is expensive due to the high cost of drilling boreholes if the design is vertical, and it requires a large land area if the design is horizontal [8]. In addition, its initial capital cost is significant, roughly 30–50% higher than traditional heating and air conditioning systems [9]. For these reasons, the GHPS is not used as widely as it could be if it were cheaper to install. To compensate for the high installation costs, researchers have conducted numerous studies to increase the GHPS’s efficiency and to make the system more reliable and affordable for household use. Ultimately, due to its high efficiency, the GHPS pays for itself in the long term.

The ground heat exchanger (GHE) is a crucial component of the GHPS and plays a key role in utilizing the earth’s heat energy by exchanging heat with the ground due to its actual contact. Most benefits from the earth’s thermal energy can be attained when the GHE is well-designed, including in terms of configuration, installation location, borehole heat exchanger capacity, and performance factors, which ensure that the GHE’s output temperature reaches the ground temperature. Thus, improvements should emphasize GHE characteristics, as any improvement will benefit from the overall GHPS efficiency and result in a lower initial capital cost.

2. GHE Performance Factors

Many factors can affect the vertical GHE’s performance, such as the soil type, backfill material thermal conductivity, pipe’s material thermal conductivity, borehole depth, the distance between two boreholes, shank spacing between the pipes, and working fluid and its velocity in the U-tube pipe, as shown Figure 1.

For instance, each soil type has a different thermal conductivity that is determined by its physical component properties. When the thermal conductivity of the soil material is high, it indicates that it is fully saturated, and its thermal resistance is therefore low, allowing heat to flow through its layers quickly. Conversely, when its thermal conductivity is low, it indicates that the soil material is unsaturated, and its thermal resistance increases, slowing the heat flow. Therefore, its thermal conductivity is dependent on its saturation level [10]. Wang et al. [11] evaluated the thermal conductivity of different soil types, including loam, clay, and sand, using the improved ‘de V-1’ model. The results revealed that the sand performed better than loam and clay soils in terms of thermal conductivity and heat diffusion. Moreover, Wang et al. [11] agreed with Tang and Hossein [12] that sand performed better than clay. Therefore, the soil’s high thermal conductivity is essential for use in either installing the GHE or as a backfill material to have a suitable design for GHE.

Following the installation of the U-tube pipe, the space between the pipe and the borehole wall is filled with conventional backfill materials, including bentonite mixed with sand and cement with the amount of water which depends on the mass ratio of each material component. These materials have the advantages of being extremely durable, swelling-capable, and impermeable, but their thermal conductivities are low. Depending on the type of surrounding soil, engineers vary the ratio of each material in the mixture to best suit the environment. Therefore, the boreholes must be completely grouted to ensure that the grout fills the borehole without any air gaps, which can cause thermal discontinuity (called contact resistance). That will create a strong thermal connection between the U-tube GHE pipe and the ground, which increases the heat transfer capacity. It also keeps the borehole wall from collapsing [13].

To maintain the GHE’s efficiency in coastal regions, Lee et al. [14] suggested utilizing bentonite that contains a substantial quantity of montmorillonite. Furthermore, Lee et al. [15] recommended using cementitious backfill material in dry soil due to its good desiccation and shrinkage properties. Shadley and Den Braven [16] tested several backfill mixtures and pure materials and then compared them to a standard bentonite backfill material. Experimental results of the backfill material test showed that the highest thermal conductivity was found to be in saturated natural sandy soil.

Table 1. Properties of some soil materials.

Parameter	Density (kg/m3)	Specific Heat Capacity (J/(kg K))	Thermal Conductivity (W/m·K)
Sand	1500	1710	1.13
Clay	1700	1800	1.2
Sandy-silt	1847	1200	1.3
Sandy-clay	1960	1200	2.1

demonstrates the soil material properties [17].

The backfill material is regarded as a heat transfer medium between the GHE and the soil; thus, it has a high thermal conductivity to fulfill heat transfer into directions. The pure backfill material will have a thermal conductivity ranging from 0.8 W/m·K to 2.4 W/m·K, but when mixed with other soils, it could be as high as 2.6 to 3.5 W/m·K [14,18,19]. Dehkordi and Schincariol [20] demonstrated that increasing the thermal conductivity of backfill material has a positive impact on heat extraction. Thus, they pointed out that by using the high thermal conductivity of the backfill material (3 W/m·K), heat extraction can be increased by more than 10% in comparison to the low thermal conductivity of the backfill material (1 W/m·K). Additionally, increasing the backfill material’s thermal conductivity could reduce the thermal gradient between the pipe and the surrounding ground, which increases the efficiency of the heat exchange rate [20]. Alberti et al. [21] compared the thermal conductivities of three different backfill materials. The results showed that the thermal resistance of the first and second boreholes backfilled with 2.3 W/m·K was lower than the thermal resistance of the third borehole backfilled with 0.7 W/m·K, with thermal resistance result values of 0.054 m·K/W, 0.054 m·K/W, and 0.1354 m·K/W, respectively. For instance, the use of high thermal conductivity backfill materials could increase the borehole surface temperature and reduce thermal borehole resistance [22]. This results in rapid heat transfer between the fluid inside the pipe and its surrounding grounds and can shorten the length of the GHE pipe. Thus, the GHPS’s initial capital and operating costs are reduced. This becomes more important in the double U-tube borehole configuration than in the single U-tube borehole configuration. Zhang et al. [23] demonstrated the influences of the backfill material’s thermal conductivity on the GHE’s four different U-tube pipe sizes and configurations, namely as double U-tube 32, double U-25, single U-32, and single U-25. They concluded that the heat thermal conductivity of the backfill material has the highest influence on reducing the GHE pipe’s length in a double U-32, followed by double U-25, single U-32, and single U-25.

The pipe acts as a heat transfer interface between two mediums, namely, the working fluid inside the pipe, and its surrounding grounds. The pipe should have durability, flexibility, and longevity to resist wear, leaks, and other damage; to be simple to install; and for cost-effectiveness [24]. Furthermore, the pipe’s material should have a high thermal conductivity to promote the heat transfer process. It is important to note that the thermal conductivity of the pipe material influences heat transfer performance, which in turn can enable heat to flow in either direction. Most geothermal ground heat pump system projects use a high-density polyethylene (HDPE) pipe because of its durability, flexibility, longevity, and low cost, despite its conductivity being lower than metal pipe. The pipe’s material thermal conductivity ranges from 0.35 W/m·K to 0.49 W/m·K, and its lifespan is up to approximately 50 years [25]. Therefore, more improvements in manufacturing characteristics are needed to increase its thermal conductivity, which could improve the GHE’s performance. For instance, Bassiouny et al. [26] mixed miniature aluminum wires with a conventional HDPE pipe material. They found that when an additional 2 mm of aluminum wire was added to the material of the pipe itself, the thermal conductivity rose 25%; when an additional 3 mm of aluminum wire was added to the material of the pipe itself, the thermal conductivity rose 150%, when compared to the original HDPE pipe. In addition, Raymond et al. [27] incorporated some additives into the polymer resin in the original HDPE pipe extrusion. The result demonstrated that the enhanced pipe material’s thermal conductivity increased by 75%, from 0.4 W/m·K to 0.7 W/m·K, as compared to the original HDPE pipe. The results also showed that this improvement reduced the thermal resistance of the borehole by 24%. Moreover, Bouhacina et al. [28] added longitudinal fins to the original HDPE pipe’s inner surface. The new design pipe was able to increase heat extraction by approximately 7% in comparison to the original pipe. Therefore, increasing the thermal conductivity of the pipe material could bring some advantages, such as reducing borehole thermal resistance, decreasing pipe length, increasing heat extraction from the ground, and either heating or cooling the water leading to lower initial installation costs.

As the borehole gets deeper, more heat transfer occurs until it reaches a point where it loses more heat than it gains. At that point, when the difference between the ground and fluid temperature decreases, the heat transfer efficiency will relatively decrease, and the borehole drilling cost will significantly rise. In the GHPS’s engineering design, increasing the borehole depth does not bring greater benefits than improving other parameters such as increasing the thermal conductivity of backfill material; this can be satisfied with the existing borehole depth if it provides adequate heating energy. Moreover, increasing the borehole depth lengthens the period of water circulation, requiring more energy to circulate the water within the loop. Given this, Song et al. [29] suggested that to limit heat loss, the borehole should not be too deep; thus, the depth for a single U-tube would be roughly 80 m. Furthermore, the depth of the borehole for double U-25 is between 80 m and 100 m, whereas the depth of the borehole for double U-32 is between 90 m and 110 m [8]. Zeng et al. [30] compared two borehole thermal resistances for double U-tube and single U-tube borehole configurations. The result showed that the double U-tube borehole had better performance than a single U-tube borehole since its thermal resistance was reduced by 30–90% when compared to the single U-tube borehole thermal resistance. Therefore, more heat transfer efficiency in the GHE can be attained when the double U-tube configuration is connected in parallel rather than in series [30]. Furthermore, the multi-U-tube borehole configuration outperforms the double and single U-tube borehole configurations [31].

It is also worth noting that the distance between two boreholes is crucial, as it has a significant impact on the borehole’s thermal resistance. Enough space must be left between them to prevent any thermal interference as it affects the performance of the GHE. The thermal interference is measured by the distance between two boreholes; if the distance is less than what is suggested for the GHE pipe installation design, there will be more thermal interference between them. However, if the distance is too large, more land area is required to install the GHE, and an additional pipe network is required for connections between boreholes pipes, which will increase the initial capital cost. Therefore, this distance should be neither too small nor too large. Gultekin et al. [32] concluded that the distance between two U-tube boreholes should be between 4.5 m and 6 m. Moreover, they agreed with another study [33] that 5 m is the optimum separation between two borehole walls.

Another factor that influences heat transfer performance in the GHE is the shank spacing between two legs of U-tube pipes, which causes a thermal short circuit phenomenon when the legs of the U-tube pipes are in contact or extremely close to each other [34]. If there is not enough space, the probability of thermal short-circuiting increases, making heat transfer inefficient. To prevent a thermal short circuit and enhance heat transfer in the GHE between the U-tube pipe and its surroundings, the U-tube pipe must be centrally located in the borehole with precise space between them. Interestingly, shank separation of less than 60 mm between U-tube pipes will cause a thermal short circuit, and the heat transfer deteriorates [35]. However, if the space is raised to 60 mm, a thermal short circuit will not occur [36]. In addition, it could be stable if the shank spacing is larger than 90 mm [12]. Furthermore, Zheng et al. [33] determined that the ideal spacing between two U-tube pipes is between 100 mm and 200 mm.

The working fluid acts as a heat transfer carrier between the heat pump and the ground; it has high thermal conductivity, heat capacity, and the lowest viscosity because it affects the convective thermal resistance (between the fluid in the U-tube and the ground) [37]. Dada and Benchatti [38] evaluated three different working fluids—water, gasoline, and glycol. According to their findings, utilizing water as a working fluid could store and recover more thermal heat energy than other working fluids. Pure water has been employed as a heat-carrying fluid in many GHPS systems since their inception for cooling and heating processes due to its high ability to absorb and store heat from the ground and low cost, but its applications have limited antifreeze properties. In cold climates, pure water could be mixed with antifreeze to prevent it from freezing, which would improve heat transfer in GHE. When antifreeze is mixed with the working fluid, the volume concentration must be exact. If the concentration is not in the proper proportion, it could affect the convective thermal resistance, which in turn can affect the GHE’s heat transfer performance. For instance, adding 33% ethanol to 67% pure water increases the kinematic viscosity, decreases the water’s thermal conductivity, and increases the convective thermal resistance, affecting the water–ground heat exchange. However, if the proportion is reduced to 20%, it could enhance the working fluid’s physical properties, maximizing the heat transfer between water and its surrounding grounds [39]. Tang and Nowamooz [12] concluded that adding 24% ethanol to the working fluid outperformed adding 20% CaCl₂, 25% propylene glycol, and 33% propylene glycol, respectively, to the working fluid. Recently, GHPS’s engineer designers have turned to using nanofluids as working fluid due to the obvious improvement in the convective heat transfer coefficient [40,41]. Narei et al. [42] concluded that employing an Al₂O₃/water nanofluid as a working fluid could decrease the borehole length of a single U-tube configuration by 1.3% when compared to pure water. Nanofluid technology is still relatively new for use as a working fluid in GHS systems; hence, more research (both theoretical and empirical) should be considered.

Other studies examined how the fluid velocity within the pipe affected the heat transfer performance of the GHE. Zhou et al. [8] assessed the influence of varying water velocities through various vertical ground heat exchanger pipes, namely, single U-32, double U-25, and double U-32, on heat transfer performance. The results showed that the appropriate water velocity ranges for single U-32 are between 0.4 and 0.6 m/s, 0.4 and 0.5 m/s for double U-25, and 0.3 and 0.4 m/s for double U-32. In addition, Li et al. [43] suggested that the water velocity range for a single U-tube should be between 0.4 m/s and 0.7 m/s for maximum efficiency. Wang et al. [44] demonstrated the impacts of water velocity on the heat transfer efficiency of deep borehole heat exchangers (DBHE) with a depth of 2000 m using a developed numerical model. They examined water velocities ranging from 0.04 m/s to 1.3 m/s and found that a higher heat exchanger rate and a maximum output temperature can be attained when the water velocity is between 0.3 m/s and 0.7 m/s. Another study by You et al. [45] used three different water velocities of 0.26 m/s, 0.51 m/s, and 1.02 m/s to demonstrate the influences of circulating water velocity on the heat exchanger rate, where the water temperature inside the pipe remained constant at 35 $℃$ through the operational test. Their results indicated that the heat exchanger rates were found to be 84 W/m, 116 W/m, and 94 W/m, when the water velocities were 0.26 m/s, 0.51 m/s, and 1.02 m/s, respectively. The maximum heat exchanger rate was achieved when the water velocity was 0.51 m/s; however, the heat exchanger rate decreased when the water velocity was raised to 1.02 m/s. Therefore, the appropriate water velocity was between 0.5 m/s and 0.6 m/s. From this, it can be concluded that the heat exchanger rate per meter increases with water velocity until it exceeds the allowable range of water velocities, at which point the heat exchanger rate decreases. Additionally, You et al. [45] agreed with another study [12] that the optimal water velocity is around 0.5 m/s. Interestingly, Casasso and Sethi [46] demonstrated that optimizing the water velocity can reduce the borehole’s thermal resistance, which can maximize the heat transfer in either direction. According to the findings of these studies, the water velocity should be in an acceptable range, that is neither too high nor too low, to promote better heat transfer and achieve the desired temperature. When the water velocity is too high, more energy is required to circulate the water, but more importantly, the exit temperature from the loop is not different enough from the entrance. When the flow velocity is too low, the heat transfer is decreased, and the water temperature reaches the ground temperature well before the exit, rendering most of the loop useless. Optimally, one would control the velocity such that the exit temperature is the same as the ground and that the entire loop is utilized equally.

Implementing the suggestions mentioned in the GHE’s performance factors section could increase the GHE’s thermal performance to reach the appropriate loop exit temperature. However, further improvements are required to increase its performance. For example, circulating fluid through a U-tube pipe consumes energy, which will increase with higher velocities; hence, controlling the fluid velocity is critical. Therefore, in this paper, we propose a mathematical model of heat transfer behavior between the fluid inside the vertical underground pipes and the surrounding ground in a ground heat exchanger for heating and cooling modes. Model predictive control (MPC) is designed to optimize the GHE’s output temperature by controlling the water velocity.

The rest of the paper is arranged as follows: Section 3 contains the heating and cooling modes. A mathematical model of heat transfer behavior between the fluid in a single vertical underground U-tube pipe and the surrounding ground is presented in Section 4. Section 5 provides a bilinear state–space model, linearization, and discretization. Section 6 presents the MPC design. An experimental GHPS case study is presented in Section 7. Section 8 presents the results and discussion. Section 9 provides a brief conclusion.

3. Heating and Cooling Modes

The heat transfer mechanism between the circulating water inside the vertical U-tube pipe and the surrounding grounds is as follows: through conduction, the heat transfers from the grout to the outer pipe wall and from the outer to the inner pipe wall, and by convection from the inner pipe wall to the circulating water. The GHPS operates in the heating mode when the ground temperature is higher than the ambient air temperature. The heat transfer profile between the fluid inside the vertical underground pipes and the surrounding ground in a ground heat exchanger for the heating mode is presented in Figure 2a. The temperature difference along the pipe is expressed as $\Delta T=T_g-T\left(x\right)$. Here, $T_g$ represents the ground temperature, while $T\left(x\right)$ is the water temperature. Based on the $\Delta T$ value, the heat transfer process can be classified as either high heat transfer, low heat transfer, or no heat transfer. At the beginning of the pipe ($x=0$), the heat transfer flows rapidly from the ground to the water inside the pipe because of the significant (∆T) temperature difference. The water temperature, $T\left(x\right)$, will continue to rise as the water flows through the pipe, and ∆T will drop. The low heat transfer occurs when the water temperature is about to approach the ground temperature, which typically occurs before the pipe’s exit. The water temperature profile, $T\left(x\right)$, barely varies by the end of the pipe ($\Delta T\approx 0$) if the water velocity is not too high and if other system parameters are in good condition for thermal conductivity. Therefore, it can be regarded as reaching the ground temperature. When in the heating mode, the system can benefit from the maximum heat in the ground when the water temperature reaches the equilibrium point. When in the cooling mode, the system works in reverse. Figure 2b demonstrated the heat transfer profile in the GHE for the cooling mode. In this case, the temperature difference $\Delta T$ was calculated as $\Delta T$ $=T\left(x\right)-T_g$. As the water reached ground temperature, the heat in the water was transferred to the ground. Further numerical example explanations are presented in Section 8.

4. The Heat Transfer Model

4.1. Model Assumptions

To study the heat transfer behavior in a vertical single U-tube configuration of GHE, some assumptions are made as follows:

▪

Steady-state.

▪

The heat transfer is one-dimensional, modeled in the x-direction only.

▪

Heat is transferred through conduction between the soil and outside the pipe wall (in cylindrical coordinates); heat is transferred through convection in the fluid inside the pipe.

▪

The area surrounding the outside pipe wall (the borehole wall and grout soils) is treated as a single medium with a single thermal resistance, $R_v$.

▪

Three thermal resistances are considered as follows:

-

$R_{ground}$: resistance of the ground.

-

$R_{wall}$: resistance of pipe wall.

-

$R_{conv}$: resistance of water in the pipe.

▪

The ground temperature (borehole wall and grout soils surrounding the out-pipe wall) is assumed to be at $11.6 ℃$ during both heating and cooling modes (the loop is below the frost line).

▪

Flowing groundwater does not affect heat transfer.

4.2. Governing Equations

A mass flow rate of water flowing in only one direction inside an underground pipe can be calculated using $\dot{m}=\rho vA$, where $\rho$ is the density of the fluid, $v$ is the velocity of the working fluid, and $A$ is the cross-sectional area of the pipe. ${\dot{Q}}_x$ is heat flow in at position x, ${\dot{Q}}_{x+\Delta x}$ is heat flow out at $x+\Delta x$, and ${\dot{Q}}_p$ is the heat coming out of the water to the ground, as shown in Figure 3.

The conservation of energy equation was used to model the heat transfer that can be conducted through the pipe from $x$ to $\Delta x$. The conservation of energy equation is

$$\frac{\partial Es}{\partial t}=\dot{Q}+{\dot{W}}_s+\Sigma \dot{m}{\left(h+\frac{1}{2}v+yz\right)}_i-\Sigma \dot{m}{\left(h+\frac{1}{2}v+yz\right)}_o,$$

(1)

where E is the system energy, $W_s$ is the system work done (all types of work other than flow work), Q is the rate of heat transfer, $\dot{m}$ is the mass flow rate of the water, $v$ is water velocity, and t is the time. The total rate of heat transfer is expressed as

$$\dot{Q}={\dot{Q}}_x-{\dot{Q}}_{x+\Delta x}-{\dot{Q}}_p$$

(2)

Substituting the total rate of heat transfer into Equation (1), implementing the assumptions in Section 4.1, and replacing $i,$ and $o$ with $x$ and $x+\Delta x$, respectively, renders

$$0={\dot{Q}}_x-{\dot{Q}}_{x+\Delta x}-{\dot{Q}}_p+\dot{m}\left(h_x-h_{x+\Delta x}\right).$$

(3)

Conduction in the direction of the flow in Equation (3) can be expressed in terms of temperature by utilizing Fourier’s model; thus,

$${\dot{Q}}_x={\left(-k\frac{\partial T}{\partial x}\right)}_x, {\dot{Q}}_{x+\Delta x}={\left(-k\frac{\partial T}{\partial x}\right)}_{x+\Delta x},$$

(4)

where $k$ is the fluid thermal conductivity. The Newton’s Law of Cooling equation is used to calculate the rate of convective heat transfer from the water to the pipe’s wall [47]. The heat exchanger ${\dot{Q}}_p$, between the water and the ground temperature is calculated as follows:

$${\dot{Q}}_p=U_0{A}_0\left[T_{\left(x\right)}-T_g\right],$$

(5)

where $A_0$ is the pipe’s surface area, which equals $A_0=P_0\Delta x$, $U_0$ is overall heat transfer coefficient, $T_{\left(x\right)}$ is the temperature, $T_g$ is the ground temperature, and $P_0$ is the outside periphery of the pipe. Substituting ${\dot{Q}}_x$, ${\dot{Q}}_{x+\Delta x}$, and ${\dot{Q}}_p$ into Equation (3) results in the following:

$$0={\dot{Q}}_x-{\dot{Q}}_{x+\Delta x}-U_0{P}_0\left[T_{\left(x\right)}-T_g\right]dx+\dot{m}\left(h_x-h_{x+\Delta x}\right).$$

(6)

Integrating Equation (6) from $x$ to $x+\Delta x$ yields

$$0={{\displaystyle \int }}_x^{x+\Delta x}{f}_{q_{\left(x\right)}}dA-{{\displaystyle \int }}_x^{x+\Delta x}{f}_{q_{ \left(x+\Delta x\right)}}dA-{{\displaystyle \int }}_x^{x+\Delta x}{U}_0{P}_0\left[T_{\left(x\right)}-T_g\right]dx+\dot{m}\left(h_x-h_{x+\Delta x}\right)$$

(7)

where

$$f{q}_x={\left(-k\frac{\partial T}{\partial x}\right)}_x,f{q}_{x+\Delta x}={\left(-k\frac{\partial T}{\partial x}\right)}_{x+\Delta x}, \mathrm{and} \left(h_x-h_{x+\Delta x}\right)\approx C_p\left(T_x-T_{x+\Delta x}\right).$$

By substituting the integration limits for these $f{q}_x$, and $f{q}_{x+\Delta x}$, we obtain $f{q}_{x+\Delta x}-f{q}_x$, as provided below:

$$0=\left\{{\left(kA\frac{\partial T}{\partial x}\right)}_{x+\Delta x}-{\left(kA\frac{\partial T}{\partial x}\right)}_x\right\}-{{\displaystyle \int }}_x^{x+\Delta x}{U}_0{P}_0\left(T_{\left(x\right)}-T_g\right)dx-\dot{m}{C}_p\left(T_x-T_{x+\Delta x}\right)$$

(8)

Dividing Equation (8) by $\Delta x$, and then substituting the mass flow rate of water, $\dot{m}=\rho vA$, and the limit is taken as $\Delta x\to 0$. These yield in the case of constant properties and cross-sectional area the following:

$$0=kA\frac{{\partial }^{2}T}{\partial x^2}-U_0{P}_0\left(T-T_g\right)-\rho vA{c}_p\frac{\partial T}{\partial x}$$

(9)

The heat equation is presented in Equation (9), a second-order ordinary differential equation. The first term is heat conduction, the second is the rate of heat transfer to the ground, and the third is fluid convection. Dividing Equation (9) by $A$ gives

$$-k\frac{{\partial }^{2}T}{\partial x^2}=-\frac{U_0{P}_0}{A}\left(T-T_g\right)-\rho v{C}_p\frac{\partial T}{\partial x}$$

(10)

where $P_o=\pi d_{o}L$, and $A=\frac{\pi }{4}{d}_i^2$. By substituting these into Equation (10), we obtain the heat transfer equation as follows:

$$-k\frac{{\partial }^{2}T}{\partial x^2}=-4 U_{0}L\left(\frac{d_{0 }}{d_i^2}\right)\left(T-T_g\right)-\rho v{C}_p\frac{\partial T}{\partial x}$$

(11)

4.3. Thermal Resistance Network

Three thermal resistances in series are considered when transferring heat from/to the ground to/from the fluid inside the pipe, respectively. As discussed in the GHE performance factors section, these resistances have an impact on GHE performance, which is dependent on the thermal conductivity of the soil, pipe wall material, and working fluid. The thermal resistances are exterior convection resistance, conduction thermal resistance, and internal convection resistance, as shown in Figure 4.

By using the thermal resistance method, the total thermal resistance equals the sum of the three thermal resistances since they are in series. Therefore, the total thermal resistance is expressed in Equation (12):

$$R_{TH}=\frac{1}{h_i{A}_i}+\frac{\ln \left(\frac{d_0}{d_i}\right)}{2\pi kL}+R_g$$

(12)

Multiplying Equation (12) by $\left(\frac{A_0}{A_0}\right)$ gives

$$R_{TH}=\left[\frac{1}{h_i}\left(\frac{A_0}{A_i}\right)\frac{1}{A_0}+\frac{A_0\ln \left(\frac{d_0}{d_i}\right)}{2\pi kL{A}_0}+R_g\frac{A_0}{A_0}\right]$$

(13)

Extracting $A_0$ from Equation (13) to obtain the R-value yields

$$R_{TH}{A}_0=\frac{1}{h_i}\left(\frac{A_0}{A_i}\right)+\frac{A_0\ln \left(\frac{d_0}{d_i}\right)}{2\pi kL}+R_g{A}_0$$

(14)

Simplifying in terms of diameter,

$$R_{TH}{A}_0=\frac{1}{h_i}\left(\frac{d_0}{d_i}\right)+\frac{d_0\ln \left(\frac{d_0}{d_i}\right)}{k}+R_{g}2\pi d_{0}L$$

(15)

Heat convection occurs when heat is transferred through the movement of a liquid or gas rather than through a solid object. When temperature differences between the water and pipe’s wall are substantial, heat convection can be enhanced. The convective heat transfer coefficient is $h_i$, and its units are $\frac{\mathrm{W}}{{\mathrm{m}}^2}°℃$. The heat transfer convection coefficient, which is the most critical parameter in the heat transfer process, is affected by the velocity of the fluid and its natural motion, and the difference in temperatures. The Dittus–Boelter Equation was used to compute the convective heat transfer coefficient as presented in Equation (16) [48].

$$N_{U_{d_i}}=0.023 P_r{}^{n}R{e}_{d_i}{}^{0.8}$$

(16)

$N_{U_{d_i}}$ is the Nusselt number that equals $\left(\frac{h_{i}di}{k_w}\right)$, $P_r$ is the Prandtl number that equals $\left(\frac{\mu Cp }{k_w}\right)$, and $R{e}_{d_i}$ is the Reynolds number of circulating fluid that equals $\left(\frac{\rho v di}{\mu }\right)$. $h_i$ is the convection heat transfer coefficient between the water and the pipe wall, $d_i$ is the pipe inner diameter, $v$ is the water velocity, $k_{water}$ is the thermal conductivity of the water (calculated at the film temperature), $\mathsf{\mu }$ is the dynamic viscosity, and $c_p$ is the specific heat. Here, $n$ depends on the weather. In the winter, when GHPS is to be used for heating, the circulating fluid temperature is lower than the ground temperature, and the exponent of heating $n=0.4$. However, in the summer, GHPS is in its cooling mode, the fluid temperature is higher than the ground temperature, and the exponent of cooling $n=0.3$. Substituting $N_{U_{d_i}}$, $P_r$, and $R{e}_{d_i}$ into Equation (16), the convective heat transfer coefficient is obtained:

$$h_i=\left(\frac{k_{water}}{di}\right)\left(0.023\right){\left(\frac{\mu Cp}{k_{water}}\right)}^n{\left(\frac{\rho v di}{\mu }\right)}^{0.8}$$

(17)

The overall heat transfer coefficient, $U_0$, measures how much heat transfers through these three series resistances. This means that when the $U_0$ is large, heat is transferred from the hot region to the cold at a higher rate. The following equation describes how to calculate the overall heat transfer coefficient $U_0$:

$$U_0={\left\{\frac{1}{h_i}\left(\frac{do}{di}\right)+\frac{do ln\left(\frac{do}{di}\right)}{k_p}+R_v\right\}}^{-1}$$

(18)

where $\frac{1}{h_i}\left(\frac{do}{di}\right)$ represents the R-value of water convection inside the pipe, $\frac{do ln\left(\frac{do}{di}\right)}{k_{pipe}}$ represents the pipe wall conduction R-value, and $R_v$ represents the ground resistance R-value. As was assumed earlier in the model assumption in Section 4.1, the ground and grout regions are treated as a single thermal medium with the same singular resistance known as $R_v$, and its units are (m²·°C)/W.

It can be seen from Equation (18) that the overall heat transfer coefficient $U_0$ is inversely proportional to the sum of R-values of the heat transfer components. Because the pipe wall has the lowest thermal conductivity in contrast to soil and water, it can act as a physical barrier that can reduce heat transfer between the ground and fluid inside the pipe, but that is mitigated with a relatively small wall thickness.

As a result, we have the steady-state heat transfer equation in its final form, which is a nonlinear second-order ordinary differential equation:

$$-k\frac{{\partial }^{2}T}{\partial x^2}=-\rho v{c}_p\frac{\partial T}{\partial x}-4\left(\frac{do}{di}\right)U_0\left(T-T_g\right)$$

(19)

5. Bilinear State Space Model

The GHE’s heat transfer process between the water inside the pipe and the surrounding ground balances the water velocity with the exit temperature, where the water velocity is an input control variable, the ground temperature is a disturbance to the system, and the exit temperature is an output of the system. Due to the flow of water inside the pipe at varying velocities, and the heat transfer dynamic between the water inside the pipe and the surrounding ground, the system has nonlinear characteristics. This type of nonlinearity is known as a bilinear system (BLS), and it occurs when the system’s dynamic behavior is affected by a process or an unmeasured disturbance. The behavior of a bilinear system can be found in a variety of applications, including chemical processes, physics, and thermal energy. For instance, in GHE, the heat transfer process is affected by the working fluid velocity, which is an inherited bilinear phenomenon. Furthermore, the bilinear system can be expressed as an input–output system or state–space model system. In general, the continuously bilinear system can be written concisely as follows:

$$\left\{\begin{array}{c}{x}^{\prime }\left(t\right)=A\left(t\right)x\left(t\right)+B\left(t\right)x\left(t\right)u\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Du\left(t\right)\end{array}\right.$$

(20)

where $x\left(t\right)$ is a state vector (i.e., the water temperature), $u\left(t\right)$ is an input control vector (i.e., the GHE system water velocity), and $y\left(t\right)$ is a control output vector (i.e., GHE’s output temperature). The A is a state matrix, B is an input matrix, C is an output matrix, and D is a direct transmission matrix for input (which directly connects between the input and output). The mathematical model of the heat transfer Equation (19) was converted to a state–space model form:

$$\left\{\begin{array}{c}\left[\begin{array}{c}{x}_1^{\prime }\left(t\right)\\ x_2^{\prime }\left(t\right)\end{array}\right]=\left[\begin{array}{cc}0& 1\\ \frac{4}{k}\left(\frac{do}{di}\right)U_0& 0\end{array}\right]\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]+\left[\begin{array}{c}0\\ \frac{\rho c_p}{k}{x}_2\end{array}\right]u\left(t\right)\\ y\left(t\right)=\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]\end{array}\right.$$

(21)

Since the GHE heat transfer model is nonlinear for water velocity (control input), designing a nonlinear controller to achieve optimal performance could be a complicated challenge with less-than-optimal performance. Therefore, by linearizing the system, designing a model predictive control (MPC) is adequate to attain optimal performance. Here, first-order Tylor series expansion was used to linearize the nonlinear model (19) at a desired operating point based on the initial condition values, which ensured it is both controllable and observable around the linearized point. The linearization result was then verified using MATLAB’s function ‘linmode,’ which returned the nonlinear model blocks built in the Simulink platform into MATLAB Command-Line in a linear state–space form as follows:

$$\left\{\begin{array}{c}\left[\begin{array}{c}{x}_1^{\prime }\left(t\right)\\ x_2^{\prime }\left(t\right)\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 69.25& 0\end{array}\right]\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]+\left[\begin{array}{c}0\\ -785.10\end{array}\right]u\left(t\right)\\ y\left(t\right)=\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{c}{x}_1\left(t\right)\\ x_2\left(t\right)\end{array}\right]\end{array}\right.$$

(22)

Then, we discretized the linearized continuous-time state–space model (22) by using Equation (23) as follows with the sample time (T_s = 0.1 s):

$$A_d=e^{A_c{T}_s}, B_d=A_c{}^{-1}\left(e^{A_c{T}_s}-I\right)B_c, C_d=C_c, \mathrm{and} D_d=D_c$$

(23)

The discretization result was then verified using MATLAB’s function ‘c2dm’, which discretizes the continuous linear state–space model so that a discrete model predictive control (MPC) could be applied.

Thus, the discrete-time state–space equation is represented below:

$$\left\{\begin{array}{c}{x}_d\left(k+1\right)=A_{d}x\left(k\right)+B_{d}u\left(k\right)\\ y_d\left(k\right)=C_{d}x\left(k\right)+D_{d}u\left(k\right)\end{array}\right.$$

(24)

6. Model Predictive Control

Model predictive control (MPC) is a modern approach that relies on the plant model and its current state to solve open-loop optimization problems. MPC can predict future control actions to obtain the control model’s sequence, which can be used to enforce the predicted output signal to track the pre-defined reference signal while the error prediction is minimized. Typically, the first control action is only computed; it then moves forward with a new state and control prediction values while the pre-defined cost function is minimized. The MPC design is based on the system model’s discrete state space matrices (A, B, and C) and its critical parameters, the prediction horizon ($N_P$) and control horizon ($N_C$), used for tuning. These parameters, $N_P$ and $N_C$, define the number of sample predictions and samples to control the system’s output, respectively. Furthermore, their values are chosen based on the following formula: $1\le N_C\le N_P$ [49]. According to MPC’s advantages and robust performance, it may be appropriate to optimize the output temperature of the ground heat exchanger (GHE).

6.1. Prediction

The future system model output is predicted based on estimated control input vectors and future state variables through the prediction horizon $N_P$. Thus, the future input control trajectory and future state variables can be expressed in Equations (25) and (26) as follows:

$$u\left(k\right), u\left(k+1\right), \dots , u\left(k+N_c-1\right),$$

(25)

$$x\left(k+1|k\right), x\left(k+2|k\right), \dots , x\left(k+N_p|k\right)$$

(26)

where $u\left(k\right)$ and $x\left(k\right)$ are input control variable and state vector at time $k$, respectively. Implementing the input control and the state vector equations through iteration yields the predicted sequence of the output $Y$, which is shown below.

$$Y=Fx\left(k_i\right)+\varnothing U$$

(27)

where $F,\varnothing$ are matrices used in the prediction Equation (27)

$$F=\left[\begin{array}{c}CA\\ C{A}^2\\ ⋮\\ C{A}^{N_P}\end{array}\right]; \varnothing =\left[\begin{array}{cccc}CB& 0& & 0\\ CAB& CB& & 0\\ C{A}^{2}B& CAB& \cdots & 0\\ ⋮& ⋮& & ⋮\\ C{A}^{N_P-1}B& C{A}^{N_P-2}B& & C{A}^{N_P-N_C}B\end{array}\right]$$

6.2. Optimization

The control design’s objective is to enforce the GHE’s predicted output temperature $\hat{T}\left(k\right)$ to track the pre-defined referenced signal $r\left(k_i\right)$ by determining the best control vector $u$ (i.e., best water velocity) to minimize the prediction error (between the pre-defined referenced and predicted output temperature) through the prediction horizon $N_P$. The pre-defined reference signal is defined below.

$$r\left(k\right)=\left[\begin{array}{c}r\left(k+1\right)\\ r\left(k+2\right)\\ \begin{array}{c}⋮\\ r\left(k+N_p\right)\end{array}\end{array}\right]$$

(28)

By controlling the water velocity $v\left(k\right)$ (input control variable), the future output temperature trajectory $\hat{T} \left(k+1\right)$ can be predicted. while the cost function is minimized. The cost function equation $J\left(k\right)$ is:

$$J\left(k\right)={{\displaystyle \sum }}_{i=1}^{N_p}Q·[{\left(\hat{T}\left(k+i|k\right)-r\left(k+i|k\right)\right]}^2+{{\displaystyle \sum }}_{i=1}^{N_C-1}R·{\left[\Delta v\left(k+i|k\right)\right]}^2+{{\displaystyle \sum }}_{i=1}^{N_P}{Q}_1·{\left[v\left(k+i|k\right)\right]}^2$$

(29)

subject to

$$\left\{\begin{array}{c}{T}_{min }\le \hat{T}\left(k+i|k\right)\le T_{max}\\ v_{min}\le v\left(k+i|k\right)\le v_{max}\\ \Delta v_{min}\le \Delta v\left(k+i|k\right)\le \Delta v_{max}\end{array}\right.$$

(30)

where $\hat{T}\left(k+i|k\right)$ and $r\left(k+i|k\right)$ are the predicted output temperature and the predefined reference temperature at the time $\left(k+i\right)$, respectively, while $v\left(k+i|k\right)$ and $\Delta v\left(k+i|k\right)$ are predicted manipulated control and predicted rate manipulated control changes at the time $\left(k+i\right)$, respectively. The $Q_1$, $Q_2$, and $R$ are the weightings of error output, input, and a variable rate of control change, respectively. Based on the system dynamics prediction, the cost function $J\left(k\right)$ is integrated along the prediction horizon N_P. It is important to note that all these positive weighting coefficients ($Q_1$, $Q_2$, $R$) should be defined in appropriate tuning ranges to ensure that the MPC operates according to the constraints, and then the following control variable sequences are computed to minimize the cost function. At each time step, the state matrices are updated to solve the quadratic problem (QP) to enhance the control performance along with the length prediction. Figure 5 shows the schematic diagram of model predictive control of heat transfer in a vertical ground heat exchanger for a geothermal heat pump system.

7. Case Study of a Geothermal Heat Pump System at Oakland University

7.1. System Description

The 400-ton geothermal heat pump system (GHPS) was installed at Oakland University’s northwest campus in August 2012 to provide space heating and cooling for the Human Health Building (HHB). The HHB is a five-story building with a total area of $159,793{ \mathrm{m}}^2$ that requires a lot of heat and cooling energy; consequently, a vertical ground heat exchanger (GHE) with 256 vertical boreholes is appropriate. The GHPS typically consists of a vertical closed-loop ground heat exchanger (GHE), a heat pump, and a distribution unit system, as shown in Figure 6. The GHE consists of 256 vertical borehole arrays that were drilled to a depth of $98 \mathrm{m}$. A $7.62 \mathrm{m}$ distance between two centered boreholes ensured that there was no interference between two U-tube GHE pipe’s legs. A high-density polyethylene (HDPE) pipe was inserted into each borehole as a single U-tube borehole configuration, with inner and outer diameters of $0.026 \mathrm{m}$ and $0.032 \mathrm{m}$, respectively. These pipes were then connected to the heat pumps in parallel circuits. Backfilling the space between the pipe’s outer wall and the borehole wall with thermally enhanced bentonite grout with a high thermal conductivity created a strong connection between them.

Table 2. HDPE piping system geometric characteristics.

Symbol	Value	Description
$k_p$	$0.5 \mathrm{W}/\mathrm{m} ·℃$	Thermal conductivity
$d_i$	$0.026 \mathrm{m}$	Inner diameter
$d_o$	$0.032 \mathrm{m}$	Outer diameter
${\mathsf{\rho }}_p$	$959.98 \mathrm{kg}/{\mathrm{m}}^3$	Density
$L$	$98 \mathrm{m}$	length

demonstrates the geometric characteristics of the HDPE pipe. The ground loop pipes were filled with water (water-to-air, which is mostly used in the U.S.) for underground circulation that can promote heat exchange between the water inside the pipe and the surrounding ground by continually circulating the water. Forty-four heat pumps were connected parallel to two circulating 60 horsepower (HP) variable-frequency drive (VFD) water pumps that circulate the water inside the pipes at varying velocities as needed; when the first pump is active, the second is in standby mode [50]. When the heat arrived from the heat pump, the distribution unit circulated air through ducts to keep room temperatures consistent.

Since the GHPS started operating in August 2012, it has been supplying efficient heating and cooling to the HHB. Leidel [50] assessed the GHPS’s performance using its input and output temperature datasets between 2013 to 2015. The results showed that the temperature fluctuation was low, meaning that the underground heat exchanger loop pipes were appropriately sized and designed to provide enough heating and cooling for the entire building with excess capacity to supply more.

To keep the GHPS operating at optimum efficiency so that it can continue to provide the HHB with adequate temperatures in the winter and summer without wasting additional energy (used for heat and water circulating pumps), controlling water velocity is necessary. The MPC controller acts as an energy management system that minimizes energy consumption by optimizing the GHE’s output temperature. This is based on the MPC’s objective to control the water velocity in the U-tube pipe so that the GHE’s output temperature can be optimized by bringing the predicted output temperature as close to the pre-defined referenced output temperature as possible.

7.2. Physical Setup

The water velocity was not constant, but it was assumed to be 0.6 $\mathrm{m}/\mathrm{s}$to determine the heat transfer convective coefficient $h_i$. The convective heat transfer coefficients for both modes were calculated based on Equation (17), and they were 2105.70 and 2182.17 $\mathrm{W}/{\mathrm{m}}^2·℃$ for the heating and cooling modes, respectively.

Once the convective heat transfer coefficient was computed, it was substituted in Equation (18) to find the overall heat transfer coefficient $U_0$. In the ground, the pipe had little room horizontally. This was defined in Equation (18) as $R_v.$The R-value was calculated as 0.1 per inch for the backfill soil material, which varies depending on composition, which in turn affects penetration depth. For example, it ranged from 5 to 15 inches, 40 to 75 inches, and 75 to 95 inches for the clay, sand, and good crushed rock, respectively [51]. In this study, the thermally enhanced bentonite grout with high thermal conductivity was used for the GHE at Oakland University. Thus, we assumed the pipe’s surrounding space to be 30 inches; therefore, $R_v=0.0385{ \mathrm{m}}^2·℃/\mathrm{W}$. The GHE parameter values used in the simulation are presented in

Table 3. Parameters of the ground heat exchanger (GHE).

Symbol	Value		Unit	Description
	Winter	Summer
µ	1.726	0.974	$\mathrm{cp}$	Dynamic viscosity
ρ	999.94	997.96	$\mathrm{kg}/{\mathrm{m}}^3$	Water density
${\mathrm{C}}_{\mathrm{p}}$	4.22	4.15	$\mathrm{kJ}/\mathrm{kg} ·℃$	Specific heat
$k_w$	$0.598$	$0.598$	$\mathrm{W}/\mathrm{m} ·℃$	Water thermal conductivity
$v$	0.35–0.45	0.35–0.45	$\mathrm{m}/\mathrm{s}$	Water velocity
$n$	0.4	0.3	-	Heating/cooling exponent
$h_i$	$2105.70$	$2182.17$	$\mathrm{W}/{\mathrm{m}}^2·℃$	Convective heat transfer coefficients
$U_0$	19.08	19.08	$\mathrm{W}/{\mathrm{m}}^2·℃$	Overall heat transfer coefficient
$R_v$	0.0385	0.0385	${\mathrm{m}}^2·℃/\mathrm{W}$	Soil thermal resistance
$T_{in}$	$2$	$22$	$℃$	Input temperature
$T_{\mathrm{g}}$	$\mathit{11.6}$	$\mathit{11.6}$	$℃$	Ground temperature

.

The overall heat transfer coefficient for both heating and cooling modes was 19.08 $\mathrm{W}/{\mathrm{m}}^2·℃$. Despite differences in convective heat transfer coefficients for heating and cooling modes, the overall heat transfer coefficients were the same.

8. Results and Discussion

8.1. Implementation of the Heat Transfer Model in a Vertical Ground Heat Exchanger for Heating and Cooling Modes

The heat transfer behavior between the water temperature inside a vertical single U-tube pipe and its surrounding ground model for heating and cooling modes, and verification of the effect of the water velocity were demonstrated through simulations. To validate this, the proposed model was applied to evaluate the heat transfer behavior in GHE at Oakland University in both modes. The proposed heat transfer model (19) then was written in MATLAB’s script and implemented using an ordinary differential equation (ODE) solver with two initial conditions: $T\left(0\right)=T_{in}$, and $T^{\prime }\left(0\right)=0$. The model’s parameters defined in Table 2 and Table 3 were substituted into the same equation. The initial temperatures were set to $T\left(0\right)=$ $2 ℃,$and $22 ℃$ for winter and summer, respectively, referring to the temperature at $x=0$ and flowing through the length of pipe$L$. The ground temperature $T_g$ remained constant at $11.6 ℃$for both modes. Additionally, the working fluid (water) flowed at a uniform rate through the pipe at different velocities; however, the water velocity did affect the heat transfer performance and the exit temperature. Thus, four different water velocities, namely, $0.35, 0.45,$ 0.9, and $1.2 \mathrm{m}/\mathrm{s}$, respectively, were applied for the heating and cooling modes, thereby demonstrating how the water velocity affects the loop temperature. Consequently, the heat transfer model was implemented, and the water temperature behavior inside a vertical single U-tube pipe for heating and cooling modes was observed (depicted in Figure 7 and Figure 8). As expected, the GHE’s output temperature matched the ground temperature, which validated the proposed model.

The water temperature behavior inside a vertical single U-tube pipe at different velocities for the heating mode is depicted in Figure 7. The temperature T(x) increased exponentially during water flow inside the pipe from $2 ℃$ until it reached the ground temperature of $11.6 ℃$ at $90 \mathrm{m}$ and $95 \mathrm{m}$ from the pipe length when the water velocity was 0.35 and 0.45 $\mathrm{m}/\mathrm{s}$, respectively. The temperature then gradually remained stable throughout the pipe’s exit. In this case, the ground heat was transferred to the water inside the pipe. In Figure 7, we observed that the output temperature $T\left(x\right)$ approached the ground temperature when the water velocity was between 0.35 and 0.45 $\mathrm{m}/\mathrm{s}$. However, it did not approach the ground temperature when the water velocity was increased to 0.9 and 1.2 $\mathrm{m}/\mathrm{s}$, despite reaching the pipe’s exit. Considerably more heat energy was wasted at these velocities.

Figure 8 demonstrates the water temperature behavior inside a vertical single U-tube pipe at different velocities for the cooling mode. According to Figure 8, the output temperature $T\left(x\right)$ decreased exponentially along the pipe from $22 ℃$ until it approached the ground temperature of $11.6 ℃$ at $90.5 \mathrm{m}$ and $95.5 \mathrm{m}$, when the water velocity was 0.35 and 0.45 $\mathrm{m}/\mathrm{s}$, respectively. That is, the water heat was transferred to the ground. However, the temperature $T\left(x\right)$ did not reach the ground temperature when the water velocity was increased from 0.9 to 1.2 $\mathrm{m}/\mathrm{s}$, as shown in Figure 8.

According to this investigation of water velocity, it affected the temperature $T\left(x\right)$, which can impact the GHE’s thermal performance, as shown in Figure 7 and Figure 8. As the water velocity increased, there was less time for the water to reach the ground temperature, resulting in significant temperature differences ($\Delta T$). This led to a larger margin of error because the water was flowing too fast to effectively reach the appropriate level. Thus, the temperature $T\left(x\right)$ only reached that of the ground when the water velocity was lowered to the point where it had time to heat or cool.

The end of the pipe, $98 \mathrm{m}$, was chosen as a representative point in the comparison process. A comparison of the heat transfer in heating and cooling modes was performed, as depicted in

Table 4. Comparison of model heat transfer performances in heating and cooling modes.

Water Velocity (m/s)	$$\begin{gathered} \mathbf{\hspace{1em}Heating} \mathbf{Mode} \\ \mathit{\Delta }\mathit{T}={\mathit{T}}_{\mathit{g}}-\mathit{T}\left(\mathit{x}\right)\left[℃\right], \end{gathered}$$ (the Pipe’s End, 98 m)	$$\begin{gathered} \mathbf{\hspace{1em}Cooling} \mathbf{Mode} \\ \mathit{\Delta }\mathit{T}=\mathit{T}\left(\mathit{x}\right)-{\mathit{T}}_{\mathit{g}}\left[℃\right], \end{gathered}$$ (the Pipe’s End, 98 m)
0.35	0	0
0.45	0	0
0.9	0.86	0.89
1.2	1.58	1.65

.

When the water velocity was 0.33 and 0.43 $\mathrm{m}/\mathrm{s}$, $\Delta T$ was 0 $℃$ for both heating and cooling modes, whereas when water velocity was 0.9 $\mathrm{m}/\mathrm{s}$, $\Delta T$ was 0.86 $℃$ and 0.89 $℃$ in heating and cooling, respectively. In addition, when the water velocity was 1.2 $\mathrm{m}/\mathrm{s}$, $\Delta T$ was 1.58 $℃$ and 1.65 $℃$ in heating and cooling, respectively. As a result, when the GHPS was set to its heating mode, the output temperature of GHE reached the ground temperature slightly faster than when it was set to its cooling mode. This indicates that heat transfer from the ground to the water inside the pipe was slightly faster than heat transfer from the water to the ground.

According to the findings, water velocities of 0.35 and 0.45 $\mathrm{m}/\mathrm{s}$ were ideal for heat transfer. On the other hand, the velocities of 0.9 $\mathrm{m}/\mathrm{s}$ and 1.2 $\mathrm{m}/\mathrm{s}$ were too fast, not giving the water enough time to sufficiently reach the ground temperature. Pushing the water at high velocities provided no benefit, as it neither delivered an appropriate temperature for use nor saved energy; this has an impact on the heat pump’s operating efficiency. However, if the water velocity is kept low, the temperature $T\left(x\right)$ will reach ground temperature before the pipe’s exit and will remain steady throughout the pipe, rendering the remaining piping loop useless. Therefore, the water velocity could be lowered sufficiently to reach the appropriate temperature without wasting any additional resources; for example, the circulation water pump needs more energy to increase its velocity. In addition, the heat pump needs to compensate for the heat that is not being delivered to the building.

To further demonstrate the performance of the proposed heat transfer model, it was used to simulate the GHE’s output temperature behavior for both U-tube De 25 $\mathrm{mm}$ (outer diameter of 25 $\mathrm{mm}$) and U-tube De 40 $\mathrm{mm}$ (outer diameter of 40 $\mathrm{mm}$) pipe diameters, as well as to identify the acceptable range of water velocities for these U-tube pipes. Note that the parameters used in this experimental study were the same as those used in the previous implementation of heat transfer behavior in the GHE at Oakland University, as shown in Table 2 and Table 3. Therefore, the acceptable ranges of water velocity were computed and are shown in

Table 5. Acceptable ranges of water velocity in varying vertical single U-tube diameters.

U-Tube Pipe Diameters	Range of Water Velocity
U-tube De 25 mm	$0.33–0.43 \mathrm{m}/\mathrm{s}$
U-tube De 32 mm	$0.35–0.45 \mathrm{m}/\mathrm{s}$
U-tube De 40 mm	$0.38–0.48 \mathrm{m}/\mathrm{s}$

.

In summary, the acceptable range of water velocity in a vertical single U-tube GHE pipe at Oakland University is between 0.35 and 0.45 $\mathrm{m}/\mathrm{s}$ to ensure that the heat pump delivers the appropriate temperature to provide HHB with comfortable temperatures regardless of the season. According to our simulation results, the reasonable range of water velocity in vertical single U-tube De 25 $\mathrm{mm}$, De 32 $\mathrm{mm}$, and De 40 $\mathrm{mm}$ pipes are between 0.33 and 0.43 $\mathrm{m}/\mathrm{s}$, 0.35 to 0.45 $\mathrm{m}/\mathrm{s}$, and 0.38 to 0.48 $\mathrm{m}/\mathrm{s}$, respectively.

8.2. MPC Implementation and Simulation Results

According to the simulation results (see Figure 7 and Figure 8), circulating fluid in the U-tube pipe affects the GHE’s output temperatures and consumes energy, which will especially increase with higher velocities. To overcome this problem and optimize the GHE’s output temperature, which can improve energy consumption, we designed and implemented an MPC that can control the water velocity to optimize the GHE’s output temperature and keep it at the appropriate temperature. The MPC controller was designed to solve the quadratic problem (QP) subjected to input constraints $-0.5\le u\left(k\right)\le 0.5$ along with the length prediction. The MPC’s objective is to control the water velocity so that it can enforce the GHE’s output temperature to track the pre-defined reference temperature as accurately as possible. This could be even more accurate because there are constraints on both the input and output signals. The MPC parameters used in this study were $N_P=12, N_C=2$, and $T_s=0.01 \mathrm{s}$. All these parameters were accurately chosen and properly tuned to achieve strong performance and to reduce the error prediction between the predicted and pre-defined reference output temperature values over the prediction horizon. Consequently, the MPC model was implemented using the MATLAB script and the Simulink platform, and the GHE’s output temperature response and control signal were observed (depicted in Figure 9 and Figure 10).

.

According to Figure 9, the GHE’s output temperature tracked the reference temperature trajectory precisely, with a slight overshoot. This result reflects the optimal performance of the control design, and its performance was influenced by water velocity fluctuations.

9. Conclusions

The heat transfer model was proposed to evaluate the heat transfer behavior in the ground heat exchanger (GHE) between the water in the vertical U-tube pipe and its surrounding grounds for both the heating and cooling modes. The proposed model was applied to a case study, Oakland University’s GHPS, to evaluate its validation and enhance the GHE’s performance. Then, the model simulation results were performed using an ODE solver in MATLAB (R2021a) software. According to our findings, the output temperature only reached the ground for both modes when the water velocity was between 0.35 $\mathrm{m}/\mathrm{s}$ and 0.45 $\mathrm{m}/\mathrm{s}$. This is, the heat transfer between two different heat sources required sufficient time to complete. The proposed model was validated when the GHE’s output temperature reached that of the ground. Furthermore, the influence of water velocity on the GHE’s output temperature was discussed. We identified that the optimum ranges of water velocity in vertical single U-tube pipes with diameters of De 25 $\mathrm{mm}$, De 32 $\mathrm{mm}$, and De 40 $\mathrm{mm}$ at a depth of $98 \mathrm{m}$ was between 0.33 and 0.43 $\mathrm{m}/\mathrm{s}$, 0.35 to 0.45 $\mathrm{m}/\mathrm{s}$, and 0.38 to 0.48 $\mathrm{m}/\mathrm{s}$, respectively (see Table 5). Interestingly, the validation of the water temperature behavior in the ground heat exchanger will be the subject of future research.

The circulating fluid in the U-tube pipe affects the GHE’s output temperatures and consumes energy, which will increase as velocities increase. Therefore, we designed and implemented a model predictive control (MPC) controller to optimize the GHE’s output temperature by controlling the water velocity. The MPC was able to predict the future GHE’s output temperature, which resulted in well-tracking performance with the reference temperature trajectory, which was proven by the simulation result. This led to an improved reduction in the energy consumption used for heat and circulating water pumps and ensured that the heat pump delivered the appropriate temperature to provide the HHB with comfortable temperatures regardless of the season. We concluded that the designed linear MPC can control nonlinear heat transfer in a vertical ground heat exchanger system around the desired equilibrium point. In the future, we anticipate controlling this nonlinear system at several operational points, controlling each point separately, and switching between these controllers using the single controller block to control it more effectively along its trajectory.