Multi‑objective optimization and long‑time simulation of a multi‑borehole ground heat exchanger system

,

GSHPs have received increasing interest because of their high efficiency, low environmental footprint, long life expectancy, low maintenance requirement, and low running costs (Self et al. 2013;Zhu et al. 2014).Despite the popularity of these systems, GSHPs have some drawbacks, such as their high installation costs, infrastructure restrictions, and design limitations, that should be overcome to make these systems more convenient (Goetzler et al. 2012).In this regard, many studies have been conducted on improving the ground heat exchanger (GHE) design and thermal performance of the GSHP systems.These studies are recognized into analytical, experimental, numerical, or optimization categories (Javadi et al. 2019).
Numerical studies are the most common method of investigation for GHEs (Javadi et al. 2019) because of the high cost of experimental studies.Several numerical studies have been conducted on the long-time simulation of GHEs (Hein et al. 2016;Huang et al. 2021;Naicker and Rees 2020;Tang and Nowamooz 2018).Long-time numerical simulations are essential for the performance assessment, maintenance, and life-cycle stability of GSHP systems (Li et al. 2018).GSHP applications generally use multiple borehole heat exchangers (BHE) to enable a higher heat transfer rate (Marcotte et al. 2010;Qian and Wang 2014); therefore, models are expected to perform well for multi-borehole systems.In addition, three-dimensional models simulate such systems more precisely by considering the temperature variations along the depth, capturing interactions between boreholes, and accurately predicting the system performance (Pu et al. 2015;Rees and He 2013).Despite these facts, few studies have conducted long-time three-dimensional numerical simulations for multi-BHE systems because of their complexity and high computational cost.Choi et al. (2018) simulated the 10-year operation of a GSHP consisting of 35 single U-tube BHEs under a load imbalance condition.However, their model was simplified by considering a single BHE instead of the multi-BHE system.Gultekin et al. (2019) investigated the influence of the arrangement and the number of boreholes on the thermal performance of a multi-BHE field.In their long-time simulation, temperature changes along the borehole direction were ignored, and the three-dimensional problem was simplified to a two-dimensional one.The present work aims to bridge the gap by numerically simulating a three-dimensional model of multi-borehole GHE for a long time.
To date, many optimization strategies have been developed for the optimal design of GSHP systems to reduce initial investment costs and improve operating performance.An overview of the recent progress in this area was presented by Ma et al. (2020).Thermodynamic and economic objectives are frequently used for the optimal design of GSHP systems.Among the thermodynamic objectives, entropy generation minimization is widely used in modeling and optimizing GHEs to reduce thermodynamic imperfections from heat transfer and fluid flow irreversibilities (Bejan 1996;Cheng and Liang 2014).The annual, life cycle and total costs are often used as economic objectives (Ma et al. 2020).Several studies (Huang et al. 2014;Kord and Jazayeri 2012;Robert and Gosselin 2014;Soltani et al. 2021) have attempted to optimize GSHP systems using a single objective, economic or thermodynamic.The main problem with the single-objective design optimization strategies is that system performance might need to be improved regarding the overlooked objective (Ma et al. 2020;Sayyaadi et al. 2009).On the other hand, in a multi-objective design optimization approach, thermodynamic and economic objectives are simultaneously optimized, which could improve the optimal design of GSHP systems.Evolutionary algorithms (EAs) in solving multi-objective optimization problems have been greatly expanded in recent years.These algorithms use a population-based approach that evolves a set of solutions in each iteration to eventually provide Pareto optimal solutions well-distributed across the Pareto frontier (Emmerich and Deutz 2018).Sayyaadi et al. (2009) presented a multi-objective optimization method for a GSHP with a vertical GHE to minimize the system's total cost and exergy destruction with eight decision variables using an evolutionary algorithm.Huang et al. (2015) conducted a multi-objective optimization to minimize entropy generation and system upfront cost using a genetic algorithm (GA).By comparing the proposed multi-objective design optimization with a single-objective optimization based on entropy generation minimization for two case studies, it was concluded that multi-objective design could reduce total installation cost by 9.5% and increase energy-saving by 6.2%.Among evolutionary algorithms, the non-dominated sorting genetic algorithm-II (NSGA-II) has widely been used for multi-objective optimization studies (Song et al. 2019;Chang et al. 2023;Mostafazadeh et al. 2023) because of its high solving speed and proper convergence (Deb et al. 2002).Other evolutionary algorithms have rarely been used in energy issues.Employing various EAs in a multi-objective optimization approach based on total cost minimization and exergetic efficiency maximization for the cooling mode of a vertical GSHP system was presented by Keshavarzzadeh et al. (2020).The NSGA-II algorithm provided more acceptable performance in their study than other evolutionary algorithms.
Considering the importance of multi-objective optimization in the appropriate design of GSHPs and the outstanding portion of the GHE in the total cost of the system, this study attempts to determine the optimum value of the multi-borehole GHE design parameters concerning entropy generation and total cost rate as two objective functions.Entropy generation is non-dimensionalized and quantified by entropy generation number (EGN) (Bejan 1996;Maheshkumar and Muraleedharan 2011).A semitwo-dimensional steady-state analytical model describes the GHE thermal behavior in the optimization process.In addition to NSGA-II algorithm, five other evolutionary algorithms are used to define the optimum design parameters: speed-constrained multiobjective PSO (SMPSO), strength Pareto evolutionary Algorithm-II (SPEA-II), Pareto envelop-based selection algorithm-II (PESA-II), multi-objective evolutionary algorithm based on decomposition (MOEA/D), and the third evolution step of generalized differential evolution (GDE3).The previous studies did not apply various evolutionary optimization algorithms to the GSHP system for simultaneous optimization of entropy generation and total cost rate.Five decision variables are considered to find the optimal design of a multi-borehole vertical GHE system.Besides the multi-objective optimization of GSHP, three-dimensional CFD simulation of GHEs assess the long-time performance of the optimized multi-borehole vertical GHE system for 3 months.To the best of the authors' knowledge, a combination of optimization and CFD simulation studies has yet to be conducted.The advantage of the three-dimensional transient CFD simulation is comparing the results with the steady-state analytical solution and getting a better assessment of the optimized GHE system by obtaining features not considered in the optimization.

