Numerical modelling of different applications in Energy Foundation Technology

. Energy foundation technology is expected to make a significant contribution to the use of renewable energy. In this context, this paper presents the use of Finite Element simulation using PLAXIS software for modelling different benchmark applications in Geothermal Foundations. An implicit fully-coupled numerical scheme with global adaptive time stepping are implemented to ensure computational efficiency and stability. Firstly, a transient simulation of thermal response tests [1], often used to estimate the thermal conductivity of ground and thermal resistance of pile, is presented. In the second part, a Thermo-Hydro-Mechanical analysis is performed to simulate the behavior of a single heat pile subject to a thermal load cycle [2]. Several ingredients (constitutive behavior, interface finite elements etc.) are employed to simulate soil-structure interactions. The obtained solutions are validated against available simulation and experimental data to demonstrate the applicability of the simulator to energy foundation analysis and design.


Introduction
Together with ground source heat pump (GSHP) and groundwater heat pump (GWHP) systems, recently underground geotechnical structures, such as deep and shallow foundations or tunnels are also being employed as energy geostructures by installing the absorber pipes directly in the structures [3]. In particular, energy foundation technology is expected to make a significant contribution to the use of renewable energy [4]. The design of heat exchanger piles requires numerical modelling in different steps, from the simulation of in situ testing to the analysis of mechanical responses and energy efficiency. As this technology becomes a vital and popular part of geotechnical engineering, the use of userfriendly computer codes would be needed to perform advanced analyses in a robust and accurate way [5] [6]. In this context, this paper presents the use of PLAXIS, a Finite Element software, for modelling different benchmark applications in Geothermal geostructures. An implicit fully-coupled numerical scheme with global adaptive time stepping is implemented to ensure computational efficiency and stability. Firstly, simulations of Thermal Response Tests [1], often used to estimate the thermal conductivity of ground and thermal resistance of pile, are presented. In the second part, a Thermo-Hydro-Mechanical analysis is performed to simulate the behavior of a single heat pile subjected to a thermal load cycle [2]. Different techniques (constitutive behavior, interface element…) are utilized to simulate soil-structure interactions. The obtained solutions are validated against available simulation and experimental data to demonstrate the applicability of the numerical approach and the capability of the simulator for energy geostructure analysis and design.

Finite element framework
PLAXIS generally considers a geomaterial as a 3-phases (solid/liquid/gas) porous medium, but gas pressure is assumed constant, which is widely adopted in geotechnical engineering. In this section, only the principal equations are briefly presented. For more details, readers may refer to [7] [8].

Balance equations
The mass continuity equation of water is written as follows: where (-) is the porosity and (-) the water saturation.
(-) is the volumetric strain of the soil skeleton and , and are the solid, water and gas densities. • denotes the divergence operator, the gradient while ⁄ the time derivative. and are the mass fluxes of liquid water and vapor, respectively.
The linear momentum balance of the whole porous medium is given by: where is the total (Cauchy) stress tensor, the gravity force vector and ρ the density of a multiphase medium, defined as: The energy balance equation for the porous medium can be written as follows: where , and are the internal energy in the water, vapour and solid phases, respectively. and are the advective internal energy flux in water and the conductive (diffusive) heat flux in the porous medium, respectively.
is the heat source term, i.e. heat generation rate per unit volume.

Basic constitutive equations
The stress and flux terms in the above balance equations can be specified by the following basic constitutive laws. The Bishop's effective stress [9] is used: where ′ is the effective stress tensor, the pore water pressure, the second order identity tensor and = the effective saturation with and being the saturated and residual saturation degrees. This effective stress applies for both unsaturated and fully saturated ( = 1) conditions. An extended Richard's model is applied to describe non isothermal unsaturated flow. The mass flux of water is defined as [10]: where k (-) is the relative permeability of the medium, the dynamic viscosity, κ is the intrinsic permeability tensor of the porous medium. In parallel, the mass flux of vapor in non-isothermal processes is defined as: where T is the local equilibrium temperature of the porous medium, is the vapour diffusion coefficient and and are the hydraulic and thermal diffusion coefficients, respectively.
The conductive heat flow (Eq. (4)) is assumed to be governed by Fourier's law: where is the thermal conductivity of the multiphase medium calculated by: where , , being the solid, water and vapor conductivities, respectively. The advective heat flux in (4), which represents the heat transported by water flow, is given by: where is the water heat capacity.

