Optimization model analysis of centralized groundwater source heat pump system in heating season

Abstract

The ground-water heat-pump system (GWHP) provides a high efficient way for heating and cooling while consuming a little electrical energy. Due to the lack of scientific guidance for operating control strategy, the coefficient of performance (COP) of the system and units are still very low. In this paper, the running strategy of GWHP was studied. First, the groundwater thermal transfer calculation under slow heat transfixion and transient heat transfixion was established by calculating the heat transfer simulation software Flow Heat and using correction factor. Next, heating parameters were calculated based on the building heat load and the terminal equipment characteristic equation. Then, the energy consumption calculation model for units and pumps were established, based on which the optimization method and constraints were established. Finally, a field test on a GWHP system in Beijing was conducted and the model was applied. The new system operation optimization idea for taking every part of the GWHP into account that put forward in this paper has an important guiding significance to the actual operation of underground water source heat pump.

Keywords

optimization model / groundwater source heat pump system / theoretical analysis / example verification / heating season

Introduction

Building energy consumption accounts for 22% of the total energy consumption in China, of which 60% are consumed by heating and air conditioning systems [ 1]. During the past few years, the renewable energy has regained attention because of fossil energy shortage and pollution. Traditional energy, such as coal and gas, is gradually substituted by renewable energy with emerging technologies. As one of the proven technologies, the groundwater source heat pump (GSHP) system, which absorbs energy from groundwater for building heating and cooling, has been widely applied in China [ 2]. According to Chinese relevant statistics, more than 400 companies and institutes manufacture, design and construct the GSHP system and over 7000 GSHP projects, most of which were centralized systems, had been completed by the end of 2009 [ 3].

Compared with the traditional central air conditioning system, the GSHP system has been proven to operate with higher efficiency and better economy. However, deviations and problems emerge frequently in projects due to the complex control processes of the GSHP system. Besides, the fact that practitioners lack theories of operation strategy also limits the further development of this promising technology. Optimization models based on theoretical analysis and numerical calculation are often used to evaluate and improve the performance of GSHP systems. A GSHP system consists of the groundwater energy source, GSHP units, water pumps, pipeline network and terminal devices. Although each component in the system has its unique properties and different optimum conditions, there exist interactions and mutual promotions between them. Therefore, the optimization model is an integration of several interconnected calculation models.

Heat transfixion occurs in the groundwater aquifer, where the recharge water causes a temperature variation of well water by seepage and conduction processes [ 4]. Heat transfixion is one of the great concerns. To avoid heat transfixion, designers need to determine the distribution of wells and distance between each well by using groundwater calculation model before the construction work starts [ 5]. In addition, for existing projects, the study of groundwater calculation models can provides optimum methods for minimizing the effect of heat transfixion and response measures when heat transfixion occurs. Other components of a GSHP system are above the ground. When building load changes, these components need to be adjusted together by changing the flow of water or the temperature difference between supply and return. In general, the phenomenon of “great temperature difference but small flow” which mainly dominated by GSHP units and water pumps, is favorable for reducing energy consumption and groundwater exploitation [ 6]. The GSHP unit consumption calculation model and water pump consumption calculation model can help the GSHP system to operate in a favorable condition with a lower consumption.

Many studies have documented the GSHP system optimization model. Li et al. [ 7] analyzed the full load intermittent running in low-load operation of GSHP which consumed large amount of extra energy. The running system was optimized by theoretical calculation. The result showed that there was an optimum running load causing the system running in lower energy consumption. Furthermore, the two-dimension space composed of thermal capacity and the ratio between dissipated heat energy and the maximum can be divided into several running regions in which the optimum value was kept invariable. Zeng and Cai [ 8] established the energy consumption model of the main equipment of the primary pump constant flow rate system with the optimal objective function of integrated energy consumption. Meantime, the dynamic programming was computed with the Matlab software to optimize the air-conditioning water system operation strategy. The optimal operation strategy was compared with a real project, and the optimal operation strategy based on the integrated energy consumption was verified to have a better effect. Researchers also conducted studies of groundwater calculation models. Nam and Ooka [ 9] simulated a 3D numerical heat-water transfer to develop the optimization method for GWHP systems and found that the GWHP system could achieve a higher coefficient of performance (COP) and save more energy by utilizing the relatively stable temperature of groundwater. The performance of the system depended on the condition of groundwater, especially temperature, which affected the performance of the heat pump and system. Pan et al. [ 10] built a coupling numerical model of underground water flow and temperature to analyze the extent of underground water heat transfixion. Besides, the change of underground water temperature for different heat transfer temperature differences was simulated. Based on the actual monitoring, the accuracy of simulated results was proved.

An optimization model was built by integrating the groundwater thermal transport calculation model and the GWHP energy consumption calculation models. In addition, a case was introduced to validate the correctness and applicability of the model. Then, results and conclusions were proposed.

Calculation models

Groundwater thermal transfer calculation model

The study of the flow field and temperature field of the groundwater is very significant for improving the efficiency of GWHP. The ground water thermal transport process is very complicated, during which the flow and heat transfer are mutual coupled. This process in groundwater aquifers can be named as the convection-diffusion heat transfer process.

Mathematical model

The flow process of groundwater always satisfies the continuity equation [ 11]. The flow field and the temperature field can be described in a mathematic way by the flow model and thermal transport model, respectively. Then reasonable boundary conditions and definite conditions are selected to solve the two models.

1) Groundwater flow model

Several basic assumptions for groundwater flow are made to simplify the model. Groundwater flow should be a saturation flow; groundwater flow should satisfy Darcy’s law; groundwater aquifer is elastic and compressible; the coordinate system is in the same direction with the aquifer anisotropy main axis; vertical z-axis upward is the positive direction; the viscosity of groundwater does not change with the temperature; and the flow caused by the density gradient can be neglected [ 12].

Darcy’s law [ 13] can be expressed as

