Investigations into thermal resistance of tunnel lining heat exchangers

. Geothermal energy is a promising and sustainable source that can reduce current dependence on conventional fuels for thermal energy production. To exploit this source of energy thermo-active geostructures such as tunnel lining heat exchangers are being investigated theoretically as well as experimentally. These geostructures are composed of concrete panels embedded with reinforcement cages fitted with absorber pipes. Several engineering projects in China, Finland and Italy have deployed such heat exchangers in tunnels. To achieve efficient energy production, characterisation of these systems require realistic models of the substructure heat exchanger. Therefore investigations into thermal resistance of the heat exchanger is vital. The present study is concerned with quantifying the thermal resistance of tunnel lining heat exchangers where the thermal boundary surfaces are applied at surfaces representing the adjacent ground and the exposed concrete, in addition to the pipe surface. Steady state temperature distribution in a two dimensional cross section of a tunnel lining heat exchanger is investigated using the boundary collocation least squares method. Design parameters including pipe and tunnel lining specifications are used as model inputs.


Introduction
As fossil fuels and other natural resources continue to decline in availability, governments around the world are looking towards sustainable energy as a viable alternative. In the UK the government has committed to reaching net zero emissions by 2050 [1]. Correspondingly, the proportion of fossil fuel use in the UK has been falling throughout the last decade. In 2019, energy from renewable generation exceeded that from fossil fuels for the first time [2]. However, the transformation to renewable and low carbon energy in the UK must continue if targets are to be met.
Over the same period the carbon intensity of UK grid electricity has halved [3]. Yet heating, one of the largest emissions sectors, remains predominantly supplied by direct burning of fossil fuels. Despite this need to rapidly decarbonise heating, consideration of geothermal energy in the UK remains relatively limited. This kind of energy can provide heating, cooling and power generation from hydrothermal resources while power outputs can remain stable and unaffected by changes in the climate [4]. This provides a situation where geothermal energy could be more effectively utilised in order to reduce CO2 emission levels.
In more recent years, the application of heat exchangers has turned towards concrete substructures; including slab foundations, bored piles and diaphragm walls; to absorb heat from the ground by embedding heat exchanger pipes within the structures themselves. Although numerous studies are carried out to investigate the performance of substructure heat exchangers less attention has been given to tunnel lining heat exchangers [5]. To that end, this research is focused on deriving a semi-analytical solution for calculating thermal resistance of tunnel lining heat exchangers.

Tunnel lining heat exchangers
There have been several case studies in the literature where tunnel heat exchanger technology is being tested and implemented. Zhang et al. [6] studied the Linchang Tunnel in China to derive an analytical solution for the heat conduction as well as investigating the thermal performance of the tunnel lining during operation. One of the predominant reasons for undergoing this study was to mitigate against damage to the tunnel caused by freezing by extracting heat from the surrounding bedrock. For economic reasons, they found that the optimal flow velocity of the carrier liquid should be 0.8 m/s to account for cost and thermal performance. In addition, the absorber pipes' (i.e. pipes that extract heat from the surrounding bedrock) geothermal energy output linearly decreased as the inlet temperature of the carrier liquid increased.
Zhang et al. [7] also investigated the aforementioned factors with the addition of the distance of the pipes. The dimensions of the pipes were an outer diameter of 25 mm and a thickness of 2.3 mm. It was concluded that the distance of the pipes, plus the flow rate and inlet temperature, all had a significant impact on heat exchange rate of the pipes themselves. Nicholson et al. [8] deals with the design of a tunnel lining GHE (referred to by the authors as a TES system) for Crossrail, UK. This paper is good at identifying where potential excess heat originates from, including the trains' motors, brakes and air conditioning units. It was proposed to use heat generated in the tunnels to heat buildings that lie within 100 m of the route alignment which varied in size and use. Wiik [9] presented preliminary calculations and similar conclusions when discussing geothermal energy utilisation within a city railway loop in Helsinki, Finland. Calderón et al. [10] studied the design, numerical modelling, and thermal behaviour of a flat heat exchanger integrated to the concrete wall of a subway tunnel in Chile and reported the dependency of heat absorption on the heat carrier flow rate, and on the spacing of pipe arrangements. They considered actual environmental conditions in the tunnel and the surrounding ground properties for their numerical analysis. Barla and Di Donna [11] provided guidance on how to proceed in developing thermal and mechanical designs of energy tunnels through effective use of computational methods. In a more recent paper, Insana and Barla [12] published their investigations on the experimental and numerical analysis of the energy performance of a real-scale energy tunnel prototype in Torino, Italy, aiming at providing quick and effective tools to quantify heat exchange in such substructures in the preliminary phase of the project. Although complex numerical simulations are an effective tool in estimating thermal performance of such heat exchangers, they are indeed costly and require careful preparation and execution. Therefore, making use of less sophisticated methods which can give good accuracy in a less time demanding manner would be greatly beneficial.
In the present manuscript we use a semi-analytical method to derive the thermal resistance of tunnel lining heat exchangers using a given boundary condition and key geometric parameters of the energy tunnel.

Thermal Resistances
In a substructure heat exchanger, thermal resistance between the fluid in the embedded pipes and the heat exchanger wall is a vital performance parameter. Thermal resistances can be used in combination with other solutions for the thermal transfer in the surrounding ground to predict the outlet temperature of a specified substructure heat exchanger or to predict the heat flow rates. This approach may be less sophisticated and contain more simplifications than the composite medium multiple line source approach of Zhang et al [6], but it has the advantage of being much simpler to adopt and hence potentially valuable to routine application.
A practical method to obtain values of thermal resistances for substructure heat exchangers is to calculate the shape factors for the two dimensional representation of their geometry in the direction of the heat flow. Shape factors, S, are a lumped parameter that depend only on the complex geometry of the heat exchanger and its thermal boundary conditions. They allow the effects of these geometric components to be separated from the thermal conductivity which is the other factor that controls the thermal resistance: Thus, by applying the appropriate boundary conditions to the two dimensional geometry of the substructure heat exchanger and the steady state heat equation it is possible to calculate thermal resistances. Shape factors are mainly dependant on the heat exchanger's geometry, thus, variations in their geometric parameters can greatly impact the values of shape factors. The initial step to derive shape factors for a given geometry is to find the steady state temperature distribution in that geometry. The approach to derive the shape factors for a tunnel lining heat exchanger is provided in the next section.

Problem definition
In this research, shape factors for a tunnel lining equipped with one set of heat exchanger pipes are investigated where all boundaries are kept at constant temperatures. The pipe is exposed to a temperature of Ti while the tunnel has a temperature of To at both extrados and intrados.
While it is accepted that there will in fact be a gradient of temperature across the lining, these boundary conditions are selected for simplicity reasons only in this initial study which will test the concept of resistance applied in this context. This will also permit a first assessment of the sensitivity of different geometrical factors on the thermal performance of the tunnel lining. A cross sectional view of the studied geometry comprising a set of equally spaced pipes in a tunnel lining is shown in Figure 1. Tunnel lining has a thickness of ( − ) and the pipes each has radius ri, a separation distance of Ps and an eccentricity of e ( Figure 2). Considering the symmetrical configuration of the geometry it can be reduced to that existing about one pipe as shown in Figure 2. The problem is then reduced to a plane region bounded by an annular inner boundary and an outer boundary where as a result part of the outer boundary becomes adiabatic. Adiabatic conditions on these sides of the lining panel reflect the fact that the heat flow is predominantly radial in terms of the overall tunnel centre. However, it is accepted that there may in fact be small contact resistances at these boundaries that we are therefore neglecting in the first instance. Ultimately, the domain can be further reduced to that in Figure 3 for its symmetry along the central horizontal line.  To derive the solution for shape factors it is crucial to find the temperature distribution in the domain illustrated in Figure 3. If a relation can be found for temperature distribution in this region, the shape factor can be calculated from integrating the normal gradient of the temperature over the inner boundary, i.e. the pipe surface, as shown in Equation (2): where Φ is the dimensionless temperature. The equation that governs steady state temperature distribution in the domain illustrated in Figure 2 and Figure 3 is the Laplace's equation in plane polar coordinates: To simplify the analysis, the following dimensionless parameters are substituted in the governing equation: where radius of the pipe, ri, is used as the characterised body length. Thus, the dimensionless formulation condenses the boundary value problem into this form: A linear arrangement of the trial functions that satisfies Equation (5) can be: ( , ) = + ln + + ln where , , , , , , , are unknown constants and to find them applying the thermal boundary conditions in Equation (6) is necessary. Applying the boundary conditions along the line where = 0 and = , i.e.: = 0 ( 7 ) = 0 will set , , and constants in Equation (6) to zero and to n=1,2,3, … . Hence, Equation (6) is reduces to: ( , ) = + ln To simplify Equation (9) boundary conditions imposed on the other surfaces need to be enforced. Using the boundary condition on the pipe surface, i.e.: Equation (11) gives an expression that is general enough for temperature distribution in the domain. Using Equation (11) in Equation (2) it can be seen that the values of the shape factors are only dependent on the value of and equals −2 . To derive these values the boundary conditions on the outer surface must be employed. The boundary condition for the right and left arcs of the outer boundary is (Figure 2 and Figure 3): where ri and rs are the distances from the pipe centre to the pipe surface, and the outer surface of the geometry respectively, as shown in Figure 2. While part of the outer boundary between the two arcs is exposed to adiabatic conditions and therefore the normal gradient of the temperature along this section of the boundary is equal to zero. Using the dimensionless variables and employing some effort the boundary condition on this part of the outer surface can take the form shown in Equation (13): The specified constant temperature boundary conditions are then used to construct a set of linear equations for solving the Fourier coefficients by attributing certain values to the angle θ. The coordinates of the collocation points are given by a combination of the angle θ and the distance = . Due to the complex configuration of the domain it is not possible to derive an analytical solution for the expression of shape factors and thus thermal resistances. In these conditions making use of less sophisticated numerical approaches can provide results with good accuracy. To that end we use a semianalytical boundary collocation least squares method, details of which are provided in the next section. If alternative boundary conditions are applied at the inner and outer surface of the tunnel, e.g. convective conditions, Equation (13) would still be valid.

Boundary collocation method
The idea behind a boundary collocation method is that the approximate solutions, to agree with the applied conditions, are imposed along the boundaries of the domain in the form of collocation points. By using a certain number of collocation nodes, M, certain values of the angle θ can be generated and the radius is expressed as a function of that angle. Subsequently, the infinite series in Equation (11) can be truncated to its first N terms and a system of linear algebraic equations are obtained for the unknowns B and An.
Frazer et. al [13] and Bickley [14] are the first to use boundary collocation method along with least squares technique and the Galerkin method to solve the unsteady heat conduction problems. Using the least squares method minimizes the residual and increases the precision of the solutions obtained. In using a least squares approach an overdetermined system of linear equations can be used so that the number of collocation nodes can exceed the number of unknowns. Details of employing boundary collocation least squares methods can be found at [15].
To avoid large errors occurring in off-point intervals of the boundary, refinement of the point spacing is required. Thus, the collocation points are not set arbitrarily but are uniformly spaced on each boundary element and the angles θ are obtained accordingly. This procedure is likely to produce a better approximation due to the even distribution of the collocation nodes. To obtain the coordinates of the selected collocation points on the outer boundary we have divided it into three regions based on the location of the eccentricity.

Results and discussions
Using the mathematical formulations introduced in the previous section and applying the corresponding boundary conditions it is possible to derive values for the tunnel lining shape factors with varying ratios of geometric parameters. In our calculations we have assumed that the pipe is located closer to the extrados (the surface of the lining which is in contact with the soil) of the tunnel lining. Table 1 details variations of the shape factors with increasing tunnel lining diameter while pipe diameter, pipe spacing and cover (distance between the outside of the pipe and the ground) values are kept constant. The pipe diameter is kept at 25 mm while the cover was set to 85 mm and the angle of pipe spacing to 5º.
Data in Table 1 demonstrates that increasing the tunnel lining diameter while keeping its thickness constant reduces the shape factor values and thus will increase their thermal resistances. Variations of shape factors with pipe diameters are shown in Table 2. These data are for a fixed ratio of Ro/Ri equal to 1.06 and tunnel diameter of 8.5 m, while pipe spacing is kept at a constant value of 85 mm. Shape factors are increasing with increasing pipe diameter for a fixed outer and inner tunnel diameter resulting in reduction of corresponding thermal resistances. However, increasing the pipe size needs to be within reasonable boundaries. It can be explained for by the fact that increasing the pipe diameter for a fixed outer boundary dimension will decrease the area of the doubly connected domain, as shown in Figure 2, and thus the corresponding thermal resistances. For a given tunnel lining with fixed inner and outer diameters, pipe diameter and pipe spacing increasing the cover will reduce the shape factor as shown in Table 3.
Variations of the shape factors with one varying geometric parameter, while the rest are kept constant, demonstrate a linear relationship.
Generally, thermal conductivity of concrete varies between 0.8 to 2.5 W/m 2 K, depending on its composition. These values can be used along with the introduced shape factors to estimate the corresponding thermal resistance of the tunnel lining heat exchangers. This suggest thermal resistance values in the range 0.1 mK/W for larger pipes and higher concrete conductivity, to 0.5 mK/W for smaller pipes and lower concrete conductivity.

Conclusions
An implementation of the least squares boundary collocation method to evaluate the conduction shape factors and thus thermal resistances of tunnel lining heat exchanger, with isothermal inner and outer boundary conditions, is reported. The data shows that variations in key design parameters including pipe diameter, tunnel thickness, tunnel internal and external diameters, cover, and concrete conductivity can affect thermal performance of such heat exchangers. Further research can be performed to investigate varying boundary conditions. Employing more sophisticated numerical analysis such as boundary integral methods may also improve the accuracy of the results.