Modelling of waste heat recovery of a biomass combustion plant through ground source heat pumps- development of an e ﬃ cient numerical framework

(cid:129) New framework has been developed for ground source heat pump investigations. (cid:129) Biomass waste heat recovery through ground source heat pump is investigated. (cid:129) Di ﬀ erent possible storage temperatures are analyzed. (cid:129) Temperature variations in ground around heat exchangers have been explored. (cid:129) E ﬀ ectiveness of the heat recovery system is presented. Development of a reliable and convenient dynamic modelling approach for ground source heat pumps remains as an important unresolved issue. As a remedy, in this work a novel, computationally-e ﬃ cient modelling framework is developed and rigorously validated. This is based upon an implicit computational modelling approach of the ground together with an empirical modelling of heat and ﬂ uid ﬂ ow inside U-tube ground heat exchangers and waste heat calculations. The coupled governing equations are solved simultaneously and the in ﬂ uences of parameters on the performance of the whole system are evaluated. The outcomes of the developed framework are, ﬁ rst, favorably compared against two di ﬀ erent existing cases in the literature. Subsequently, the under- ground storage and recovery process of the waste heat through ﬂ ue gases generated by a biomass combustion plant are modelled numerically. This reveals the history of temperature distributions in the ground under dif- ferent con ﬁ gurations of the system. The results show that for a biomass combustion plant generating ﬂ ue gases at 485.9K as waste heat with the mass ﬂ ow rate of 0.773kg/s, the extracted heat from the ground is increase by 7.6%, 14.4% and 23.7% per unit length of the borehole corresponding to 40°C, 50°C and 60°C storage temperatures. It is further shown that the proposed storage system can recover a signi ﬁ cant fraction of the thermal energy otherwise wasted to the atmosphere. Hence, it practically o ﬀ ers a sizable reduction in greenhouse gas emissions.