(1) v={\displaystyle \frac{Q}{A}}=K{\displaystyle \frac{h_1-h_2}{l}}={\displaystyle \frac{k\rho g}{\mu }}J,

where v is the seepage velocity, which is not the actual velocity of the fluid, m/d; Q is the seepage discharge per unit time, m³/h; A is the area of the cross section, m²; K is seepage coefficient, m/d; J={\displaystyle \frac{h_{\mathrm{1}}-h_{\mathrm{2}}}{l}} is the hydraulic gradient; k is the permeability; ρ is the density of the fluid, kg/m²; g is gravitational acceleration, m/s²; and μ is the dynamic viscosity of the fluid, kg/(m·s).

The relationship between the actual velocity of the fluid and the seepage velocity is expressed as

(2) V={\displaystyle \frac{v}{n}},

where V is the actual velocity of the fluid, m/s; and n is porosity, %.

For an aquifer of uniform thickness, the continuity equation of groundwater can be obtained by analyzing a micro unit dx dy dz.

(3) {\displaystyle \frac{\partial \rho }{\partial t}}+{\displaystyle \frac{\partial \rho v_x}{\partial x}}+{\displaystyle \frac{\partial \rho v_y}{\partial y}}+{\displaystyle \frac{\partial \rho v_z}{\partial z}}=\mathrm{0.}

Furthermore, the seepage equation of groundwater can be established by taking the origin flow into account and writing the seepage velocity as Darcy’s law shows.

(4) {\displaystyle \frac{\partial (n\rho )}{\partial t}}=\nabla \left[{\displaystyle \frac{\rho K}{\mu }}(\nabla P+\rho g\nabla z)\right]+{\rho }_{s}q,

where P is fluid pressure, Pa; t is time, s; ρ_s is the density of the origin flow, kg/m³; q is the volume flux intensity of the origin flow, m³/ (m²⋅s); and \nabla is Hamilton operator.

As a common variable in seepage process, the fluid pressure is usually used for establishing the initial condition and the boundary conditions.

(5) {P|}_{t=\mathrm{0}}=P(x,y,z,\mathrm{0})(x,y,z\in \mathrm{\Gamma }),

(6) {P|}_{\mathrm{\Gamma }}=P(x,y,z,t)(x,y,z\in \mathrm{\Gamma }),

{K_x{\displaystyle \frac{\partial P}{\partial z}}|}_{\mathrm{\Gamma }}=q_x(x,y,z,t)(x,y,z\in \mathrm{\Gamma }),

(7) {K_y{\displaystyle \frac{\partial P}{\partial z}}|}_{\mathrm{\Gamma }}=q_y(x,y,z,t)(x,y,z\in \mathrm{\Gamma }),

{K_z{\displaystyle \frac{\partial P}{\partial z}}|}_{\mathrm{\Gamma }}=q_z(x,y,z,t)(x,y,z\in \mathrm{\Gamma }),

(8) {\displaystyle \frac{\partial P}{\partial n}}+f_{1}P=f_2,

where Г is the area of the calculation; n is the outer normal direction of the boundary; and f₁ and f₂ are known functions.

The groundwater flow model includes the seepage equation of groundwater Eq. (4), the initial condition Eq. (5) and the boundary conditions Eqs. (6)-(8).

2) Thermal transport model

Several basic assumptions for thermal transport are also made to simplify the model. The kinetic energy and the thermal radiation of the fluid are ignored; the thermal diffusion is similar to mass diffusion; the heat transfer happening in the fluid and in the porous media is simultaneous; the coordinate system is in the same direction with the aquifer anisotropy main axis; specific heat and thermal conductivity are constants; liquid and solid matter is always in thermal equilibrium state; and the enthalpy change caused by pressure, the fluid viscosity change and porous media deformation caused by temperature can be neglected [ 14].

Based on the above assumptions, the heat transfer equation, where the temperature is set as the variable, is expressed as

(9) {\displaystyle \frac{\partial [n\rho c_f+(\mathrm{1}-n){\rho }_m{c}_m]T}{\partial t}}=\nabla [(n{\lambda }_f+(\mathrm{1}-n){\lambda }_m)I\nabla T]+\nabla [n{D}_H\nabla T]-\nabla [\rho V{c}_{f}T]+c_{f}q{\rho }_s{T}_s+q_H,

where T is groundwater temperature, °C; T_s is the temperature of the origin flow, °C; ρ_m is the density of the solid skeleton medium, kg/m³; c_f and c_m are the specific heat of groundwater and the solid skeleton, respectively, J/(kg·°C); λ_f and λ_m are the thermal conductivity of fluid and solid skeleton, respectively, W/(m·°C); I is 3D unit vector; D_H is the mechanical dispersion coefficient tensor of the thermal, W/(m·°C); and q_H is the intensity of the heat source, W/m³.

(10) D_{H,ij}=\rho c_f{D}_{s,ij}=\rho c_f\left[({\alpha }_L-{\alpha }_T){\displaystyle \frac{V_i{V}_j}{V}}+{\alpha }_{T}V{\delta }_{ij}\right],

where D_s,ij is the coefficient of mass diffusion, m²/s; α_L and α_T are the longitudinal and transverse thermal dispersion degree, respectively, m.

The thermal diffusivity is introduced for that direct thermal conduction and diffusion can be combined. Because direct thermal conduction and diffusion can be combined, the coefficient of thermal diffusion is introduced. The heat transfer equation can be changed as

(11) {\displaystyle \frac{\partial [n\rho c_f+(\mathrm{1}-n){\rho }_m{c}_m]T}{\partial t}}=\nabla [D\nabla T]-\nabla [\rho V{c}_{f}T]+c_{f}q{\rho }_s{T}_s+q_H,

(12) D_{ij}=n{D}_{H,ij}+[n{\lambda }_f+(\mathrm{1}-n){\lambda }_m]{\delta }_{ij},