Numerical implementation
The above equations are discretized using standard Bubnov-Galerkin Finite Element method (FEM). Apart from soil triangular elements (either 6-noded (quadratic) or 15-noded (fourth order)), special elements, such as interface and plate, are also implemented to simulate realistically soil-structure interactions. A fully coupled scheme was implemented so that all primary field variables (displacements, temperature and water pore pressure) are solved within the same time step. This is physically realistic as it does not need any further "decoupling" assumptions within each step.
An implicit scheme is used for time discretization. The time marching scheme is automatic as the time step is adjusted according to convergence rate and accuracy. This relaxes step size restriction required in explicit discretization schemes to ensure stability.
Beside built-in geomaterial models, advanced models can also be defined using the well-known API of User-Defined Soil Model (UDSM) and User-Defined Flow Model (UDFM) module. Interested readers may refer to [8] for more detail. In the following, the numerical approach presented above will apply to different applications in Energy Foundation technology.

Simulation of Thermal Response Test
One of the key parameters in the design of ground sourced heating and cooling systems is the heat exchange capacity of the soil formations. For borehole heat exchanger systems, it is common to carry out an in situ Thermal Response Test (TRT) to determine both thermal conductivity of the ground and thermal resistance of the borehole, which are the key factors for the design of the system and will allow appropriate sizing of the borehole field to meet specific energy demands. As energy piles are used more and more frequently, it is becoming important to develop testing procedures and analysis methods to evaluate the thermal conductivity of such systems which are geometrically different than boreholes. For this purpose, Loveridge [1] used numerical modelling and case study data to evaluate the applicability of the standard TRT to pile heat exchangers. In the following, we reproduced the numerical tests performed by Loveridge [1] to demonstrate the applicability of our simulator to TRT simulation.

Model definition
Several 2D numerical models were built for different pile heat exchangers geometries, as listed in Table 1. The short nature of the piles allows to assess in the short term (timescale of thermal response test) no significant variation of the undisturbed ground temperature over the depth of the pile [1]. Thus, simple 2D models (a horizontal cross-section of the pile) are used. They include a concrete pile with a concrete cover c from the edge, plastic pipes and surrounding soil (Figure 1).
Different scenarios corresponding to different pile diameter as well as different number of pipes are tested (Table 1). For each scenario, we performed 3 different simulations corresponding to 3 values of the ratio between the thermal conductivity of soil and concrete (respectively = 1 , 2 0.5). Similar to [1], the pipes and circulating fluid are not modelled, instead thermal boundary conditions are applied on each pipe. Furthermore, the surrounding ground in the model is extended far enough to avoid any boundary effects. All the external boundaries are set to be fully insulated, while constant flux is applied at all the pipe edges. The finite element mesh is refined around those edges to enhance numerical accuracy ( Figure 1).

Fig. 1. Model geometry and FE mesh
An initial temperature of 11°C is assumed for both the pile and soil. The main purpose of this study is to produce the temperature response function for the ground at pipesoil interface. Therefore, transient thermal analyses are performed and both mechanical and hydraulic behavior are neglected. Material properties are listed in Table 2.

Simulation results
The purpose of the analysis is to compare the obtained results with existing numerical data from [1] as well as analytical methods generally used in heat exchanger design.
In the analytical solutions available in the literature, a dimensionless temperature Φ = 2 ∆ / , where Tg is the ground temperature and q the heat flux, is often expressed as a function of Fourier Number = /( ), where t is the time, c the soil heat capacity, the soil thermal conductivity and rb the pile radius. For instance, assuming a constant heat injection , the temperature change in the ground computed at a radial distance = can be described by the following closed form expressions: • Line source heat model [11] ∆ = ln − where = 0.5772 is Euler's constant • Cylindrical source heat model [12]: These three solutions are compared to the numerical results obtained. Figure 2 shows a typical temperature distribution for the case where pile diameter is 600 mm with 4 central pipes. Note that to compare with the analytical solutions, the value of ∆ from the numerical results is averaged from non-uniform temperature along pipe circumference). Some typical numerical/analytical comparisons are presented in Figure 3 for the scenarios 1, 3, and 6 (see Table 1). A close agreement between PLAXIS results and the numerical data reported in [1] is observed. As expected, the analytical and numerical solutions follow the same trends but do not coincide. The difference may come from the simplifying assumptions of the analytical approach, where the pipe geometries and positions are not considered.
As aforementioned, design methods that consider a line source may be valid for boreholes heat exchangers