usually increase by 10-30% due to the reduced heat transfer efficiency between the high-temperature flue gas and heat exchanges [6]. It follows that there is a significant potential to recover the waste heat of the flue gases from industrial boilers [6]. This is of higher importance in the case of a biomass-based boilers as the waste heat is essentially carbonfree. More importantly, the limitations in biomass resources further necessitate saving biomass and therefore minimizing heat losses and recovery of waste energy.
There have been already studies on the waste heat recovery from flue gases. The early investigations proposed to increase the heat exchanging surface areas of air or water preheating, but it is limited to space constraints or high cost. Another method for the flue gas heat recovery involves applying gas/gas heat exchangers (GGHs) in which cleaned flue gas is heated by the uncleaned flue gas, decreases the temperature of the uncleaned flue gas while maintains the temperature of the cleaned flue gas for venting. The technology however does not have any effect on energy saving [7]. Organic Rankine Cycle (ORC) [8,9] is also a practical way to recover the exhaust waste heat. The lowgrade energy in the exhaust flue gas is used to generate high-grade energy with ORC system, which improves the combined system efficiency. However, the high cost and complex operation of the integrated systems are the disadvantages [8]. As a result, the methods of waste heat recovery from combustion systems are still under development and constantly call for more research.
Ground has been found to be an excellent medium for storing heat for a long time with a relatively low cost due to its proper heat capacity [10,11]. Thermal energy can be stored in the ground with the applications of Ground Heat Exchangers (GHEs) made of high-density polyethylene (HDPE) pipes with different shapes such as boreholes with U-tube pipes, slinky, spiral etc. As a general rule, polyethylene (PE) pipe for pressure applications can be safely used for temperatures as low as −40°C and as high as 60°C [12]. For non-pressure service, the allowable temperature range widens up to 82°C. There are a few PE piping materials that have qualified for a pressure rating at 82°C [12]. In the borehole thermal energy storage method, heat is transferred from a high temperature medium to the ground through vertical boreholes with U-tube pipes and then stored in the ground [10]. The volume of the storage area is not exactly separated. Thermal properties of the ground play an important role in specifying the thermal capacity of the storage medium [10]. Ground allows for a clean, sustainable and costeffective heat storage [13,14]. Further, due to the unlimited boundaries of the ground, it leads to a stable temperature distribution mostly unaffected by the seasonal variations [15]. A geothermal energy extraction system can be applied to transfer ground thermal energy to buildings [16] with little cost compared to other sources [17]. Ground source heat pump (GSHP) system is one of the most promising applications in supplying heating/cooling to buildings due to its high efficiency and environmental friendliness [18]. Compared with traditional air conditioning systems, GSHP system has 25-50% less power consumption by utilizing the stable ground energy [15].
Applying GSHP system for waste heat storage and recovery from industrial boilers, especially biomass plants have great potential for waste heat usage. The technology used for the storage process is typically dependent on the temperature of waste heat [19]. However, the performance of GSHP system, which can be predicted numerically or analytically, need to be further investigated. In doing so, the ground and heat pump should be analyzed separately. Several studies have modelled GSHP systems of waste heat utilization. These are reviewed briefly in the following. Recently, industrial waste heat storage process using large scale heat storage medium was explored by Moser et al. [20]. Their research focused on the case study of the industrial city of Linz (Austria) and advantages and disadvantages of seasonal heat storage were discussed vastly. The results indicate that the number of annual cycles is crucial for a seasonal heat storage. In Dehghan's work [10], waste heat from micro gas turbine exhaust gases was stored in ground through spiral (helical) GHEs in the ground and then recovered the GSHP system. The process was simulated by COMSOL and the results show that amount of extracted heat from the ground is considerably increased after the waste heat storage process. Furthermore, for the storage process, the optimum distance between GHEs was calculated to be 7 m. Kim et al. [15] improved the efficiency of GSHP system by coupling with a heat storage tank. The analyses on the operational behavior of heat storage tank and the variations of heat pump energy performance due to the connection with heat storage tank were performed in detail. They found that, compared with traditional air conditioning systems, the energy saving for cooling and heating are An integrated simulation framework about GSHP system was developed by Rui et al. [21]. They modelled all parts of GSHP including heat pump components, ground and GHEs separately and integrated simulations of each part to introduce a model for the whole system. A finite element model was developed which could predict temperature variations around GHEs. Esen et al. [22] analyzed thermal performance of a GSHP system based on an adaptive neuro-fuzzy inference system with a fuzzy weighted pre-processing method in MATLAB. Their method was shown to be faster and simpler than the existing ones. Further, to investigate thermal performance of the GSHP systems, different approaches such as artificial neural networks [23], vector machine [24] have been presented in the literature. A new numerical modelling approach for the transient temperature distribution in soil around GHEs was presented by Kayaci et al. [25]. Numerical results were validated against experimental data and the long-term transient soil temperature profile was obtained numerically with realistic boundary conditions using meteorological data. Their results indicate that for small burial depths (< 1 m), the effects of surface conditions can be seen earlier and therefore the soil temperatures are higher than for burial depth of 2 m and 3 m. Also, when the GSHP system operates for a long heating period, the average soil temperature in the solution domain and soil temperature around the pipe decrease by approximately 2.5°C, 0.7°C when Q max = 21 W/m, respectively. Schmidt et al. [26], investigated central solar heating plants with seasonal storage in Germany. They described seasonal heat storage system as well as solar heating plants in detail and concluded that, through integration of the seasonal heat storage, more than 50% of the annual heating demand for space heating and domestic hot water can be supplied by solar energy. Liu et al. [27] designed an experimental demonstration of solar seasonal storage coupling with GSHP system. In their heat storage experiment, a system with 1500 m 2 solar thermal collectors and 580 sets of 120 m deep GHEs was considered. The results showed that the system could improve the efficiency of solar energy utilization up to 50.2%.
Central to the wide application of underground heat storage, is analysis and simulation of the storage medium to aid the design process. However, simulation of heat storage in the ground by using computational packages such as COMSOL [28][29][30] and ANSYS or through analytical methods such as Green's function [31,32] method can be time-consuming and often involves complex procedures. Therefore, novel, efficient and reliable framework should be developed to model thermal performance of GSHP systems. To achieve these goals, ground storage of the thermal energy of exhaust, gasses generated by combustion of biomass, is modelled numerically. The heat system under investigation includes U-tube GHEs and storage during warm season followed by heat extraction in cold seasons. The current work puts forward a novel accurate and yet simple numerical framework for evaluation of GSHP system performance by using Engineering Equation Solver software (EES). In comparison with other simulation techniques, the developed framework can solve complex problems of GSHP modelling faster and more accurately.