where D is the coefficient of thermal diffusion tensor, W/(m·°C); D_ij is the component of coefficient of thermal diffusion tensor; and δ is dirac symbol, δ_ij = 1, i = j; δ_ij = 0, i≠j.

Numerical solution method

It is difficult to solve the equations mentioned above in the 3D space. So the numerical solution is adopted to get the pumping water temperature. The groundwater thermal migration process during the operation of the GWHP system was simulated by the Flow Heat software, and the results showed that the pumping water temperature can be fitted to a liner function of time.

(13) t_{g,1}=A\tau +B,

where t_g,1 is pumping temperature, °C; τ is time step, h; and A and B are fitting coefficients.

The fitting coefficients were affected by recharge water temperature and flow.

(14) A=f_{\mathrm{1}}(G_g,t_{g,2}),

(15) B=f_{\mathrm{2}}(G_g,t_{g,\mathrm{2}}),

where G_g is recharge water flow, m³/h; and t_g,2 is recharge water temperature, °C.

Besides, the transient heat transfixion in engineering occurs in a short time, which is different from the slow thermal diffusion process. The cause of transient heat transfixion is that high permeability layers distribute around the wells or the block is lax [ 15]. As a result, the recharge water flows directly into the pumping region, not leaving enough time for heat transfer. The pumping water temperature is affected by both the transient heat transfixion and slow heat diffusion process. In addition, the former, of which the characteristic is fast response, has a huge influence on the variation of pumping water temperature than the latter. When the transient heat transfixion happens, the groundwater is short-circuited by the recharge water of which the temperature changes with the heat load [ 16]. Thus, the pumping water temperature cannot get the exact solution through numerical simulation and analytical calculation.

The correction coefficient method was adopted. Since transient heat was greatly influenced by the heat load, the correction coefficient could be expressed as a function of recharge water temperature and flow, as expressed in Eq. (16).The value of the correction coefficient can be fitted using the data obtained from the existing project or those from the project of the same kind.

(16) t_{g,1}^{\prime }=(A_{\mathrm{0}}+A_{\mathrm{1}}{G}_g+A_{\mathrm{2}}{t}_{g,\mathrm{2}})[f_{\mathrm{1}}(G_g,t_{g,\mathrm{2}})\tau +f_{\mathrm{2}}(G_g,t_{g,\mathrm{2}})],

where A₀, A₁, and A₂ are transient heat transfixion correction coefficients.

GSHP energy consumption calculation model

Building load and heating parameters

The building heat load includes the retaining structure heat loss and the infiltrating loss under the condition of ignoring the heat produced by indoor person and equipment. The equation of building heat load can be expressed in Eqs. (17) to (19). For a particular building, there is a linear relationship between the building heat load and the temperature difference.

(17) Q=q_{\mathrm{1}}+q_{\mathrm{2}}=(aKF+\mathrm{0.278}{L}_a\rho c)\times (t_n-t_w)=C(t_n-t_w),

(18) q_{\mathrm{1}}=aKF(t_n-t_w),

(19) q_{\mathrm{2}}=\mathrm{0.278}{L}_a\rho c(t_n-t_w),

L_a=k{l}_{a}l,

where q₁ is retaining structure heat loss, W; q₂ is infiltrating loss, W; a is correction coefficient; K is the heat transfer coefficient of palisade structure, W/(m²·°C); F is the heat transfer area of palisade structure, m²; t_n is indoor design temperature, 18°C; t_w is outdoor temperature, °C; ρ is outdoor air density, kg/m³; c is the specific heat of air, kJ/(kg·°C); l_a is the amount of air infiltrating through per meter aperture, m³/m; l is the crack length calculation for doors and windows, m; and k is direction coefficient.

Various calculation software including DOE-2, ESP-r, DeST, Energy Plus, etc., which use the reaction coefficient method, finite difference method, state space method or other methods, can get the dynamic building heat load by entering building location, meteorological data, basic information of the retaining structures and internal usage (person, lighting and equipment) information. The equation of dynamic building heat load and temperature difference is obtained by linear fitting.

The heating parameters of the GWHP system are mainly affected by the heat supply quantity and the terminal equipment. A balance exists between the building heat load and the heat release of the terminal equipment in an ideal condition.

(20) Q_{\mathrm{1}}=Q_{\mathrm{2}},

where Q₁ is design building heat load, Q₂ is the heat release of terminal equipment at design outdoor temperature. For different types of terminal equipment, Q₂ is a function of the temperature difference between the heat medium temperature and indoor temperature.

(21) Q_{\mathrm{2}}=f(tm-tr),

where t_m is heat medium temperature, °C; t_r is indoor temperature, °C.

For centralized GWHP systems, heat loss exists in the pipe network. Assuming that the heat loss of the supply pipe and the return pipe is equal and the thermal efficiency of pipe network is constant at different water temperature, the calculation equation of actual heating quantity of pipe network is [ 17]

(22) Q_{\mathrm{sys}}=Q_{\mathrm{1}}{(\mathrm{1}-{\eta }_{\mathrm{DL}})}^{-\mathrm{1}},

where {\eta }_{\mathrm{DL}} is distribution loss rate, %.

In summary, when the design building heat load has been calculated and heating parameters have been established, the heating quantity and the supply and return temperature of systems using different terminal types can be obtained by solving Eqs. (20)-(22).

In fact, the real building heat load is always lower than the design value which is calculated in the most unfavorable condition. Thus, the conversion between Q₁ and Q₂ is a ratio rather than the original equal value. The real building heat loads at different outdoor temperatures need to be optimized based on the calculation value. The ratio can be fitted using the observation data of the existing project or of the project in the same mode.

GSHP unit calculation model

The steady and concentrating method is used in establishing calculation models, which can simplify the calculation process and enhance the stability of the system. Steady means that the influence of the time variation can be ignored, while concentrating means that the inner complex and nonlinear characteristics of the GWHP unit are ignored. The model only takes into consideration the import and export fluid parameters.