System description
GeoCool GSHP system located at the Universitat Politènica de València in Spain (Corberan et al. 2011;Montagud et al. 2011) is studied as illustrated in Fig. 1.GeoCool plant is a reversible water-to-water GSHP that provides heating and cooling for a set of spaces with a total area of 250 m 2 .The GHE of this system is made up of six vertical boreholes connected in parallel and arranged in a 2 × 3 rectangular grid over an area of 18 m 2 .Each borehole contains a single high-density polyethylene (HDPE) U-tube, is filled with sand, and is finished with a bentonite layer on top.The present study only considers the cooling operation of this GSHP system.The specifications and design conditions of the system are provided in Table 1.
Three case studies are considered from the optimization of the system named A, B, and C; these cases represent the optimized system from thermodynamic, multi-objective, and economic points of view, respectively.The details of the three case studies will be provided in Sect."Material properties and geometry".

Entropy generation modeling
Entropy generation leads to exergy destruction and loss of available work.Therefore, entropy generation minimization (EGM) is an effective optimization tool to minimize losses and maximize system performance (Bejan 1996;Cheng and Liang 2017).In a vertical GHE, entropy is generated by fluid friction along the tube and finite temperature difference between the fluid and the borehole, resulting in pressure drop and heat transfer losses, respectively (Bejan 1996).Therefore, minimizing entropy generation improves the system's hydraulic and thermal performances.
The following introduces the relations to obtain borehole wall temperature, outlet, and mean fluid temperatures.Next, the equations for pressure drop calculation, total entropy generation, and the economic model are presented.

Heat transfer analysis
The schematic of a single U-tube vertical GHE, along with its main geometric parameters, including borehole radius ( r b ), inner and outer tube radius ( r i , r o ), and borehole shank spacing (2D), is presented in Fig. 2. Two heat transfer processes occur in a BHE: heat transfer between the ground and the grout, and heat transfer between the grout and the fluid through the tube wall.
The analytical model representing the heat transfer process of the vertical GHE consists of the one-dimensional infinite line source model (ILSM) and the two-dimensional thermal resistance of the boreholes acquired from the line source approximation (Lamarche et al. 2010).The ILSM can correctly predict heat transfer in GHE if the operating time is higher than 20r b 2 /α s , which is 54 here (Yang et al. 2010).Also, the operat- ing time is one year.