Simulation of THM response of a single heat pile
When piles are used as heat exchanger, additional loads may result from heating-cooling cycles. Generated excess pore pressures, stresses, and deformation should be considered in the design process. Finite element THM simulation is believed to be a useful and reliable tool to assess correctly this information [10]. In this section, an in-situ experiment for the behaviour of a single heat pile [11] will be simulated to illustrate the capability of our numerical framework.

Model definition
The numerical model consists of an axisymetric representation of the model with a concrete pile (length of 26m and diameter of 1m) embedded in a stratified ground. The ground is extended far enough to ensure an imperturbation at far-field. The corresponding mesh is shown in Figure 4.  Plastic pipes and circulating fluid are not modelled and only a thermal boundary condition is applied to the pile-soil interface. The initial, as well as the far-field (unaffected) temperatures are fixed at 11°C. Radial displacements and flow are fixed at the symmetry axis, while the head of the pile is not restrained because no mechanical loads are considered in this test. According to [2], a constant temperature of 11° is applied to the top surface to simulate the atmosphere temperature and a E3S Web of Conferences 205, 06004 (2020) ICEGT 2020 https://doi.org/10.1051/e3sconf/202020506004 heating-cooling cycle ( Figure 5) is applied to the pile shaft boundary to simulate heat exchange by the pile, according to the experimental conditions [11]. The properties of the soil layers are based on the data reported in [2] and [14] (Table 3a,b). The concrete is supposed to behave elastically (LE) while the soil has an elasto-plastic Mohr Coulomb behavior (MC) (see [7] for more detail).
To simulate realistic soil -structure interaction, interface elements are used at the pile -soil contact. In each geological layer, interfaces behave similarly to soil material of that layer, but with a strength reduction ratio Rinter. As the concrete pile is assumed to be non-porous, these interfaces are set to be impermeable for water flow. For thermal behavior, we assume that the interfaces are infinitely conductive (without any thermal resistance). The interface element parameters are also presented in Table 3a,b.
"Staged construction" concept is used to model this in situ test with different construction phases, including a pile installation phase and thermal loading phases. The pile installation phase is particularly useful to ensure an equilibrium stress state prior to the thermal phases. Although heating and cooling can be applied easily within one unique phase, we decided to split heating and cooling processes into 2 phases to capture precisely the time evolution of each period the test. In the following, we will only present the results for the thermal loading phases where experimental measurements are available for numerical comparisons.  Figure 6 shows the temperature distribution at the end of the cooling period. As expected, hot temperature concentrates around the pile shaft due to the imposed boundary conditions, while at far-field the soil is not thermally affected (T=11°C). Because of the delay effect of heat conduction, the hottest zone, however, is not at the pile edge where the temperature is applied, but somewhere near that location. Figure 7 shows the temporal vertical pile displacement. A good agreement is observed between the obtained numerical results and the experimental ones measured by different techniques reported in [2], as well as with other numerical results reported in [15]. As expected, pile movement follows thermal loading path due to soil and concrete thermal expansion.

Simulation results
The imposed thermal field generates strains and displacement in the pile due to thermal expansion. The profiles of vertical strains along the pile at the end of the heating and cooling periods are shown in Figure 8. Once again, we observe a reatively good agreement between the numerical results and the experimental data reported in [2] [15]. The imposed thermal field generates strains which distribute differently for each soil layer

Conclusions
This paper illustrates the applicability of THM coupled numerical approach for modelling of energy foundation systems. Experimental and numerical data available in the literature are successfully reproduced using the PLAXIS software. This helps to enhance the confidence on modelling capability for complex multiphysics problems encountered in this topic, which is generally beyond classic design in geotechnics