1)Single unit calculation model

In actual projects, external parameters, including import and export water temperature and flow rate of the condenser and the evaporator, are changing dynamically during operation. However, the water temperature and flow rate of supply water often deviate from the rated condition and cause the performance of the unit to change. Therefore, the performance of the unit under full load condition must be analyzed in order to obtain the unit energy consumption in variable condition.

Keeping the compressor refrigerating capacity constant, there are mathematical relationships between the evaporating or condensing temperature and the refrigerating capacity, as well as the input power. The evaporator heat transfer is a function of evaporation temperature and evaporator import water temperature, while the condenser heat transfer is a function of condensation temperature and condenser export water temperature.

(23) Q_{\mathrm{cp}}=f_Q(t_{\mathrm{0}},t_k),

(24) P_{\mathrm{cp}}=f_P(t_{\mathrm{0}},t_k),

(25) Q_0=F_R{K}_0{A}_{\mathrm{0}}(t_{\mathrm{evp},\mathrm{1}}-t_{\mathrm{0}}),

(26) Q_K=F_R{K}_C{A}_C(t_k-t_{\mathrm{con},\mathrm{1}}),

where Q_cp is compressor heating capacity, kW; t₀ is evaporating temperature, °C; t_k is condensing temperature, °C; P_cp is compressor input power, kW; Q₀ and Q_K are the quantity of evaporator and condenser heat transfer, respectively, kW; F_R is the coefficient related to heat transfer coefficient and flow rate, which is a constant in a certain load range; K₀ and K_C are the heat transfer coefficient of the evaporator and condenser, respectively, W/(m²·K); t_evp,1, and t_con,1 are the import water temperature of the evaporator and condenser, respectively, °C; and A₀ and A_C are the heat transfer area of the evaporator and condenser, respectively, m².

The inner parameters (evaporating temperature and condensing temperature) can be eliminated. The mathematical equations of Q (GWHP unit heating capacity) and P (input power) can be expressed with the outer parameters (import water temperature of evaporator and condenser).

(27) Q=f(t_{con,1},t_{evp,1}),

(28) P=f(t_{con,1},t_{evp,1}).

The quartic polynomial of temperature is recommended to express P by the ASHRAE Handbook.

(29) P=P_{\mathrm{rated}}{\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{2}}{\displaystyle \sum _{j=\mathrm{0}}^{\mathrm{2}}{A}_{ij}{\left(t_{\mathrm{con},\mathrm{1}}-\overline{t_{\mathrm{con},\mathrm{1}}}\right)}^i{\left(t_{evp,1}-\overline{t_{\mathrm{evp},\mathrm{1}}}\right)}^j}},

where P_rated is rated input power, kW; \overline{t_{con,1}} is the regression mean of condenser import water temperature, °C; \overline{t_{evp,1}} is the regression mean of evaporator import water temperature, °C; and A_ij is regression coefficient, determined by the unit actual performance.

The above expression of P can be fitted using the observation data of the existing project or the performance chart of GWHP unit sample.

2) Multi unit calculation model

Most of the time, the GWHP unit operates under partial load condition in actual projects for the reason that building heat load is calculated under the most unfavorable condition and the units are selected with a certain degree of surplus.

(30) P_{\mathrm{pl}}=P·f(\mathrm{PLR})

where P is the input power under full load condition, kW; P_pl is the input power under partial load condition, kW; f (PLR) is correction function; and PLR is the partial load rate of the GWHP unit, %.

For the fact that the heating capacity of the control strategy of the GWHP unit is usually to add or reduce the operated number of the compressor, which is nonlinear, the PLR can be expressed in a form of quadratic polynomial.

(31) f(\mathrm{PLR})=a+b\mathrm{PLR}+c{\mathrm{PLR}}^{\mathrm{2}},

where a, b, c are the correction coefficients of partial load, which can be fitted based on the performance of the GWHP unit under partial load condition.

The survey found that there was usually more than one unit operating in parallel in a GWHP system, which makes the capacity of the GWHP units very huge. Therefore, it is also necessary to distribute the load rate of units to achieve the minimum energy consumption. There are two kinds of strategy for the joint operation of multiple units: working one by one and undertaking load rate evenly. Assume that there is more than one unit of the same type operating in the GWHP system, the distribution of the load rate can be solved by seeking the extreme value of the mathematical model