Biomass combustion plant
The investigated combustion system includes a small biomass thermal plant for the purpose of water heating. Biomass fuel is mixed with air in the burner and combustion takes place with the assistance of ignition devices in the combustion chamber [4]. Part of the heat generated by biomass combustion is transferred to the water through heat exchangers, whilst the rest are wasted chiefly through the flue gases. The hot water produced by the boiler is delivered to various thermal equipment including hot water coils in air handling units, service hot water heating equipment and terminal units. The mass flow rate ṁf uel (kg/s) and temperature T f (K) of flue gas are determined by the boiler thermal output P (kW), boiler thermal efficiency η t , flue gas heat loss percentage η f , excess air coefficient x air and the lower heating value of the fuel LHV (kJ/kg), which can be obtained from an ultimate analysis.
Therefore, the heat of combustion and boiler output is expressed by the following equations [33], The flue gas heat loss percentage is defined by η f and therefore the energy balance for exhaust flue gas can be expressed by [33], In Eq.
(2), T 0 is the temperature of air at the boiler inlet and c p flue . is the specific heat of flue gas.
According to the mass conservation of chemical reactions and boiler system, the mass flow rate of exhaust flue gas is the sum of feed fuel and supplied air, which can be expressed as follow, Here, m (̇) air s is the required air flow rate for stoichiometric combustion of biomass fuel burnt.

Integration of underground thermal energy storage and heat recovery systems
The waste heat from the flue gases can be utilized beneficially by storing it in the ground and recovering that through a GSHP system during cold seasons. The recovered heat can be used to meet the thermal demands of buildings. Given the carbon neutrality of the biomass combustion as the source of energy, the resultant thermal technology will be a low-carbon one and is therefore environmentally benign. The key parameters for finding the heat transfer and losses involved in the storage process are thermal properties of the storage medium, duration and temperature of storage as well as storage geometry and volume [10,34]. High temperature waste heat coming from biomass thermal plant is cooled in a heat exchanger by an external fluid loop. Controllers 1 and 2 control temperature of circulated fluid, which should not exceed 90°C (storage temperature, T s ) at the entrance of GHEs. Fluid is pumped into the ground through 50 m depth U-tube GHEs at constant T s for three months continuously for waste heat storage. In this study, U-tube GHEs, as shown in Fig. 1, are employed to investigate the effects of waste heat storage on the performance of GSHP system.
After completing the storage process, the stored heat can be recovered by GSHP system as shown in Fig. 2. The stored heat is also extracted from the ground by using U-tube GHEs. In the heat recovery process, the initial temperature of the ground surrounding U-tube GHEs is not homogeneous because of the thermal interactions during the storage process. In this case, more heat is expected to be extracted from the ground as the temperature difference between U-tube GHEs and ground is large. Therefore, more heat can be delivered to the building. As shown in Figs. 1 and 2, control system governs the temperature of the circulated fluid and their volumetric flow rates. It optimizes and governs the inlet temperature of heat pump and prevents its inlet temperature exceeding 40°C (due to the protection of heat pump) by passing the extra heat to the inlet of GHE again (through line × in Fig. 2) [29,35].

Numerical modelling approach
Thermal performance of U-tube GHEs and efficiency of GSHP is greatly influenced by the distance between GHEs [30], shank space, borehole vertical length, major diameter of borehole [36]. All these parameters should be properly optimized to achieve maximum efficiency. However, dynamic modelling of this storage system using computational software such as COMSOL and ANSYS-FLUENT [10,[28][29][30]35] or analytical methods [31,32,37] is a difficult and timeconsuming task. To resolve this issue, a novel numerical modelling framework using EES is developed in this section. The technique can solve complicated problems of GHEs modellings in short time (less than 20 s using single CPU). The numerical modelling consists of two main parts of modelling of the ground and fluid flow inside the U-tube. After separate investigations of each part, they will be coupled with each other to form a unified framework.

Implicit modelling of ground
Time dependent boundary conditions are a frequently encountered problem in transient modelling of thermal systems. The transient boundary conditions and the complicated geometry of GHEs hinder investigation of GSHP system performance by mathematical modelling except for the high-fidelity numerical techniques [31,32,37]. An implicit method is employed in this study due to some of its advantages. In explicit methods, the stability of the calculations is governed by the selection of x and t, however no such restriction is imposed on the solution of equations in implicit methods. This means that larger time increments can be selected to speed up the calculations. Although in implicit methods the number of iterations is generally large, problem can be still solved very fast with no restriction [38].
Inside the U-tube GHEs, heat is mainly transferred by convection and conduction. As shown in Fig. 3, the borehole and ground are divided into nodes in i and j directions.
In point 1 heat is transferred by convection since there is fluid flow inside pipes, and in later points (2, 3, 4 …) heat is transferred by conduction by assuming that the borehole is situated above the water table. Transient behavior of each point is investigated by numerical solutions of one-dimensional unsteady condition and convection problems. 1 denotes the average value of inlet and outlet fluid temperatures ( = where T i is assumed to be constant value [10,35]. Temperatures of points 2, 3, 4 … (m) which vary as time passes are evaluated by using the following implicit equation: where m and n are horizontal and vertical numbers of nodes, α is thermal diffusivity of the ground, t is time and x Δ and y Δ are horizontal and vertical distances between nodes (shown in Fig. 3).
In the first part of the current work, initial temperature of all nodes is assumed to be the same is the ground temperature (T g ). That is Eq. (4) shows that, temperature of any location in the ground (T i j can also be affected by the ground surface temperature (T i,0 1 ) which is usually the annual average temperature of the storage site.
Borehole wall temperature (T W ) is another significant parameter, which needs to be considered in the numerical approach. T W varies in different time steps and wall temperature of each segment (T W,j ) is assumed to be the average value of all horizontal nodes of that segment. Expectedly, increasing the number of nodes leads to more accurate model. T W is determined by: where n is the number of nodes. In the current study 240 nodes with 0.1 m distance are employed, based on a previously published work [28,30]. Further, diffusivity of heat is decreased significantly in far distance from the borehole wall (specially more than 6 m) [28,30].

Empirical modelling of the thermo-hydraulics of the flow inside pipes
Experimental results of the fluid flow and heat transfer inside the Utube GHEs are usually expressed in the form of empirical correlations. A large number of empirical relations for pipe and tube flows under different flow configurations can be found in the literature [38]. In this study, an accurate empirical formula, expressed by Eq. (7) [38][39][40], is used.
where f is friction factor and is calculated by; and Re d is Reynold number given by: where k f is thermal conductivity of fluid and d is the internal radius of the pipe.

Calculation of heat transfer rate in the U-tube
Ground implicit model and empirical model of fluid flow inside the pipes should be coupled to investigate dynamic thermal behavior of GHE and calculate heat transfer rate. There are two ways to calculate the heat transfer rate inside the U-tube. Firstly, since there is no phase change in the fluid, the heat transfer rate can be determined by the following thermodynamic equation; Secondly, it can also be calculated by accounting for convective heat transfer through the following equation; and T w,j is the borehole wall temperature in each j segment and is calculated by Eq. (6). Total heat transfer rate of U-tube GHE is sum of heat transfer rates of each j segment (Eq. (15)).
Heat transfer rate per unit vertical length of GHEs is defined by; where L [m] is the vertical length of U-tube GHE.
As shown in Fig. 4, to increase the accuracy of numerical solution Utube GHE is discretized into the 2n segments such that the outlet of each segment is the inlet of the next one. This method highly increases the precision of Eqs. (7) and (12).
In order to investigate the thermal behavior of GSHP system, all discussed Eqs. (4)- (15) need to be solved simultaneously in the same network. They are coupled with each other using Engineering Equation Software (EES) which has thermo-physical databank for different types of working fluid and is one of the most useful environments for solving thermodynamic and heat transfer problems. Framework developed in EES for this study consists of more than 9000 variables and equations which are solved simultaneously. Table 1 summarizes the boundary and working conditions and the numerical modelling approach.

Validation of the newly developed numerical model
The accuracy of the developed framework is investigated by comparing against the results from the literature [31,36].

Model validation, Case 1
Long term performance of a borehole with a single 1U tubes has been investigated by Aydin and Sisman (2015) [36]. Physical properties and operating conditions reported by them are given in Table 2.
By applying the same operating conditions presented in Table 2,  • It is the average initial temperature of the application area.
• It will be changed after GSHP start working.
• It is equal to ground initial temperature (Tg) • Temperature distribution is uniform.
Real temperature of various points in the ground (T i j , 1 ) • It shows real temperature of different points in the ground (i,j) shown in Fig. 3.
• It is variable and will be changed after start of the procedure.
• It should be calculated and updated by using Eq. (4) in each time step.
Ground surface temperature (Tsurface=T i,0 1 ) • It is the average of weather temperature in installation site.
• It is assumed to be constant.
• Referring to Fig. 3, T i,0 1 is constant and equal to T surface.
Borehole wall temperature (Tw j , ) • It is the average wall temperature of the borehole in each row (j).
• It is variable and will be changed after start of the procedure.
• It is assumed to be equal to the average value of real temperature of the ground (T i j • Eq. (7).
Total heat transfer rate of the borehole (Q t ) • It can be simply calculated by inlet and outlet fluid temperature differences (Eq. (12)) • Inlet fluid temperature (T inlet ) is assumed to be constant • Outlet fluid temperature is variable and has to be calculated at each time step • It should be calculated by Eq. (13) • To calculate heat transfer coefficient of the fluid (h f ), Eqs. (9)-(11) should be solved together.
• To increase the accuracy of the calculations, U-tube GHE is divided into the segments (Fig. 4) Temperature of the fluid in each segment (T j ) • It is variable and changed during the operation • It is shown in Fig. 4  long-term performance of a single 1U borehole GHE is evaluated using the newly developed equations network framework. Fig. 5 depicts that the results of the numerical model are in excellent agreement with those of Aydin and Sisman [36].

Model validation, Case 2 and 3
The experimental and numerical results presented in Ref. [31] are also used to validate the developed framework. Dehghan B. and Kukrer [31,32], developed 1D and 2D analytical solutions based on Green's function. These works showed that although Green's function based methods are useful, they are complicated and time-consuming. The operating conditions and physical properties are as shown in Table 3 and the results of validation are shown in Fig. 6. Fig. 6 shows that the results are in good agreement with each other. The compared results clearly demonstrate that the developed numerical framework can be an accurate, reliable and fast alternative to the computational and complex analytical solutions. Thermal performance of GSHP system can be determined as well as temperature of any point inside the ground, while the outlet fluid temperature of GHE can be evaluated at different time steps. Furthermore, this framework has the possibility of adding more GHEs and investigating the effects of various parameters such as number of GHEs and the distance between them upon the efficiency of GSHP system. It also gives the option of extending the equations network to 2D and 3D analyses, which highly reduces the simulation time.

Comparison between different approaches
To further demonstrate the relative merits of the developed numerical approach, numerical results of case 2 are compared with the ones obtained by Green's function method [31,32,37], g function approach [37] and COMSOL environment [10,[28][29][30]35] under the same working conditions. As already discussed, the developed numerical method is simple, reliable and fast. There is no need for high speed computers and it can be used in open sources software such as python. However, the complexity of 1D and 2D Green's function and g-function approaches can be seen in Ref. [31,32] and the work of Man et al. [37]. COMSOL [10,28,29,35], TRNSYS and ANSYS Fluent are useful, however, their processing time is long and they are often expensive. Fig. 7 shows the thermal performance of a single U-tube GHE (case 2) calculated by two different methods. Also, the processing time of each method is shown and R 2 (goodness of fit) values are calculated by referring to experimentally validated numerical results [31]. As shown, the developed numerical approach is the faster with an acceptable R 2 value (0.985). It is important to note that, to avoid any deviation from the main objectives of this research, only results of different approaches have been given and more details can be found in the previously mentioned Refs. [28][29][30][31][32][35][36][37] 4. Numerical modelling of biomass waste heat storage process

Calculation of the biomass combustion plant
Here, the biomass fuel is Pine pellets for which the physical and chemical properties are listed in Table 4. Based on the ultimate analysis (dry ash free based), the chemical formula of the biomass fuel can be expressed as CH O 0.89 0.94 . Therefore, the stoichiometric combustion reaction is given by: In practice combustion systems operate with excess air to ensure complete combustion and avoid formation of pollutants. The excess air coefficient is x air , the chemical reaction becomes: Considering the operational parameters of the biomass boiler, taking that the excess air coefficient is 1.4, the boiler thermal efficiency is 0.65, and the flue gas heat loss percentage is 0.29 [41,42]. Therefore, the temperature and mass flow rate of flue gas were obtained under different boiler output powers by applying Eqs. (1)-(3). The results are summarized in Table 5.
In the rest of this study, the boiler with 1000 kWt output is considered.

Waste heat storage process and analysis
Section 4.1 implied that in a biomass combustion system large amount of high temperature heat is wasted into the environment. This wasted high quality thermal energy can be used efficiently by supplying heating demands of buildings. The performance of a GSHP system can be significantly increased through waste heat recovery from the biomass thermal plant. In the considered biomass combustion plant, heat is wasted into the environment at 485.9 K and flow rate of 0.773 kg/s and the goal is to cool the flue gas to 300 K. In this study, water is chosen as the coolant fluid and as shown in Fig. 1, storage temperature should not exceed 90°C. To ensure about this, controllers have been implemented.
Operating conditions as well as boundary conditions given in Table 6 are applied to simulate storage process in the ground [10]. For the specific case study and based on the calculations presented in 4.1, to be able to cool exit waste heat temperature down to 300 K, at least 4 units of GHEs are needed (n = 4) which should be placed 7 m apart from each other. Different case studies with various storage temperatures, flow rates, etc. can be investigated by using the framework presented in this work. Fig. 8 indicates multiple U-tube GHEs which are used for heat injection and extraction in this work. During the storage process, fluid (water) is pumped into the ground at a constant temperature 90°C with constant flow rate of 0.355 kg/s [10,28]. Temperature of the ground around GHEs is homogeneous and the average value is 18°C. By pumping water at 60°C into the ground, heat is transferred by conduction from the borehole wall into the soil (ground). Therefore, temperature of the soil around GHEs starts increasing. Investigations show that, the ground temperature is not highly affected by the distance between U-tubes when it exceeds 6 m during the storage process. Therefore, second, third and fourth etc. GHE can be placed with 6 m distance between them as shown in Fig. 8. When the number of GHEs is more than one, the related boundary conditions should also be revised in the equation network. Fig. 9 shows the maximum amount of heat that can be stored in the ground versus time (3 months storage) under three different storage temperatures (90°C, 50°C and 40°C). It is obvious that, a significant amount of heat can be stored and recovered. During 3 months non-stop storage process, average stored heat in the ground (q s ) is 83.15, 75.77 and 70.40 W/m when storage temperatures are 60°C, 50°C and 40°C, respectively. As expected, the amount heat storage decreases as the time passes due to the thermal interactions between the borehole and its surrounded soil. Temperature of the ground around the borehole increases which leads to a non-uniform temperature distribution around boreholes. Fig. 10 indicates time dependent temperature variations along line x-x (middle depth) between two adjacent boreholes corresponding to 60°C storage temperature. As expected, temperature of the soil around the boreholes increases as the time passes. Here, 60 nodes with 0.1 m distance between them are applied (distance between boreholes is 60 × 0.1 = 6 m). Decreasing distance between boreholes (d) leads to strong thermal interactions between them and this issue causes performance losses in the system [28,30]. Temperature variation between boreholes 2, 3 and 4 is the same as that shown in Fig. 10 due to the symmetric behavior of the application area.

Waste heat recovery through ground source heat pump
In this section, thermal capacity of the application area which consists of 4 U-tube GHEs after completing storage process discussed is investigated. Stored heat is recovered and delivered to the building to supply heating demands through GSHP system. For this purpose, evaporator of the heat pump should be heated by an external GHE loop and condenser connects to the building in heating mode. To extract stored heat, fluid (in this study, water) is pumped into the ground at approximately 1°C under the most critical operating conditions (3 months non-stop operation) and then extracted heat is delivered to the building by the heat pump. Under these conditions, temperature of the ground around the boreholes is not homogeneous which means Eq. (5) is not valid in that area. To be able to analyze thermal process of the ground after storage, the average initial temperature of each node shown in Fig. 8 around GHEs is defined to be equal to their last value at the end of the storage process (at the end of 2160 h operation), as shown Fig. 10. Table 7 presents the working conditions of the storage and extraction systems.
In the equations network the total extracted heat from ground is defined as; where n is the number of boreholes. Fig. 11 illustrates the total extracted heat transfer rate from the ground (q e ) by four borehole GHEs (shown in Fig. 8) per unit length of borehole (W/m) in heating mode (in heating mode it is assumed that average T inlet shown in Fig. 4 is about 1°C). For the current case study, the amount of heat extracted from the ground increases by 7.6%, 14.4% and 23.7% per unit length of the borehole when the storage temperature is 40°C, 50°C and 60°C, respectively. These values are evaluated   under the most critical working condition (3 months non-stop operation) and the real performance is better than the results given in Fig. 11 due to the intermittent operation of GSHP system. Amount of heat that can be delivered to the building is depended on the heating Coefficient of Performance (COP) of the heat pump and the vertical length of GHEs. Based on the buildings' need and maximum available budget, these parameters can be optimized. Temperature distributions in the ground are essential for calculations of heat losses of the buildings to the ground and design of relevant equipment [43]. Fig. 12 presents the temperature distribution in the ground along line x-x at the beginning and at the end of extraction operation. This figure shows that the initial temperature of the ground (t = 0) after the storage process is not homogeneous and it depends on different parameters such as storage temperature, distance between nodes (x), distance between boreholes (d), etc.
As discussed in the previous section and in Fig. 4, temperature of fluid flowing through the pipe is different in each segment. Fig. 13 illustrates fluid temperature variations inside the pipe in different segments (T j j = 1,n). The pipe has been divided into 40 segments and fluid temperature variations are shown versus three different operation times. As the time passes, temperature differences between each two segments declines due to the saturation of surrounded ground heat capacity. Accordingly, borehole wall temperature variations are shown in Fig. 14. Based on the results, borehole wall temperatures calculated by Eq. (6) are different at each segment.    Another important parameter which can affect the performance of GHEs, is the average ground surface temperature. To investigate the effects of ground surface temperature on the thermal performance of GHEs, the amount of possible extracted heat from the ground is calculated for the range of −10, 20°C average ambient temperature. Based on the results given in Fig. 15, thermal performance of GHEs will not be significantly affected by the ground surface temperature when L > 50 m. When vertical length of U-tube GHE is 100 m, the amount of extracted heat from the ground varies for less than 2% over −10, 20°C range of ambient air temperature. Results also show that, surface ground temperature can play an important role on the thermal performance of shallow GHEs such as helical and slinky.
For further clarification of the numerical approach, solution steps are listed below: (1) For storage process, at the beginning system is under equilibrium and temperature distribution of the application area is uniform (Eq. and 10) (6) For extraction and recovery of stored heat, temperature distribution in the ground is not uniform anymore and it has to be updated (Fig. 12).  16) and (19) should be calculated simultaneously again by applying new initial and boundary conditions. (9) Amount of extracted heat from ground (q e ) should be determined   (Eq. (19)) as well as heat that will be delivered to the building (Fig. 11).

Heat recovery performance and CO 2 reduction
Results of this study show that a considerable amount of waste heat can be recovered and beneficially used for meeting the heating demands of buildings. In Fig. 16 the amount of recoverable waste heat of different biomass combustion plant with different installed capacity (0.5-10 MW) versus two different numbers of U-tube GHE (n = 1 and 4) has been shown. Waste heat recovery ratio is the amount of recovered heat divided by the total wasted heat of the biomass plant.
In almost all combustion systems a fraction of heat generated by burning fuel is wasted by flue gases to the atmosphere. Burning fossil fuels results in emission of CO 2 and thus wasting thermal energy intensifies the emission of greenhouse gases. For those systems which burn biomass, waste of heat is essentially waste of renewable fuel and is therefore an environmental burden. The underground storage system introduced in this work offers an efficient way of storing a significant fraction of the heat that is normally wasted from the chimney. The preceding analyses showed that for the 1 MW combustion system under investigation between 30% and 100% of the waste heat can be successfully recovered and delivered to buildings for space heating purposes. This range can be made even wider by implementing larger number of boreholes. Assuming an average recovery rate of 65%, it can be readily shown that the recovered heat saves emission of almost 1ton of CO 2 per day, in comparison with the case of burning natural gas for supplying heat to the buildings. If the heat is to be supplied by biomass combustion approximately 700 kg of biomass should be burned per day. Clearly, higher storage capacities applied to bigger combustion plants will result in larger heat recovery and can further reduce the CO 2 emissions and save biomass. It is essential to note that the calculations throughout this work were conservative and improvements in heat exchangers efficiency and the specifications of boreholes can readily increase the recovery rate.

Conclusion
Storage of heat wasted by a biomass combustion plant and heat recovery by using ground source heat pump were investigated. Results showed that a considerable amount of waste heat can be recovered and delivered to buildings. A novel, fast, highly accurate and yet simple numerical modelling framework developed in EES environment could solve the complicated problems of GHEs modellings in a short time (in less than 20 s). Results showed that for 3 months continuous storage processes the average rate of heat storage in the ground (q s ) is 83.15, 75.77 and 70.40 W/m when the storage temperature is 60°C, 50°C and 40°C, respectively. For the considered biomass combustion plant, the amount of extracted heat from ground increases by 7.6%, 14.4% and 23.7% per unit length of the borehole corresponding to 40°C, 50°C and 60°C storage temperatures. Based on the results, it is recommended to   use 60°C as storage temperature for which at least 23.7% more heat can be extracted.

Declaration of Competing Interest
The authors declared that there is no conflict of interest.