Design conditions Value
Cooling load (kW) 14 Indoor temperature (℃) 24 Outdoor temperature (℃) 31 Undisturbed ground temperature (℃) 19.6 By considering a long time-scale and constant heat flux in the vicinity of the borehole, Eq. (1) approximates the borehole wall temperature ( T b ) based on ILSM (Huang et al.

2014):
T s,0 , k s , ρ s , c s , and α s are the undisturbed temperature, thermal conductivity, density, specific heat, and thermal diffusivity of soil (i.e., ground), respectively.r b is the borehole radius, t is the time, and E 1 (x) is the exponential integral function.q′ is the heat transfer rate per unit length: where Q is the design load of the GSHP system.N and L b are the number and depth of the boreholes.
The thermal resistance of the borehole determines the heat transfer rate between the circulating fluid and the grout material.The borehole thermal resistance ( R b ) in a single U-tube BHE has three components (Lamarche et al. 2010): (1) ( Convection resistance ( R conv ) and conduction resistance ( R cond ) are defined as Lamarche et al. (2010): and in which, r i and r o are the inner and outer radii of the U-tube.h f is the convection coef- ficient, and k t is the tube's thermal conductivity.Grout's thermal resistance ( R g ) is deter- mined by the line source approximation (Lamarche et al. 2010) where D is the half-shank spacing and k g is the grout's thermal conductivity.
The coefficient of convection heat transfer ( h f ) is obtained through the Nusselt number definitions (Bejan 2013) for fully developed laminar and turbulent flows using the following relationships: and Equation ( 8) is valid for 0.7 < Pr < 160 and Re > 10 4 .Prandtl number and Reynolds number for the flow in the U-tube are obtained from and in which, c f , k f , and µ are the fluid's specific heat, thermal conductivity, and dynamic vis- cosity, respectively.ṁf is the mass flow rate of the fluid in a single BHE.
The thermal resistance between the borehole wall and the soil ( R s ), and the thermal resist- ance for multiple boreholes connected in parallel ( R s,N ) are determined from Lamarche et al. (2010): for uniform heat flux 3.66, for uniform wall temperature (for laminar flows) (8) Nu = 2h f r i k f = 0.023Re 0.8 Pr 0.3 (for turbulent flows). (9) and where B denotes the distance between the boreholes.The arithmetic mean of the fluid's temperature in a BHE ( T f ) is related to the undis- turbed soil temperature, heat transfer rate, and thermal resistances through the following equations: for a single BHE, and through: for a multi-BHE system.
T f is defined as: where T f,1 and T f,2 represent the inlet and outlet temperatures, respectively, for the circu- lating fluid in BHE.

Pressure drop calculation
Pressure drop through a single U-tube is approximated as Li and Lai (2013): where f is the friction factor (Li and Lai 2013): and

Entropy generation rate calculation
The entropy generation rate caused by temperature difference ( Ṡgen, T ) is expressed as (Bejan 1996): in which, T is the difference between T f,m and T b .The parameter χ is the dimensionless temperature difference defined as �T /T f,m , and T f,m denotes the mean temperature of the fluid (Bejan 1996): The entropy generation rate caused by fluid friction ( Ṡgen, P ) is expressed by Eq. ( 21) for an incompressible fluid under non-adiabatic conditions (Bejan 1996): where ṁf,tot is the total mass flow rate through the multi-borehole GHE system.The total entropy generation rate in the GHE ( Ṡgen ) is the sum of the hydrodynamic and thermal entropy generation rates components:

Economic modeling
For economic modeling of the GSHP system, the total revenue requirement (TRR) method presented by Sayyaadi et al. (2009) is employed.Based on the method, the levelized cost rate related to the capital investment and the operation and maintenance expenses of the kth system component ( Żk ), and the levelized cost rate of the electricity provided for the overall system ( Żelec ) are: where τ is defined as the annual operating hours in the cooling mode, and PEC k is the purchase equipment cost for the kth system component.CC L , OMC L , and FC L represent levelized annual carrying charges, levelized annual operation and maintenance costs, and levelized annual fuel costs, respectively; for more details see Sayyaadi et al. (2009).The purchase equipment cost for the GHE is determined by Eq. ( 25) (Sayyaadi et al. 2009), and the purchase equipment cost for the heat pump's components is provided in Table S1 of the Supplementary Material: L t is the U-tube length.C t is the cost of the U-tube per meter, set to be 1.6 €/m.C b , the drilling and grouting cost per meter is 50

Optimization
Figure 3 shows a block diagram of the optimization strategy, including three main steps used for the multi-borehole GHE systems.After the energy and economic analysis of the system, the decision-making variables are selected among all the design parameters based on sensitivity analysis.The decision-making variables affect the objective functions more strongly than the other design parameters.The next step is formulating the optimization method.After developing the mathematical model of the GSHP system, the objective functions, decision variables, and optimization constraints are defined.The decision variables are then optimized through a variety of multi-objective evolutionary algorithms, which search for a set of Pareto optimal solutions.

Objective functions, decision variables and constraints
Two objective functions are considered: EGN and total cost.The optimization algorithms are used to minimize both objective functions simultaneously.EGN ( N s ) describes system irreversibilities is defined as: The total cost of the GSHP system ( Ċtot ) is expressed as: Q .

Fig. 3 Outline of the optimization methodology
Five influential design parameters, including borehole depth ( L b ), borehole radius ( r b ), number of boreholes (N), circulating fluid mass flow rate per U-tube ( ṁf ), and U-tube outer radius ( r o ) have been selected as the decision variables in this study based on the results of a global sensitivity analysis conducted by Huang et al. (2014); these five decision variables are determined based on their sensitivity indices in terms of the EGN.It is worth mentioning that GHE-related costs are a significant part of the GSHP total cost.L b and N are the most influencing parameters on the cost of the GHE/GSHP system (Robert and Gosselin 2014).
A reasonable range should be used for the values of the decision variables to avoid unrealistic solutions.These constraints are determined based on the recommendations obtained from practical engineering projects (ASHRAE 2011; Banks 2012), shown in Table 2.The values of the low-sensitive parameters, which are considered constant in the optimization process, site-related parameters, and grout thermal conductivity are presented in Table 3.
Two other constraints are used for the heat transfer rate per unit of length ( q′ ) and outlet temperature.The acceptable range for q′ is between 30 and 130 W/m proposed by Robert and Gosselin (Robert and Gosselin 2014).The recommended values for the minimum outlet temperature ( T f,2,min ) and maximum outlet temperature ( T f,2,max ) are deter- mined by Eqs. ( 29) and (30) for the heating and cooling modes, respectively (ASHRAE 2011):  where T s,min and T s,max are the minimum and maximum soil temperatures over a year.

Multi-objective evolutionary algorithms
NSGA-II is a multi-objective evolutionary algorithm based on genetic algorithm.
In addition to the NSGA-II algorithm, five other evolutionary algorithms, including SMPSO, SPEA-II, PESA-II, MOEA/D, and GDE3, are applied and compared.SMPSO incorporates particle swarm optimization (PSO) and Pareto dominance to handle multi-objective optimization problems.This algorithm can store non-dominated solutions found in a search in an external repository.It could enhance the exploratory ability of PSO by utilizing a time-dependent mutation operator (Nebro et al. 2009).SPEA-II is an improved version of SPEA.This algorithm behaves similarly to NSGA-II (Zitzler et al. 2001).PESA-II is a developed version of PESA.This algorithm uses a new selection method based on regions instead of individuals.The mechanism of this method is very similar to multi-objective particle swarm optimization by the difference that PESA-II uses genetic operators instead of PSO operators.PESA-II has an external archive for non-dominated solutions; parents and mutants are selected from the created regions in this archive (Corne et al. 2001).The structure and mechanism of MOEA/D algorithm are different from other algorithms.It decomposes multi-objective optimization problems into some scalar optimization sub-problem and optimizes them simultaneously.The number of sub-problems is usually equal to the size of the population.The method assigns a weight vector to each sub-problem.The evaluation is performed based on this vector.Finally, each sub-problem theoretically provides a Pareto Front solution at the end of the search (Zhang and Li 2007).GDE-3 is an extension of differential evolution that can support an arbitrary number of objectives and constraints.It provides a better distribution solution than the previous GDE version (Kukkonen and Lampinen 2005).

Model validation
The BHE was optimized using the entropy generation minimization method and genetic algorithm.A comparison between the results of Huang et al. (2014) and those from the present study is shown in Table 4.The aim of the optimization is to enhance the thermodynamic efficiency of the original system design (considered an unoptimized design) by minimizing EGN through modifying decision variables.The EGN, the optimization (30)

Numerical simulations
From the optimization results of the studied GSHP system (to be presented in Sect."Optimization results"), three prominent points labeled as points A, B, and C in Fig. 8 will be selected for CFD simulation as case studies.Transient 3D models of the selected case studies are simulated long-term.The model geometry is divided into four subdomains: the circulating fluid, the U-tube walls, the grout, and the ground.The following assumptions are used for the simulations: -The ground is an infinite homogeneous medium for heat transfer.
-The thermophysical properties of the ground, grout, U-tube walls, and fluid are constant and independent of temperature.-The effects of underground water flows are neglected.
-The circulating fluid is Newtonian and incompressible.
-Heat losses due to radiation, convection, and evaporation from the ground surface are neglected.-The temperature of the ground surface is considered uniform and constant.
Transient fluid flow and heat transfer in the fluid and solid subdomains are simulated by solving the three-dimensional and unsteady governing equations: conservation of mass, conservation of momentum, and conservation of energy.
conservation of energy for the fluid subdomain (Pletcher et al. 2012) is where e is the internal energy per unit mass, Q is the heat produced per unit volume by external agencies and ϕ is the viscous dissipation term.Conservation of energy for the solid subdomains (Pletcher et al. 2012) is in which, ρ i , c i and k i are the density, specific heat, and thermal conductivity of one of the solid subdomains, respectively.The solid subdomains include the ground, grout, and U-tube walls.
The Reynolds numbers in the present study are higher than 10 4 ; therefore, the realizable k−ε model is selected as a turbulence model (Serageldin et al. 2018).
The finite-volume discretization method is used to approximate the governing equations.The numerical calculations are carried out using a double-precision and pressurebased solver.The SIMPLE algorithm is selected for pressure-velocity coupling, and a second-order discretization scheme is used to discretize pressure, momentum, energy, (31) k and ε equations.The convergence criteria are set to be 10 −8 for energy equation, and 10 −5 for the rest.The time step is considered 300 s with a total number of 26,574 steps.

Material properties and geometry
The simulated geometry consists of four subdomains: one fluid subdomain, which includes the water flowing inside the U-tube, and three solid subdomains, including the U-tube walls, grout, and the ground.A high-density polyethylene (HDPE) U-tube is placed inside the borehole.The space between the tube and the ground is filled with grout material, combining sand and bentonite.The thermal properties of the materials used in the fluid and solid subdomains are shown in Table 5.The geometric parameters and the size of the computational domains for cases A, B, and C are shown in Table 6.It should be noted that cases A, B, and C represent three prominent points of A, B, and C on the Pareto frontier presented in Fig. 8 and detailed in Table 10.Therefore, the values of the design parameters, i.e., N , L b , r b , and r o , for cases A, C, and B have been deter- mined by optimization to minimize entropy generation number, total cost rate, or a combination of both targets, respectively.Other parameters not subject to optimization, such as D , B , and t , are constant across three cases.In all three cases, the boreholes are connected in parallel, with a distance of 3 m from each other.These arrangements closely align with the original design described in Sect."System description".Figure 4 illustrates the arrangements of boreholes in cases A, B, and C. As an example, the computational domain for Case B and the essential geometric parameters are shown in Fig. 5.The outermost boundaries of the soil subdomain are considered far enough from the boreholes to prevent the boundary conditions from affecting the results.

Mesh generation
The computational domain is divided into subdomains with two types of threedimensional meshes to reduce the number of cells, shown in Fig. 6 for Case B. Other cases have a similar mesh structure.An unstructured mesh is considered.Tetrahedral cells are applied for meshing the 180° bend of the U-tube and the base of the grout and the ground.The sweep meshing method is used for the vertical parts of the tubes, plus the grout and the surrounding ground.A two-dimensional horizontal cross-section of the domain was firstly discretized into the triangular cells, swept in the vertical direction to generate a three-dimensional mesh, and divided into hundred slices to capture variations in the vertical direction.The cell sizes gradually decrease as moving from the outer boundaries toward the center.The smallest cell size of 2 mm is used for the 180° U-tube bend, where strong temperature gradients are expected.An optimal mesh size is obtained by performing a mesh independency test for the three case studies under steady-state operating conditions.The outlet temperatures of Case B are presented in Table 7.The optimum is mesh#3 with 1.65 × 10 6 cells because the temperature's relative deviation for a finer mesh (mesh #4) is only about 0.1%.With a similar method, the optimum number of cells for Case A and Case C are determined to be 4.42 × 10 6 and 0.73 × 10 6 , respectively.

Initial and boundary conditions
The initial temperature of the whole computational domain is equal to the undisturbed ground temperature, i.e., 19.6 • C .In Table 8, a Dirichlet boundary condition with a constant temperature of 19.6 • C is assigned to the upper and bottom surfaces of the ground.In contrast, the adiabatic condition is used for the ground's lateral surfaces.Constant-mass-flow and constant-pressure boundary conditions are used for the inlets and outlets of the fluid subdomain, respectively.An additional time-varying Dirichlet boundary condition is also used for the inlets, which describes the incoming fluid's temperature changes as a function of the heat transfer between the circulating fluid and the heat pump.The following equation expresses the time-varying inlettemperature boundary condition: The heat rate ( Qb ) assigned to each borehole, for cooling operation: (33)   Ground lateral surfaces Adiabatic ( q = 0)

Flow inlet
Constant mass flow rate (different for each case study) and time-varying Dirichlet (Eq.33) Flow outlet Pressure-outlet ( P gauge = 0)

Tube wall
No-slip COP is the average coefficient of the heat pump's performance in the cooling mode, equal to 5 (Nam et al. 2008).
According to Eqs. ( 33) and ( 34), the inlet fluid temperature at each time step is a function of the cooling load and the outlet fluid temperature from the previous time step.

Model validation
An experimental dataset from the 5-year operation of GeoCool GSHP (Ruiz-Calvo and Montagud 2014), Table 9, has been applied to validate the long-term transient simulation.Despite a few quantitative differences, this experimental study (Ruiz-Calvo and Montagud 2014) closely matches the physical processes, setup, and operating parameters used in the present numerical model, making it a suitable case for validation.
The system was monitored from February 2005, working 15 h daily and five days weekly.Data were collected every minute.The numerical model uses collected data for each (34)  borehole's inlet water temperature and mass flow rate as boundary conditions.The measured data for outlet water temperature is used for validating the model.Figure 7 compares the average outlet flow temperature of six boreholes during April 2005 in the heating mood.The model predicts the outlet temperature with a maximum error of about 4.5%, showing the acceptable performance of the model.The numerical and experimental outlet flow temperature values for each borehole are compared in Fig. S1 of the Supplementary Material.The difference in outlet temperature between the numerical and experimental studies can be attributed to a 20% uncertainty in estimating the ground thermal conductivity, as reported in Ruiz-Calvo and Montagud (2014).Additional contributing factors include experimental errors due to environmental conditions, as well as the accuracy of the data acquisition system and the temperature sensors.

Results and discussion
This study presents a multi-objective optimization of the GSHP system used in an institutional building in Spain and a 3D CFD simulation of its multi-borehole GHE.The optimization results, presented in Sect."Optimization results", are obtained based on a steady-state analytical model using NSGA-II and five other evolutionary algorithms.The simulation results of the same system, presented in Sect."Numerical simulation results", study the outlet fluid temperature distribution, the borehole wall temperature distribution, and the transient entropy generation have been studied for three cases over a long time.Finally, a comparison between the optimized and simulation results is presented.

NSGA-II optimization
Figure 8 presents the Pareto frontier solution of the proposed GSHP system obtained by minimizing EGN (Eq.27) and total cost rate (Eq.28) simultaneously, using NSGA-II multi-objective optimization algorithm.An increase in total cost accompanies a decrease in the EGN.As illustrated in Fig. 8, the minimum EGN occurs at point A (0.040), where the total cost rate is maximum (37.776 €/h).On the other hand, point C has the highest value of the EGN (0.192), while its total cost rate has the smallest value (3.329 €/h).
Points A and C represent the optimal design points in single-objective optimization, in which EGN and total cost rate are the objective functions, respectively.In multi-objective optimization, selecting the desired optimal point among a set of optimal solutions on the Pareto frontier requires a decision-making procedure based on potential design limitations or financial considerations.One of the procedures for determining the desired optimum point is to use a hypothetical point indicated as the ideal point, Fig. 8 (Navidbakhsh et al. 2013).At this point, each objective has its optimal value and is independent of the other one; the ideal point cannot be located on the Pareto frontier because it is impossible to reach the optimum values of both objectives at the same time.Therefore, the closest point on the Pareto frontier to this ideal point is selected as the desired optimal point (Huang et al. 2015;Keshavarzzadeh et al. 2020;Ghanadi Arab et al. 2013).This point, presented as point B in Fig. 8, has an EGN of 0.079 and a total cost rate of 13.462 €/h.Point B is considered a trade-off between EGN and total cost rate.The calculated results for the three prominent points (points A, B, and C) and the original design data are summarized in Table 10.EGN in points A and B are 77.4% and 55.4% less than the EGN in the original system; whereas for point C, EGN is 8.47% higher than that for the original design.From an economic point of view, the total cost rate for points A and B is 4.6 and 1.6 times more than the original case, respectively.However, the total cost rate in point C is 59% less than that of the original system.Therefore, the original system design is close to point C, which is the optimal design point based on the economic single-objective optimization of the GSHP.However, the multi-objective design optimization shows that by increasing the total cost rate of the system by just 5 €/h, a significant reduction in the entropy generation number is achieved, which improves the thermal performance of the GHE.
Figure 9 presents scatter distributions of five decision parameters along the Pareto frontier, which illustrates the trends of each parameter and their impact on the GHE performance.It can also aid designers in the decision-making process.Figure 9a indicates that the borehole radius is nearly constant and located near its lower bound for all Pareto optimal points because an increase in the borehole radius increases the thermal resistance between the fluid and the ground.In Fig. 9b, the distribution of the number of boreholes is scattered uniformly between 3 and 10 because of the direct impact of the number of boreholes on both objective functions.Although increasing the number of boreholes of the GHE is desirable from a thermodynamic point of view, it leads to an increase in the total cost rate of the system similar to borehole depth, Fig. 9c.Although the spread of the optimal values of borehole depth is relatively uniform, they are mainly located near the upper and lower bounds.According to Fig. 9d, the optimum mass flow rate of the fluid in this study is around 0.4.In Fig. 9e, the distribution of the U-tube outer radius is located near the upper bound for almost all optimum solution points.Larger tubes reduce the overall thermal resistance of the GHE and improve the heat transfer rate between the fluid and the ground.

Comparison of evolutionary algorithms
The Platypus library of Python was employed to compare the performance of different EAs in multi-objective optimization of the GHE.Python allows running several EAs simultaneously with the same random number generator seeds while saving CPU time with parallelization.For all algorithms, the direction of the objective functions is to minimize both EGN and total cost rates.The selection of the best optimization algorithm among various options is based on its ability to achieve a uniform distribution of optimal points and ensure diversity along the Pareto frontier (Keshavarzzadeh et al. 2020(Keshavarzzadeh et al. , 2019;;Chen et al. 2010;Li et al. 2014).Figure 10a shows the NSGA-II Pareto frontier discussed in the previous subsection, demonstrating excellent diversity and continuity in the distribution of optimal points.The Pareto frontier of GDE-3, Fig. 10b, is similar to that of NSGA-II, and there is not any considerable difference between them.The PESA-II Pareto frontier, Fig. 10c, is almost identical to the previous two algorithms, although it is not as uniformly distributed as those generated by GDE-3 and NSGA-II.SPEA-II algorithm yielded a discontinuous Pareto frontier, Fig. 10d, particularly for EGNs above 0.12 and Similarly, the SMPSO algorithm results in a discontinuous Pareto frontier, as shown in Fig. 10e, with some data missing in two subranges.Finally, Fig. 10f clearly shows that MOEA/D algorithm provides an unstable behavior, failing to provide a well-defined Pareto frontier.Considering the notably diversity and uniform dispersion of optimal solutions along the Pareto frontier achieved by the NSGA-II and GDE-3 algorithms in contrast to the other algorithms, these two algorithms exhibit superior performance and reliability.These observations are consistent with earlier studies (Keshavarzzadeh et al. 2020(Keshavarzzadeh et al. , 2019) ) comparing evolutionary algorithms where the NSGA-II algorithm provided acceptable results.The comparison emphasizes the capability of NSGA-II as a dependable and efficient algorithm for addressing the complexities of multi-objective optimization in energy issues.

Temperature distributions
Variations of the average fluid temperature at the outlet of the GHE during three months of cooling operation for the three case studies corresponding to points A, B, and C on the NSGA-II Pareto frontier (Fig. 8) are shown in Fig. 11.The cooling load for all three case studies was 14 kW.Case C has the highest outlet fluid temperatures, followed by Case B. This difference between outlet temperatures is mainly due to the difference between the total borehole lengths (= borehole number × borehole depth).According to Table 6, the total borehole length for Case B is about one-third of the corresponding value for Case A and 3.9 times more than that for Case C. Shorter total borehole lengths lead to an increase in the heat flux per meter per borehole, which increases the outlet temperature of the circulating fluid in the GHE.Increasing the temperature of the fluid entering the heat pump system in the building increases the energy consumption of the heat pump and decreases its exergetic efficiency.resistance for the total length ( R * b ) in cases A, B, and C are calculated as 9.52 × 10 −5 , 13.69 × 10 −5 and 27.38 × 10 −5 K/W , respectively.According to the simulation results of the temperature distributions in the GHEs, case study A demonstrates the best thermal performance, and case study C gives the weakest thermal behavior.These observations align entirely with the optimization results provided in Sect."Optimization results", Fig. 8; among all the optimal points across the Pareto frontier, point A emerges as the most optimal point from a thermodynamic point of view, and point C presents the weakest point from this aspect:

Transient entropy generation
Transient entropy generation is calculated using Eq. ( 22) as a function of time for the three case studies.From Sect."Analytical modeling for multi-borehole GHE", the total entropy generation rate ( Ṡgen ) in a GHE is the sum of the thermal ( Ṡgen, T ) and hydro- dynamic ( Ṡgen, P ) entropy generation rates caused by heat transfer and fluid friction.Variations of all entropy generation rates during the 3-month operation of the system for case studies A, B, and C are demonstrated in Fig. 14a-c, respectively.In all three studies, Ṡgen, P remains almost constant and has a small contribution to the total entropy generation rate in the GHE.The average values of Ṡgen, P for cases A, B, and C are 4.69, 1.48, and 0.34 W/K , respectively.According to Eq. ( 16), the pressure drop along the U-tube ( P ) is directly related to the borehole length.Therefore, reducing the borehole This reduction in pressure drop, combined with reducing total mass flow rate ( N × ṁf ) from Case A to Case C, reduces the entropy generation rate caused by friction based on Eq. ( 21).According to Fig. 14, the thermal entropy generation component has the major contribution to the total entropy generation rate for all cases during the entire system operation.The trends of the Ṡgen, T variations are relatively similar for all three case studies.Ṡgen, T increases for a short period after the system starts operating because of the increase of T .With the continuing growth of T f,m while T is constant, Ṡgen, T slowly decreases with different slopes.The variations of T f,m and T over simulation time for case studies A to C are shown in Figs.S2 and S3 of the Supplementary Material.Due to the highest T in Case C, the Ṡgen, T during the entire system operation has the highest amount, followed by Case B. So, with the decrease in the total borehole length from cases A to C, despite the reduction of the Ṡgen, P value, the value of Ṡgen, T is increased.Since the temperature difference factor has a more substantial effect than the fluid friction factor on the total entropy generation rate, Ṡgen is also increased by decreasing the total borehole length.Figure 15 compares the transient EGN variations for the three cases.After a slight increase of EGN for the early times of running, EGN gradually decreases and converges to a constant value.EGN is highest for Case C, with values ranging from about 0.11 to   11.The trend of EGN variations for the three case studies is the same with both methods.The result acquired for case studies A and B agree to be within about 9.0% and 7.6%, respectively.There is a 29.7% difference between the two data sets for Case C.These quantitative differences are expected because the two methods are fundamentally different.

Conclusion
The design of a ground-source heat pump with multi-borehole ground heat exchangers (GHE) for an educational building located in Spain is optimized using a variety of multi-objective evolutionary algorithms to minimize entropy generation number (EGN) and total cost rate simultaneously.In addition, the optimized multi-borehole GHEs are numerically simulated using a three-dimensional transient model for three months to study the variations of their thermal performance over time.
• Based on the NSGA-II Pareto frontier, the desired optimum values of EGN and total cost rate are 0.079 and 13.46 €/h, respectively (point B).Absolute EGN optimization yields the optimum solution at 0.040, where the total cost rate is at its highest value, 37.77 €/h (point A).Additionally, absolute economic optimization gives the optimum solution at 3.33 €/h, where the EGN is at its highest value, 0.192 (point C).The multi-objective solution achieves a remarkable 55.4% reduction in ENG while experiencing a 65.1% increase in the total cost rate compared to the original case • The selection among points A, B, and C as the final optimized solution depends on the priorities and specific requirements of decision-makers and designers.If the primary concern is maximizing the thermodynamic performance of the GHE regardless of costs, then point A is the optimal choice.Conversely, point C is most suitable if the focus is just on minimizing costs.For scenarios where a balance between thermodynamic efficiency and cost-effectiveness is desired, point B presents the ideal one.• NSGA-II and GDE-3 evolutionary algorithms perform superior to PESA-II, SPEA-II, SMPSO, and MOEA/D algorithms in generating the Pareto-optimized solutions.MOEA/D algorithm is not an acceptable algorithm in this case.• The optimal values of the number and depth of boreholes are scattered uniformly between the specified lower and upper bounds; the optimum state is obtainable through various values between the defined bounds.The lower borehole radius and the higher tube radius contribute to attaining the optimized condition of the system.• The trend of entropy generation variations obtained from the 3D transient numerical simulation is similar to the semi-2D steady-state analytical model, which affirms the efficiency of the optimization approach in determining the optimal configuration of BHE.Consequently, this approach provides designers with a reliable tool for choosing appropriate design parameters to achieve an optimum design.• Simulation results show that reducing the total borehole length lessens the entropy generation caused by fluid friction while increasing the entropy genera-tion caused by temperature differences.The temperature difference has a more substantial effect than the fluid friction on the entropy generation in the BHE, so the total entropy generation increases by decreasing the total borehole length.

Fig. 1
Fig. 1 Schematic of GeoCool GSHP.The heat pump is connected to an internal and external circuit.The internal circuit consists of 12 parallel-connected fan coils, a water storage tank, a circulation pump, and an internal hydraulic loop.The external circuit includes a multi-borehole GHE, a circulation pump, and an external hydraulic loop

Fig. 2
Fig. 2 Schematic of a single U-tube BHE

Fig. 5
Fig. 5 Computational domain and geometric properties of Case B: a entire GHE model consisting of 5 parallel boreholes; b upper view of each borehole, cut in half; c bottom view of each borehole, cut in half

Fig. 8
Fig.8The NSGA-II Pareto frontier; points A and C represent the optimal solutions of single-objective optimization with the objective function of EGN and total cost rate, respectively.Point B represents a trade-off between the two objective functions in multi-objective optimization

Fig. 9
Fig. 9 Scatter distribution of design parameters along the Pareto frontier: a borehole radius, b borehole number, c borehole depth, d fluid mass flow rate, e U-tube outer radius.The distributions of the optimum borehole number and depth are scattered, indicating that the optimum situation could be obtained with various values of these parameters, depending on the designer's decision

Fig. 10
Fig. 10 Pareto optimal frontier of different evolutionary algorithms: a NSGA-II, b GDE-3, c PESA-II, d SPEA-II, e SMPSO, f MOEA/D.NSGA-II and GDE-s algorithms perform best in obtaining Pareto optimal solutions, thus selected as the most suitable algorithms

Fig. 11
Fig. 11 Average outlet fluid temperatures over time for three cases: Case A: single-objective thermodynamic optimum design, Case C: single-objective economic optimum design, Case B: multi-objective thermo-economic optimum design

Figure 12
Figure 12 demonstrates the variations in the average temperature of the borehole wall.Similar to the results for outlet temperature, Case C has the highest borehole wall temperature, followed by Case B. Reducing the total borehole length from Case A to Case C increases the heat flux per meter of the borehole.The difference between the arithmetic average fluid temperature ( T f ) and the average borehole wall temper- ature is plotted in Fig. 13.After about three days, the temperature of the borehole wall and circulating fluid rises at the same rate, and the difference between them reaches an almost constant value of 1.6, 2.3 and 4.6 ℃ for cases A, B, and C, respectively.The highest temperature difference for Case C indicates the most significant thermal resistance between the fluid and the boreholes, leading to the lowest heat transfer rate.Using numerical results and Eq.(35), the values of the overall borehole

Fig. 12
Fig. 12 Average borehole wall temperatures for three cases over time: Case A: single-objective thermodynamic optimum design, Case C: single-objective economic optimum design, Case B: multi-objective thermo-economic optimum design

Fig. 14
Fig. 14 Thermal, hydrodynamic, and total entropy generation rates in the GHE during three months of system cooling operation: a Case A, b Case B, c Case C

Fig. 15
Fig. 15 Transient EGN for three cases: Case A: single-objective thermodynamic optimum design, Case C: single-objective economic optimum design, Case B: multi-objective thermo-economic optimum design

Table 4
Huang et al. (2014) et al. 2014)optimization results ofHuang et al. (Huang et al. 2014)objective function, is compared for validation purposes.The comparison indicates that the method used herein provides a lower EGN, by 2.38%, compared to that obtained byHuang et al. (2014).The optimized values of the decision variables in this study and that ofHuang et al. (2014)are slightly different, which stem from the different optimization models utilized in these studies.

Table 6
Geometric properties and computational domain size for three cases Cases A, B, and C correspond to points A, B, and C presented in Fig.8and Table10, reflecting a distinct optimization perspective.The values of the design parameters, including N , L b , r b , and r o , for each case, have been determined from the optimization results (to be presented in Sect."Optimization results")

Table 7
Mesh-independency results for Case B with five boreholes

Table 8
Boundary conditions

Table 10
Design parameters and objective functions for the original design and points A, B, and C on the Pareto frontier

Table 11
Comparison of the EGN values between simulation and optimization studies

study EGN at the end of the 3-month simulation EGN obtained from optimization
, and lowest for Case A, with values below about 0.05.EGN values for Case B lie between 0.04 and 0.08.A comparison between the EGN values at the end of the 3-month simulation and the EGN values obtained from optimization for three cases is presented in Table