(32) \{\begin{array}{l}{\displaystyle \sum _{i=1}^N{P}_{pl,i}}={\displaystyle \sum _{i=1}^{N}P·f({\mathrm{PLR}}_i)}={\displaystyle \sum _{i=1}^{N}P(a{\mathrm{PLR}}_i^2+b{\mathrm{PLR}}_i+c),}\\ {\displaystyle \sum Q·{\mathrm{PLR}}_i}=Q_{\mathrm{sys}}\phi ,\\ \mathrm{0}\le {\mathrm{PLR}}_i\le \mathrm{1},\end{array}

where ϕ is the total load rate of the system.

The extreme value

{\mathrm{PLR}}_i={\displaystyle \frac{Q_{\mathrm{sys}}\phi }{NQ}}(Q_{\mathrm{sys}}\phi ＞NQ),

or

{\mathrm{PLR}}_{\mathrm{1}}=\cdots ={\mathrm{PLR}}_M={\displaystyle \frac{Q_{\mathrm{sys}}\phi }{MQ},}

{\mathrm{PLR}}_{M+\mathrm{1}}=\cdots ={\mathrm{PLR}}_N=\mathrm{0}(MQ＜Q_{\mathrm{sys}}\phi \le NQ) (M is integer less than N).

Water pumps calculation model

The groundwater circulation system and supply water circulation system are completely independent and have their own water pumps. Frequency conversion adjustment is one of the most energy saving and economic adjustments. The similarity law concerning the parameters of the water pumps is expressed as

(33) {\displaystyle \frac{n_{\mathrm{1}}}{n_{\mathrm{2}}}}={\displaystyle \frac{G_{\mathrm{1}}}{G_{\mathrm{2}}}}={\left({\displaystyle \frac{H_{\mathrm{1}}}{H_{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{2}}={\left({\displaystyle \frac{P_{\mathrm{1}}}{P_{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{3}},

where n is the rotating speed of pumps, r/min; G is pump flow, m³/h; H is pump head, m; and P is pump input power, W.

Assume that the water pumps are of the same type and the network characteristics remain constant when the operation strategy of the pump changes. The working point of the pump is shown in Fig. 1.

Curve 1 is the H-G curve of the pump under the rated condition. According to the sample performance curve provided by the manufacturer, quadratic equations can be fitted to describe curve 1,

(34) H=a_{\mathrm{0}}+a_{1}G+a_2{G}^{\mathrm{2}},

(35) \eta =b_{\mathrm{0}}+b_{1}G+b_2{G}^2,

where η is the efficiency of water pump, %; and a₀, a₁, a₂, b₀, b₁ and b₂ are the least-square fitting coefficients.

Curve 2 is the H-G curve of the pump during frequency conversion adjustment, whose quadratic equation can be established based on Eqs. (33) and (34).

(36) H=c_{\mathrm{0}}{n}^2+c_{1}nG+c_2{G}^2,

where c₀, c₁, c₂, d₀, d₁ and d₂ are the coefficients corresponding to rotating speed.

Curve 4 is the H-G curve of the pumps in parallel. When N pumps work in parallel, the quadratic equation of curve 4 can be expressed as

(37) H=c_{\mathrm{0}}{n}^{\mathrm{2}}+c_{\mathrm{1}}n\left({\displaystyle \frac{G}{N}}\right)+c_{\mathrm{2}}{\left({\displaystyle \frac{G}{N}}\right)}^{\mathrm{2}}.

Curve 3 is the characteristic curve of the pipe network, whose expression is expressed as

(38) H=H_{\mathrm{0}}+S_1{G}^2,

where H₀ is pump suction height, m, for closed type water pump, H₀ = 0 m; and S₁ is the impedance of the pipe network, s²/m⁵.

According to the pump similarity law, curve 5, a quadratic parabola passing through the origin, is the curve of similar working conditions of the pump. The expression of curve 5 is shown as

(39) H=S_2{G}^2,

where S₂ is the equivalent impedance of the pipe network, s²/m⁵.

In frequency conversion, point A is the working point of the pumps in parallel while point B is the working point of a single pump. The flow rate and head of the pumps at point A can be obtained by combining Eqs. (38) and (39).The head at point B is the same as that at point A. Thus the flow rate and efficiency of the pump at point B can be obtained from the combination of Eqs. (36) and (37). The energy consumption of a single water pump can be calculated using

(40) P={\beta }_{\mathrm{1}}{\beta }_{\mathrm{2}}{\displaystyle \frac{G(H_{\mathrm{0}}+S{G}^{\mathrm{2}})}{\eta }},

where P is the input power of a single pump, kW; β₁ is the flow reserve coefficient, for a single pump β₁ = 1.1, for two pumps in parallel β₁ = 1.2 and β₂ is head reserve coefficient, β₂ = 1.1- 1.2.

Optimization method and constraints

The optimization model method is a kind of method which seeks for the extreme. During the process of optimization, the fact that the GWHP system is a unified organic entirety must be fully realized. Groundwater, building heat load and GWHP equipment are interacted and restricted.

The total energy consumption of the GWHP system is driven by the units, the user side circulating water pumps and well pumps, etc. The terminal equipment energy consumption, which is unstable and low, is not taken into consideration. The optimization process is to solve the minimum value of the objective Eq. (41) at time τ,

(41) E_t=E_{\mathrm{hp},t}+E_{\mathrm{cp},t}+E_{\mathrm{wp},t},

where E_t is the total energy consumption of the GWHP system at time τ, when the outdoor temperature is t, kWh; E_hpt is the total energy consumption of the units at time τ, when the outdoor temperature is t, kWh; E_cp,t is the total energy consumption of user side water pumps at time τ, when the outdoor temperature is t, kWh; and E_wp,t is the total energy consumption of well pumps at time τ, when the outdoor temperature is t, kWh.

In the above optimization model, the variables that can be controlled include PLR, user-side water temperature t_con,2, the frequency of the user side pumps and well pumps. The outer disturbance factors include outdoor temperature t_w and time τ, etc. The building heat load and the heating parameters change with the outdoor temperature. Besides, the groundwater temperature is not only a function of the building heat load, but also a function of the running time. Thus, the optimization process must last throughout the heating period.

The flowchart of the optimization process of the GWHP system is illustrated in Fig. 2.

The establishment of the constraints which is based on the actual situation ensures the operability and safety of the GWHP system. In actual operation, the flow rate must be limited within the prescribed range. So the specific speed of the variable frequency pumps is limited by the inequality constraint as shown below.

(42) \mathrm{0.5}\le V_{\mathrm{p}\mathrm{u}\mathrm{m}\mathrm{p}}\le \mathrm{1.0.}

To ensure the stable and efficient operation of the units, there is a suitable heating capacity adjustment range. The minimum part load rate is determined as 25% for that the different types of units have different ability of capacity adjustment. This optimization model must satisfy the inequality constraints as shown in Eq. (43).

(43)

Therefore, the optimization process of the GWHP system ultimately comes down to solving the nonlinear equations above with the equality constraints and inequality constraints. The optimization variables include the PLR, the system supply water temperature, the groundwater flow rate and the system circulating water flow rate.

Validation and discussion

Overview of project

The project is located in Changping District in Beijing. The building, whose gross area is 40883 m², is a commercial building for car exhibition and spot sale. The GWHP system consists of three central screw GWHP units and fan-coil unit terminals. The designed supply and return water temperature is 45°C /40°C. The type and parameters of the units, pumps and terminals are tabulated in Table 1, Table 2 and Table 3, respectively. The exploration and pumping test indicate that the groundwater buried depth is 22 m. The yield of a single well is 80m³/h and the recharge capacity is 40 m³/h. Six wells support the GWHP system, of which Well 1 and 2 are used as the pumping wells for cooling and heating, respectively. Well 3, 5 and 6 are used as the recharge wells. Well 4 is a pumping and recharge well (PRW), which is used as the subsidiary pumping well when the building load is large. Figure 3 shows the distribution and the distance between each well. Moreover, the distance between Well 1 and 4 is 70 m. Well 2 is 50 m away from Well 4.

The actual measurements of the GWHP system are listed in Table 4.

The total energy consumption of the GWHP system (excluding terminal equipment) during heating season is depicted in Fig. 4. There is a great fluctuation between January 29 and February 13.The reason for this is that the Chinese Spring Festival in 2013 lasted from February 7 to 13 in which the operation time of GWHP system was short and the daily energy consumption very low. The daily COP of the system and the hourly COP of the unit are displayed in Figs. 5 and 6, respectively. In addition to the Chinese Spring Festival and the end of the heating season, the daily energy consumption is more than 2500kW. The COP of the system and of the unit both fluctuate greatly. The average COP of the system and unit is 2.59 and 4.27, respectively, which means that the energy efficiency is comparatively low. The reason for this is that the system lacks operation strategies. Although the outdoor air temperature changes in winter, the pump always adopts a constant flow operation, which leads to a great transport energy consumption.

Figure 7 shows the variations of supply and the return water temperature of the unit system and the variations of pumping and recharge water temperature in a typical day. The fact that groundwater temperature increases approximately 3°C at noon indicates a greater impact of transient heat transfixion. In addition, the temperature difference between supply and return water in the system is only 2.5°C-3.5°C, which is not conducive to the efficient operation of the GWSH system.

Calculation of groundwater

The simulation region, with a scale of 400 m × 400 m × 90 m, was meshed in Flow Heat, as shown in Fig. 8. Six wells were established in the simulation region in accordance with the actual distribution. The X, Y and Z directions were divided into 61, 63 and 6 layers, respectively, and the grids were dense in the vicinity of the wells.

The flow rate was set based on the measured data. During the test, the flow rate of each well was monitored for a week. The results showed that in order to ensure the reinjection, the flow rate remained stable most of the time. The flow rate of Well 1 and 4 was 80 m³/h and 40 m³/h, respectively. The flow rates of remaining four wells were 30 m³/h. The intermittent operation mode was used in this commercial building, so the simulation time of the wells was set from 8:30 am to 16:30 pm, accordingly. The groundwater temperature was obtained as presented in Fig. 9.

Figure 9 indicates that the variation of groundwater temperature is very small (less than 1°C) in the intermittent operation mode for the reason that the recharge water has finished the heat exchanger process with the skeleton of underground aquifer at night completely. Thus, the pumping temperature of groundwater in Eq. (13) can be considered as a constant which is only affected by the transient heat transfixion.

The transient heat transfixion correction coefficients are determined based on the analysis of the observed data. In this project, Well 4 is a PRW and mainly used as a subsidiary one, which is operated from about 8:40 am to 12:30 pm and determined by the building heat load. The measured pumping temperature of two successive typical days is demonstrated in Fig. 10.

The transient heat transfixion correction coefficients were fitted according to the measured data showed in Fig. 10, thus further ground water temperature calculation model was obtained.

(44) t_{g,1}=(\mathrm{0.952578}-\mathrm{0.00934}{G}_g+\mathrm{0.0181}{t}_{g,\mathrm{2}})t_{\mathrm{0}},

where t₀ is the groundwater temperature without the effect of transient heat transfixion, °C.

When the heat load of pumping water (Q_g) remains a constant, there is a mathematical relationship between the pumping and recharge water temperature.

(45) t_{g,\mathrm{1}}-t_{g,\mathrm{2}}={\displaystyle \frac{Q_g}{\rho C{G}_g}},

(46) t_{g,\mathrm{1}}=t_{\mathrm{0}}\times (A_{\mathrm{0}}+A_{\mathrm{1}}{G}_g+A_{\mathrm{2}}{t}_{g,\mathrm{2}}).

Therefore, the calculation equation for ground water temperature can be obtained.

(47) t_{g,\mathrm{1}}={\displaystyle \frac{A_{\mathrm{0}}+A_{\mathrm{1}}{G}_g-A_{\mathrm{2}}{Q}_g/\rho C{G}_g}{\mathrm{1}/t_0-A_{\mathrm{2}}}}.

Calculation of GSHP energy consumption

Building load and heating parameters

The simulation using the eQUEST software indicates that the building heat load is between 20.35 W/m² and 225.75 W/m². The heating season in Beijing lasts from November 15 to March 15 of the next year. The frequency distribution of the daily average outdoor temperature during the heating season in a typical meteorological year is shown in Fig. 11. Based on the building heat load and the daily average outdoor temperature, the relationship between the building heat load and the daily average outdoor temperature can be obtained, as shown in Fig. 12.

There is a significant linear relationship between the building heat load and the daily average outdoor temperature. By using the least square method, the equation between them can be expressed as

(48) Q=\mathrm{151.08}(t_n-t_w)-\mathrm{677.61.}

When the outdoor temperature is the design value, the building heat load is 5920.54 kW, which is calculated with the measured data, while the total heat dissipation of the fan coil units is 7637.5 kW as shown in Table 3. Therefore, the ratio of the total heat dissipation of the terminal equipment and the actual building heat load is 1.29 in design condition. The heat dissipation of the fan coil units is affected by the supply and return water temperature, flow rate, wind speed and indoor temperature [ 18, 19]. The equation for calculating the heat dissipation of the fan coil units in varying conditions is expressed as

(49) Q_{\mathrm{2}}=A_{\mathrm{2}}{a}_{\mathrm{2}}{V}_{\mathrm{y}}^{\mathrm{m}}{V}_{\mathrm{w}}^{\mathrm{n}}(t_m-t_r),

where A₂ is the heat transfer area of the fan coil units, m²; a₂ is the characteristic parameter of the fan coil units; V_y is the heat exchange wind speed of the fan coil units, m/s; and V_w is the flow rate in the fan coil units, m/s.

The t_m at different outdoor temperatures can be established based on Eqs. (48)–(49) and the ratio. The results are presented in Fig. 13.

Figure 13 shows that t_m decreases when the outdoor temperature increases. The average temperature of both supply and return water and heating supply of the system significantly declines because of the oversized area of the terminals and indoor heat gaining. The average temperature of supply and return water drops from the design value of 42.5°C to 40.5°C at the design outdoor temperature (-7°C).

Consumption GSHP units

When the load rate of the units is 100%, the selected measured data of the import and export temperature of the evaporator and condenser, heat capacity and input power are listed in Table 5.

The value of A_ij in Eq. (29) is fitted and shown in Table 6 by using Matlab software to fit the data in Table 5 with the least square method.

When the working environment (inlet water temperature of the evaporator and condenser, etc.) remains constant, the input power of the unit only changes with the load rate. By calculating the ratio of input power under different load rates and under the full load condition, the calculation equation of PLR can be fitted according to Eq. (31).

(50) f(\mathrm{PLR})=\mathrm{0.2326}+\mathrm{0.1015}\mathrm{PLR}+\mathrm{0.8484}{\mathrm{PLR}}^{\mathrm{2}}.

Therefore, the GWHP unit consumption calculation model under variable conditions can be obtained. Besides, there are three identical units in parallel operation in this project and the distribution of the load rate is tabulated in Table 7.

Consumption pumps

There are 4 fixed frequency pumps in parallel operation in the system and the flow rate in the pipe network can be adjusted by operating or switching off the pumps. Three frequency conversion submersible pumps are arranged in Well 1, 2 and 4, respectively. The pump in Well 4 is operating annually. The pump in Well 1 operates mainly in summer and the pump in Well 2 operates in winter. The H-Q and η-Q curve of pumps under the rated condition provided by the manufacturer are shown in Figs. 14 and 15.

The performance equation of each kind of pump can be fitted on the curve by the least square method. The R² are more than 0.96.

Circulating pump:

(51) H=-\mathrm{0.0002}{G}^2+\mathrm{0.}0267G+\mathrm{33.964},

(52) \eta =-\mathrm{0.}0015G^2+\mathrm{0.6707}G+\mathrm{5.0047.}

Submersible pump:

(53) H=-\mathrm{0.00008}{G}^2+\mathrm{0.0164}G+27.989,

(54) \eta =-\mathrm{0.}0006G^2+\mathrm{0.352}G+\mathrm{34.922.}

Based on Eqs. (51) to (54) and the model calculation method, the input power of the pump under the frequency conversion condition can be obtained.

Optimization method

The objective of the optimization model is to get the lowest energy consumption. The optimization variables contain supply water temperature, load rate and frequency of each pump. According to the optimization method and constraints discussed in Subsection 2.3, the GWHP consumption calculation model based on pumps, GWHP units, groundwater temperature calculation model and building heat load forecasting equation is established.

(55) \{\begin{array}{c}{P}_{\tau }=P_{\mathrm{hp,}\tau }+P_{\mathrm{cp,}\tau }+P_{\mathrm{wp,}\tau },\hfill \\ Q_{\mathrm{b,}\tau }=P_{\mathrm{hp,}\tau }+Q_{\mathrm{g,}\tau },\hfill \\ Q_{\mathrm{g,}\tau }=C\rho G_{\mathrm{g,}\tau }\left(t_{\mathrm{g1,}\tau }-t_{\mathrm{g2,}\tau }\right),\hfill \\ Q_{\mathrm{b,}\tau }=C\rho G_{\mathrm{b,}\tau }\left(t_{\mathrm{b1,}\tau }-t_{\mathrm{b2,}\tau }\right),\hfill \end{array}

where P_τ is the total input power of system at time τ, kW; P_hp,τ is the total input power of GWHP units at time τ, kW; P_cp,τ is the total input power of circulating pumps at time τ, kW; P_wp,τ is the total input power of submersible pumps at time τ, kW; Q_b,τ is the building heat load at time τ, kW; Q_g,τ is the groundwater heat load at time τ, kW; G_g,τ and G_b,τ are the flow rate of pumping water and circulating water at time τ, respectively, m³/h; t_g1,τ and t_g2,τ are the pumping and recharge water temperature at time τ, °C; and t_b1,τ and t_b2,τ are the supply and return water temperature at time τ, °C.

Operation strategies should be different when building load is different, and the optimization of consumption calculation model is aimed at the load in each hour. The simulation result of the eQUEST software has a linear relationship with the outdoor temperature which changes each hour; therefore, it can be used for the consumption calculation model.

Optimization results

The building heat load is determined by the outdoor temperature. There is a critical value of outdoor temperature below which the number of operating units is reduced. At this time, the energy consumption of the GWHP units dropped abruptly. The input power of the GWHP unit when the circulating water flow rate is 400m³/h is shown in Fig. 16.

When the outdoor temperature is constant, the flow rate increases, the temperature difference decreases and the temperature of return water rises. So the condensation temperature rises which causes a decline in the energy efficiency ratio. As a result, the energy consumption increases. When the outdoor temperature is 1°C, the input power of the GWHP units changes with the variation of the flow rate, as shown in Fig. 17.

Figure 18 indicates that when the flow rate is lower, the consumption of circulating pumps is lower. However, the consumption of the GWHP units increases.

According to Eq. (47), the underground water temperature can be expressed as a single value function of groundwater flow rate. Meanwhile, with the increase of flow rate, the underground water temperature decreases first and then increases. Therefore, there is a critical value in the flow rate. When the flow rate is less than the critical one, the groundwater temperature decreases, and vice versa. This critical value can be obtained by taking a derivative with respect to the flow rate in Eq. (47).

(56) G_g=\sqrt{{\displaystyle \frac{A_{\mathrm{2}}{Q}_g}{-A_{\mathrm{1}}\rho C}}}.

Based on the actual parameters of the project, the groundwater thermal transfer calculation model, the GSHP energy consumption calculation model, and Eq. (55), the operation modes and parameters of the GWHP system at different outdoor temperatures are determined as shown in Table 8.

According to the measured outdoor temperature during the test, the outdoor temperature hour frequency is added up, as shown in Fig. 19.

Based on the data in Fig. 19 and Table 8, the total energy consumption during the heating season can be worked out, which is 604411.2 kWh, a 7.4% energy saving than before.

Conclusions

In this paper, a thorough research on the operation strategy of a GSHP system was conducted and the optimization model of the GSHP system was obtained, which had a guiding significance on the scientific operation of the GSHP system. The following conclusions can be reached from the analysis in the paper.

1) The groundwater thermal transfer model was built and the temperature varying trend of groundwater in the thermal diffusion process was found by numerical simulation. In addition, the effect of transient heat transfixion was analyzed and correction coefficient based on the measured data of transient heat transfixion was proposed. Finally, the groundwater temperature calculation model at different building heat loads was established.

2) Based on the actual operation data, on the one hand, the GWHP consumption calculation model under variable condition was established; on the other hand, the solution to the water pump consumption under the frequency conversion condition was analyzed.

3) The optimization model of the GWHP system was proposed consideration each component in the GWHP system.

4) By combining with case analysis and comparing with the measured value, an optimization result of 7.4% of energy saving than before was obtained, validating the optimization model of the GSHP system.

References

[1] Diao N R, Cui P, Liu J H, Fang Z H. R&D of the ground-coupled heat pump technology in China. Frontiers of Energy and Power Engineering in China, 2010, 4(1): 47–54

[2] Zhang P M. The comprehensive evaluation of techno-economical and environmental benefit of ground source heat pump system. Dissertation for the Master’s Degree. Beijing: Tsinghua University, 2011 (in Chinese)

[3] Xu W. The development of ground source heat pump technology in domestic and international. Construction Science and Technology, 2010, (18): 14–18 (in Chinese)

[4] Ma J C, Zhao J, Li P L. The Researches of the couple numerical model of groundwater flow, heat transferring and solute movement and the application of energy storage in brackish aquifers. Advances in Intelligent and Soft Computing, 2012, 112: 573–586

[5] Xue Y, Xie C, Li Q. Aquifer thermal energy storage: a numerical simulation of field experiments in China. Water Resources Research, 1990, 26(10): 2365–2375

[6] Lu S L, Cai Z, Zhang L, Li Y R. Evaluation of the performance of a centralized ground-water heat pump system in cold climate region. Frontiers in Energy, 2014, 8(3): 394–402

[7] Li B X, Wang Y B, Jiang B C, Obara S. Optimum operating of ground source heat pump in low-load. Taiyang Neng Xuebao, 2003, 24(5): 713–717 (in Chinese)

[8] Zeng Q X, Cai L J. Optimal analysis on air-conditioning water system operation strategy based on the integrated energy consumption. Heating Ventilation & Equipment, 2010, 38(29): 27–35 (in Chinese)

[9] Nam Y J, Ooka R. Numerical simulation of ground heat and water transfer for groundwater heat pump system based on real-scale experiment. Energy and Building, 2010, 42(1): 69–75

[10] Pan J, Zang H Y, Gao W C, Zuo Y. Application of numerical simulation of underground water temperature to groundwater source heat pump engineering.  HV&AC,  2010,  40(10):  67–70  (in Chinese)

[11] Yuan J W. Energy abstraction on geo-temperature field of groundwater heat pump. Dissertation for the Master’s Degree. Beijing: Beijing Institute of Civil Engineering and Architecture, 2008 (in Chinese)

[12] Liu S Y. Study on energy storage characteristic in confined aquifer of double-well groundwater heat pump. Dissertation for the Master’s Degree. Dalian: Dalian University of Technology, 2012. 5 (in Chinese)

[13] Zhe C. Operation Optimization Analysis for GWHP systems in heating season. Dissertation for the Master’s Degree. Tianjin: Tianjin University, 2013 (in Chinese)

[14] Hu J H. Research on hydraulics of ground water source heat pump system and temperature field. Dissertation for the Master’s Degree. Changchun: Jilin University, 2009 (in Chinese)

[15] Ni L, Ma Z L, Xu S H, Sun J. Field test of a groundwater source heat pump project with pumping and recharging in the same well in Beijing. HV&AC, 2006, 36(10): 86–92 (in Chinese)

[16] Tsang C F, Buscheck T, Doughty C. Aquifer thermal energy storage: a numerical simulation of Auburn University field experiments. Water Resources Research, 1981, 17(3): 647–658

[17] Heller A J. Heat-load modeling for large systems. Applied Energy, 2002, 72(1): 371–387

[18] Zhu J L, Yao T, Zhang W. Empiric formula study of the fan-coil unit performance on low-temperature geothermal systems. Acta Energiae Solaris Sinica [Tai Yang Neng Xue Bao], 2003, 24(6): 822–826 (in Chinese)

[19] Wu X Z, Zhao J N, Wei J M. Study of heat release calculation method of heating-only fan coil unit. HV&AC, 2010, 40(7): 63–66 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg