A Finite Element Model for Investigating Unsteady-State Temperature Distribution and Thermomechanical Behavior of Underground Energy Piles

Abstract

The underground energy geostructure represented by the energy pile is one of the key paths for the cooperative development of underground space and geothermal energy. Because of its advantages of low cost, high efficiency and no extra occupation of underground space, it has become a feasible alternative to the borehole heat exchanger. The change in the temperature field of the energy pile and its surrounding ground not only affects the geological environment but also influences the thermomechanical performance and the durability of the structure. However, the temporal and spatial unsteady-state temperature distribution of piles and surrounding rock under typical intermittent and unbalanced thermal load conditions is still unclear. In this paper, a finite element model was applied to analyze the unsteady-state temperature distribution, and the thermomechanical behavior of the energy pile group was developed and verified. The temperature field distribution of pile and surrounding rock under typical intermittent working and unbalanced thermal load conditions were determined. Moreover, the thermomechanical behavior characteristics of the energy pile group were investigated. Finally, the influences of pile layout on the thermomechanical behavior of the energy pile group were identified by designing six different scenarios. The results indicate that under typical intermittent operation conditions, the temperature of the energy pile and surrounding ground near the heat exchange pipe varies periodically. For areas with unbalanced cooling and heating loads, long-term operation of energy piles leads to thermal accumulation, and the maximum temperature of energy piles occurs in the first daily cycle. In summer/winter working conditions, the increase/decrease in pile temperature induces axial compression/tensile stress. When the pile group is partially used as the energy pile, the non-energy pile acts as the “anchor pile”, and it generates the added tensile stress.

1. Introduction

In the coming decades, increased use of renewable energy is needed to help reduce carbon dioxide emissions [1]. Shallow geothermal energy is one of the renewable energy that can solve the problem of heat and cold energy consumption in buildings. At present, shallow geothermal is mainly extracted by a borehole heat exchanger (BHE), which requires additional drilling and construction acreage (high initial costs) [2]. An energy pile is a set of heat exchange pipes buried in the foundation, which is combined with the building structure. By using the relatively stable temperature underground, the circulating fluid in a closed loop (through heat exchanger pipes) can bring heat/cold energy from the underground to the building [1,3]. It is a dual-purpose foundation system that incorporates heat exchangers inside the pile and provides foundation support for the building, thus making it possible to use shallow geothermal while withstanding the load of the superstructure.

As a practical alternative to the BHE, energy pile has been used more and more globally since they were first adopted in the 1980s [3]. Although an energy pile may provide less heat/cold energy than BHE, it can cover the basic heating/cooling needs of a building and show many advantages over BHE [4,5]. The energy pile is buried in a pile foundation to avoid the high cost of drilling and backfilling. Moreover, in the case of scarce urban underground resources, the energy pile does not need additional space to install, saving valuable space resources, and has outstanding applicability in densely built-up cities. Furthermore, the energy pile heat exchange system is composed of concrete and steel with high thermal conductivity, which can increase the energy transfer between the system and the underground. Hence, compared with BHE, an energy pile can provide a higher heat transfer capacity [6,7]. Moreover, in the past few decades, considerable efforts have been focused on energy piles through full-scale tests [3,8], model tests [9,10] and numerical modeling [11,12].

On the one hand, the heat transfer performance of energy piles and their influencing factors have always been a concern. Bandos et al. [13] proposed a cylindrical source model to evaluate the heat transfer performance of energy piles, which can consider the heat capacity in the borehole. Park et al. [14,15] compared the heat transfer performance of different types of energy piles through in situ tests. Yang et al. [16] investigated the influencing factors of heat transfer performance of energy piles, such as fluid inlet temperature, operation mode, and structural materials. On the other hand, the thermomechanical behavior of energy piles has recently become a new research hotspot. Zhang et al. [17] studied the variation law of peak lateral resistance of energy pile in the process of external force loading/unloading at different temperatures. It was found that the heat and moisture migration occurred around the pile during the temperature increase, which led to the decrease in water content around the pile. Yavari et al. [18] conducted cyclic thermal response tests on energy piles and found that the heating/cooling process can induce reversible rise/settle of the pile top when the pile top load is less than 30% of the ultimate bearing capacity, respectively. Knellwolf et al. [19] found that the constraint of the superstructure on the pile top significantly affected the distribution of axial stress and the temperature-induced radial deformation of the energy pile, as well as the lateral friction resistance and the bearing performance of the pile foundation.

There are some limitations in the in situ experiment and analytical method, and the reliable numerical model has excellent versatility in the study of the energy pile [20]. Compared with other numerical methods, the finite element method can efficiently characterize transient effects; this ensures the specified engineering precision and can deal with complex conditions (such as different types of materials) [20]. Therefore, increasing attention has been paid to the study of the thermomechanical behavior of energy piles using the finite element method. Park et al. [21] combined the finite element model with the heat transfer theory, analyzed the heat transfer mechanism and heat conduction characteristics of energy piles, and subsequently studied the influence of different media on the thermal conductivity of energy piles. Ozudogru et al. [20] analyzed the mechanical performance of energy pile under the influence of different soil strength parameters; the results show that the axial elongation and contraction of energy pile under the constraint condition is unfavorable to the stability of pile foundation. Gashti et al. [22] investigated the temperature-induced stress and found that the dynamic lateral friction (induced by temperature) has little effect on the bearing capacity of end-bearing piles. Olgun et al. [23] investigated the reasons for the increase in bearing performance with temperature growth and found that the contact pressure between pile and soil is generally smaller than 15 kPa.

The unsteady-state temperature distribution of energy piles is the basis for the study of both thermomechanical behavior and heat transfer performance of energy piles [4,24,25]. This is because the generation of thermal stress during the heat transfer process may have a significant effect on the stress distribution of the pile foundation structure [1,22]. However, the temperature field around energy piles is less studied under typical intermittent operation and unbalanced thermal load conditions. There are few reports on the thermal behavior of the energy pile group, especially the influence of different arrangement modes on the thermal performance of the energy pile group. This is because each pile in the group is connected together through the platform, the load acting upon the pile is borne by both the pile and the foundation, and the load transfer mechanism of the pile group is different from that of a single pile and is more complex [26].

In this paper, a finite element numerical model was proposed and verified in order to identify the temperature field distribution around the energy pile and the thermomechanical behavior of the group. First, the heat conservation control differential and thermomechanical equation of energy pile were proposed, and a finite element numerical model, which enables the evaluation of the temperature field of energy pile and the thermomechanical behavior of the pile group, was developed. A full-scale in situ test of a large diameter over-length energy pile (1000 mm-diameter, 44 m-long) was carried out to verify the model. Second, the distribution of temperature field in pile body and surrounding underground under typical intermittent working conditions (heating load and cooling load unbalance conditions) were determined. The conceptual model of thermal load under typical intermittent operation conditions was proposed. Thermomechanical behavior characteristics of the energy pile group were determined, and the influences of arrangement on it were analyzed. This study provides a reliable and convenient method for the analysis of the temperature field of a single energy pile and the thermomechanical behavior of the pile group.

2. Materials and Methods

As a heat transfer component, the heat transfer process of the energy pile can be divided into (Figure 1): heat convection between circulating liquid and the inner wall of the pipe; heat conduction inside the pipe; heat conduction between the outer wall of the pipe and the pile concrete; heat conduction inside the pile concrete; heat conduction between the pile concrete and surrounding ground; heat convection/conduction in the surrounding ground. Based on the above heat transfer process and the geometric characteristics of the large aspect ratio of energy pile, the following assumptions were made in our finite element modeling: (1) the heat exchange pipe, concrete and surrounding underground are treated as homogeneous isotropic materials; (2) the influence of surface temperature change is not considered; (3) the underground around the energy pile is completely saturated without considering the seepage of groundwater; (4) the thermal strain of the pile in the thermal cycle is thermoplastic.

2.1. Temperature Field Governing Equation

The thermal behaviour of an energy pile is characterized as a transient process. The differential equation that governs heat transfer in the saturated ground is:

$$\rho c_{\mathrm{p}}\frac{\partial T}{\partial t}-{\rho }_w{c}_{w}k\nabla \Psi \nabla T=\nabla •\left(\lambda \nabla T\right)$$

(1)

Low permeable and hydraulic gradients are more favorable to ensuring the seasonal energy balance of the energy pile [3]. Therefore, groundwater-induced thermal convection is neglected, and the heat transfer can be regarded as a pure heat conduction process. Equation (1) can be simplified as follows:

$$\rho c_{\mathrm{p}}\frac{\partial T}{\partial t}=\nabla •\left(\lambda \nabla T\right)$$

(2)

In the three-dimensional Cartesian coordinate system and at a constant thermal conductivity, Equation (2) can be transformed into the following Equation (3):

$$\frac{\partial T}{\partial t}=\frac{\lambda }{\rho c_{\mathrm{p}}}\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)$$

(3)

where T is the transient temperature, °C; t is time, s; ρ is the concrete density, kg/m³; c_p is the mass-specific heat capacity at constant pressure, J/(kg·°C). The concrete thermal conductivity k, W/(m·K), can be assumed constant because it is not substantially affected by the temperature variation within the studied range.

In order to solve the differential equation above, boundary conditions need to be determined. There are three types of boundary conditions: (1) Dirichlet boundary condition, i.e., the temperature on the boundary is known or is a function of time. It is suitable for the situation if the specific temperature on the periphery of the surrounding ground can be determined after running for a period of time. (2) Neumann boundary condition, i.e., the heat flux on the boundary is known or is a function of time. For the heat convection between the heat exchange liquid and the inner wall of the pipe, it can be simplified (as shown in Equation (4)). (3) Cauchy boundary conditions, i.e., the temperature and heat transfer coefficient of the fluid medium in contact with the boundary, are known, as shown in Equation (5).

$$-k{\frac{\partial T}{\partial n}|}_{\Gamma }=q_2$$

(4)

where q₂ is the heat flux, W/m²; the outflow boundary q₂ is positive, which means the direction is from hot to cold; and the outflow boundary q₂ is negative.

$$-k{\frac{\partial T}{\partial n}|}_{\Gamma }=\alpha (T-T_{\mathrm{f}})$$

(5)

where α is the heat transfer coefficient, W/(m²°C); T_f is the fluid medium temperature, °C.

The transient temperature field under such complex geometric boundary conditions is unlikely to be solved using the analytical method; the finite element method is an alternative choice [20,22,24,25], which adopts the weighted margin to solve the mathematical problem:

$$\frac{\partial J^D}{\partial T_l}={\displaystyle \underset{D}{\iint }\left[k(\frac{\partial w_l}{\partial x}\frac{\partial T}{\partial x}+\frac{\partial w_l}{\partial y}\frac{\partial T}{\partial y})+\rho c_{\mathrm{p}}{w}_l\frac{\partial T}{\partial t}\right]}dxdy-{\displaystyle {\oint }_{\Gamma }k{w}_l}\frac{\partial T}{\partial y}{d}_{\mathrm{s}}=0$$

(6)

where w_l is the weighted function, n is the direction vector of the normal outer line at any boundary surface.

The boundary size can be estimated according to the far boundary theory:

$$R_{\infty }\ge 4\sqrt{{\alpha }_{\mathrm{s}}t}$$

(7)

where R_ꝏ is the influence radius, m; α_s is the average heat transfer coefficient, W/(m²°C).

The average heat flux of the inner wall of the heat exchange pipe can be converted according to the measured amount of heat transfer:

$$f=\frac{1000Q}{H\cdot n\cdot \pi \cdot d_{\mathrm{in}}}$$

(8)

where f is the average heat flux of the inner wall of the heat exchange pipe, W/m²; H is the vertical depth of the heat exchange pipe, m; n is the number of pipes in the borehole; d_in is the inner diameter of the pipe, mm; Q is the total heat transfer of a single energy pile, W, which can be calculated by Equation (9):

$$Q=\frac{{\rho }_{\mathrm{w}}V{c}_{\mathrm{w}}\Delta T}{3.6}$$

(9)

where ρ_w is the density of water, kg/m³; V is the total water flow at the ground source side, m³/h; c_w is the specific heat capacity of water at constant pressure, J/(kg°C); ΔT is the difference of water temperature inlet and outlet of the pipe, °C.

In order to consider the different heat transfer between the inlet and outlet pipes, the heat transfer coefficient of the inlet side (η) is introduced, which represents the ratio of the total heat transfer amount of the inlet pipe to the total heat transfer. The heat flux of the inlet and outlet pipes is obtained as follows:

$$f_{\mathrm{in}}=2f\eta$$

(10)

$$f_{\mathrm{out}}=2f\left(1-\eta \right)$$

(11)

where f_in and f_out are the heat flux of the inlet and outlet pipes, kW/m².

2.2. Thermomechanical Governing Equation

The thermoplastic constitutive equation can be expressed as:

$$\left\{\epsilon \right\}={\left[D\right]}^{-1}\left\{\sigma \right\}\Delta T$$

(12)

$$S={\left\{\alpha \right\}}^T\left\{\sigma \right\}+\frac{\rho c_{\mathrm{p}}}{T_0}\Delta T$$

(13)

where {ε} is the total strain tensor; [D] is the elastic stiffness matrix of the material; {σ} is the stress tensor; S is entropy density; {α} is the temperature expansion coefficient vector; T₀ is the initial temperature, °C.

According to the second law of thermomechanical:

$$Q=T_{0}S$$

(14)

Equation (12) can be rewritten as follows:

$$\left\{\sigma \right\}=\left[D\right]\left\{\epsilon \right\}-\left\{\beta \right\}\Delta T$$

(15)

where {β} is the thermoplastic coefficient tensor.

Substituting Equation (13) into Equation (14):

$$Q=T_0{\left\{\beta \right\}}^T\left\{\epsilon \right\}+\rho C_{\mathrm{v}}\Delta T$$

(16)

Substituting Equation (16) into the first law of thermomechanical can be obtained:

$$\frac{\partial Q}{\partial t}=T_0{\left\{\beta \right\}}^T\frac{\partial \left\{\epsilon \right\}}{\partial t}+\rho C_{\mathrm{v}}\frac{\partial \left(\Delta T\right)}{\partial t}-\left[K\right]{\nabla }^{2}T$$

(17)

where [K] is the thermal conductivity matrix of the material, expressed as:

$$[K]=\left[\begin{array}{ccc}{K}_{xx}& 0& 0\\ 0& K_{yy}& 0\\ 0& 0& K_{zz}\end{array}\right]$$

(18)

where K_xx, K_yy and K_zz are the thermal conductivities of ground and concrete in 3 different directions. In fact, K_xx=K_yy=K_zz (isotropic case) because of the little practical differences be-tween the directional conductivities.

The initial temperature is 20 °C, which applies to the lower boundary of the model, while adiabatic boundary conditions are applied to the far boundary of the upper and horizontal of the model. The initial stress field is gravitational. The pile top load is the external force load that corresponds to the safety factor K = 2. The pile is subject to a tem-perature of 40 °C under summer conditions and 0 °C during winter..

2.3. Modeling Process

The thermomechanical behavior of the energy pile is a one-way mechanism. The temperature field affects the stress field, whereas the change in the stress field does not affect the temperature field. In this work, the ANSYS APDL parameter design language was used for modeling (a scripting language); it can build models using parameterized variables that allow for complex data inputs, such as the size of the model, the performance of the material, the load, the position of the boundary conditions, and the density of the mesh. In the simulation of the thermomechanical coupling characteristics under vertical loads, the superstructure constraint is simplified as an equivalent elastic body with a diameter equal to the diameter of the pile, a unit length, a small mass and a certain stiffness. The specific parameters are as follows: thermal conductivity is 1.80 W/(m∙K), density is 1.0 kg/m³, the elastic modulus is 0–140,500 MPa, Poisson’s ratio is 0.1 and volume expansion coefficient is 1 × 10⁻¹⁰.

Under the action of temperature load, the pile top rises or settles so that the equivalent elastic body is stretched or compressed, and the load on the pile top changes. Therefore, it is necessary to ensure that the temperature field of the equivalent elastomer does not change during the simulation; otherwise, the function of the equivalent elastomer will be affected. In this study, the equivalent elastomer was modeled separately from the pile. When the temperature load was analyzed, the pile-surrounding rock and soil medium system were analyzed separately. When the thermodynamic coupling load was analyzed, the joint at the bottom of the equivalent elastic body and the joint at the top of the pile were coupled together so that they have the same degree of freedom so as to ensure that the temperature field of the equivalent elastic body does not change. Under vertical load, the geometric shape, temperature and external load of the energy pile are axisymmetric, and the model can be simplified as a two-dimensional axisymmetric model. For temperature field and structure analysis, a 2D four-node solid temperature unit (PLANE55) and a 2D four-node solid structural unit (PLANE182) were used, respectively, which were divided into 59,136 units and 59,925 nodes. The model process is shown in Figure 2.

2.4. Model Verification

In order to verify the reliability of the model, a full-scale in situ test of a large diameter over length energy pile (1000 mm-diameter, 44 m-long) was conducted in Nanjing, China. The heat exchange pipe is high-density polyethylene (HDPE) pipe connected in with 3U formalization, and the dimensions are 25 mm × 2.3 mm. Multiple sensors were fixed to different parts of the test pile to capture relevant parameters (Figure 3).

The finite element model of the test pile was established, and Figure 4 compares the simulated and measured results of four different feature locations of the reinforcement cage, 1.2 m, 2.1 m and 3.0 m from the center of the pile. Although the temperature values of different positions are slightly different, the simulated and measured results show very similar variation trends, which reflects the rationality of the model. The reasons for the deviation are: (1) the heat exchange pipe in the energy pile is loose and not completely fixed on the steel cage, which leads to the deviation of the position of the heat exchange pipe; (2) the soil layer is saturated, and the influence of heat convection in free water in soil pores is not considered. Therefore, the finite element model can reflect the spatial and temporal evolution characteristics of energy pile operation without considering the factors of groundwater heat convection.

Figure 5 reports the relationship between the pile axial force calculated according to the model and the vertical load on the pile top. It can be found that the vertical load exerted by the pile top continuously spreads to each layer through the pile body, and the axial force of the pile body decreases along the depth. It is noted that the pile bottom load is almost zero; even if the load reaches 1.5 times the ultimate capacity (8000 kN), the pile tip still does not increase significantly. It coincides with the law of pile load transfer [27], which indicates that the numerical model is reliable.

3. Results

3.1. The Evolution of Energy Pile Temperature Field under Typical Working Conditions

3.1.1. The Temperature Field Evolution under Typical Intermittent Condition

In a complete operation cycle, the intermittent operation is used to maintain the heat balance between the heat taken from underground and the heat dissipating, thus benefiting the long-term, efficient and stable operation of the system [28,29]. Furthermore, from the perspective of environmental friendliness and efficient utilization, the daily intermittent operation is a better choice [30,31]. Therefore, the daily thermal load is simplified, as shown in Figure 6.

Based on the actual data from the in situ tests, the simulation condition is set as follows: the initial ground temperature field is 20 °C; the simplified model of thermal load is adopted; the average heat transfer power in summer and winter conditions are 171 W/m² and −171 W/m², respectively; the heat transfer coefficient of the inlet side (η) is 0.6; and the heat flux of the inlet and the outlet pipe are 205.2 W/m² and 136.8 W/m², respectively. Figure 7 shows temperature field evolution over time under summer conditions. It can be seen that heat gradually diffused outward through heat exchange liquid and heat exchange pipe, energy pile and surrounding ground. The temperature of the pile and surrounding ground gradually increased. Additionally, the thermal influence area gradually expands with time, but the thermal influence area does not expand infinitely, and the final thermal influence radius is about 4 m.

Figure 8a exhibits the temperature changes over time at different locations under intermittent working conditions in summer. The temperature at different locations changes periodically and increases. The inner wall temperature is greater than the external wall, and a lower temperature can be observed if the test point is far from the center of the pile. Moreover, the temperature in the ground more than 1 m away from the center of the pile can no longer reflect the characteristics of intermittent thermal load, and its temperature slowly rises. Figure 8b reports the average temperature of the pile under intermittent working conditions in summer. It can be seen from the results that the average temperature of the pile gradually increases with time, while the difference between the highest and lowest temperature of the pile is relatively stable, varying between 0 and 4.2 °C, and the maximum temperature difference of the pile occurs when the first cycle operation stops.

Figure 9 reveals the temperature distribution of the pile body after intermittent operation for 90 days (outage ratio = 0.5) and continuous operation for 45 days under the condition of the same total heat transfer. Intermittent operation-induced temperature change is less than that of continuous operation; specifically, the temperature is 5.6 °C lower after 90 days of intermittent operation than after 45 days of continuous operation. Furthermore, intermittent operation can reduce the maximum temperature and reduce the temperature inhomogeneity of the pile.

3.1.2. The Temperature Field Evolution under Typical Heating Load and Cooling Load Imbalance Condition

The increase (or decrease) in the ground temperature field affects the heat transfer performance of the system and eventually inhibits energy saving or even causes the failure of the heat transfer system [32]. In order to investigate the influence of different geographical regions and the imbalance rate of cooling and heating load on the ground temperature field, three kinds of scenarios were designed. Scenario 1 operated for 120 and 60 days in summer and winter, respectively; Scenario 2 operated for 90 and 90 days in summer and winter, respectively; and Scenario 3 operated for 60 and 120 days in summer and winter, respectively.

Figure 10a shows the temperature of the pile–soil contact surface under different scenarios. It can be found that the proportion of cooling and heating load in an annual cycle significantly affects the change in the temperature field. Specifically, under the condition of the cooling and heating load balance (scenario 2), the periodic change in ground temperature within a certain range of initial ground temperature does not affect the long-term operation of the heat transfer system. Under the condition of heating load dominance (scenario 1), the ground temperature cannot be restored to the original ground temperature within an annual cycle, leading to a continuous rise in temperature, and at this time, the temperature of the inlet and outlet of the heat exchange pipe is increased. Under the condition of cooling load dominance (scenario 3), the opposite is true. Figure 10b reports the ground temperature recovery at different years under the three scenarios. It can be found that the ground temperature can be completely recovered under the condition of balance, but it cannot be completely recovered under the condition of heating load (or cooling load) dominance, which leads to cold and hot accumulation after the long-term operation (10y). In practical engineering, the imbalance of heating and cooling load is widespread. If not intervened, it would not only reduce the heat transfer efficiency of the system but also seriously affect sustainable utilization [33]. Generally, when the heating load is dominant, auxiliary measures should be adopted, such as adding cooling towers; when the cooling load dominates, joint heating measures such as solar energy can be taken.

3.2. Thermomechanical Behavior of Energy Pile Group under Vertical Load

In this work, nine piles of low platform pile groups were taken as an example; the specific dimensions and finite element mesh are shown in Figure 11. The pile and the ground near the pile are areas where high strain may occur, so a finer mesh is used. The center of the pile group platform was subjected to a rectangular distribution load of 1.6 m× 1.6 m. During the simulation of pile thermomechanical behavior under vertical load, the constraint of the superstructure is simplified to an equivalent elastic body with a certain stiffness whose diameter is equal to that of the pile. Under the action of temperature load, the pile top was uplifted or settled, and the equivalent elastic body was thus stretched or compressed, which changes the load on the pile top. In ANSYS simulation, there was still a problem that needed to be solved: we needed to ensure that the temperature of the equivalent elastomer did not change in the temperature field simulation; otherwise, the function of the equivalent elastomer would be affected. In order to solve this issue, the equivalent elastomer and pile were simulated separately. When the temperature field was analyzed, the pile-ground system was analyzed separately. The node at the lower end of the equivalent elastic body was coupled with the node at the top of the pile to investigate thermomechanical behavior. Moreover, under vertical load, the geometric shape, temperature and external force load of the energy pile were all axisymmetric, so the model could be simplified to a two-dimensional axisymmetric model. Figure 11 is a schematic diagram of the simulated energy pile group. According to the position of the pile under the platform, the energy pile can be divided into three types: center pile (pile 5), corner pile (pile 1/3/7/9), and side pile (pile 2/4/6/8).

Figure 12a,b shows the temperature axial force of the pile body when all the nine piles are used as energy piles. In the summer condition, the temperature of the pile increases, and the axial compressive force generates; in the winter condition, the temperature of the pile decreases, and the axial tensile force occurs. The magnitude and distribution of the temperature axial force are related to the position of the energy pile. The center pile, side pile and corner pile show different characteristics. The temperature axial force generated by the center pile (pile 5) is evenly distributed along the pile. The temperature axial force of the corner pile (1/3/7/9) and side pile (2/4/6/8) is small at both ends and large in the middle, and the maximum value is 35 m below the pile top. Under the same temperature load condition, a smaller degree of freedom (DOF) can generate a greater temperature force, and the DOF is related to the constraints of the pile tip and surrounding ground. The DOF of the pile group under winter and summer conditions are shown in Figure 12c,d. It can be found that the variation in pile temperature does not significantly change the size and distribution of DOF. Specifically, the upper part of the corner pile (pile 1/3/7/9) has the largest DOF, followed by the side pile, and the center pile has the smallest DOF, indicating that different positions of piles under the platform have different DOF. The DOF of each pile varies in the range of 0.985~1, suggesting that when the pile group consists of an energy pile, the constraint between piles is very small, thus inducing smaller additional stress. Moreover, the pile group is heated or cooled at the same time, and the geological conditions around the pile are relatively consistent, resulting in relatively consistent displacement, as shown in Figure 12e,f.

3.3. Influence of Layout on Thermomechanical Behavior of Energy Pile Group

If only one pile in the group is heated, the DOF ranges from 0.4 to 0.5 [19,34]. This indicates that the number and distribution of energy piles can affect the DOF of energy piles and non-energy piles, making the constraint degree of each pile complicated. Additionally, there is a certain damage rate during the construction of energy piles, which may change the number and distribution of energy piles. In order to investigate the influence of layout mode or damage rate on the internal forces of energy pile and non-energy pile, six different scenarios are designed based on the pile group foundation of nine piles, as shown in Figure 13.

Figure 14a shows the axial temperature stress distribution of the center pile under different layout scenarios. When the center pile is used as an energy pile, the compressive stress generates, and the less used energy piles result in higher compressive stress. Especially in scenario 6 (only the center pile is the energy pile), while the temperature of the other eight piles increased, the compressive stress is smaller than that of the energy pile. The other eight piles became “anchor piles”, with the maximum temperature compressive stress at the top of the center pile reaching 1 MPa, accounting for about 3.33% of the compressive strength of concrete (C30). Additionally, when the center pile is used as a non-energy pile (scenario 3/4/5), tensile stress is generated in scenarios 3 and 5, which is distributed evenly along the pile body. If more energy piles are applied, greater tensile stress can be expected. The maximum tensile stress generated in scenario 3 is about 20 kPa, accounting for about 0.67% of the tensile strength.

Figure 14b reports the axial temperature stress distribution of the corner pile under different layout scenarios. When the corner pile is an energy pile, the compressive stress occurs; when it is a non-energy pile, the tensile stress can be expected. In scenario 1, the temperature stress is the smallest, and the stress at both ends is almost zero. The maximum temperature stress (97.8 kPa) occurs at about 36 m below the pile. The maximum compressive stress is generated in scenario 5, which reaches 822 kPa (4.1% of the compressive strength of C30 concrete) at about 2 m from the pile top. The corner pile (1/3/7/9) in scenario 4 and the non-energy corner pile (3/7) in scenario 5 induced the maximum temperature tensile stress, which reached 448 kPa and 443 kPa at the top of the pile 1 m, respectively. Moreover, the corresponding total stress of the pile was 7.7 MPa and 7.7 MPa, respectively. Therefore, the stress induced by the corner pile does not destroy the pile or reduce its durability of the pile under various scenarios.

Figure 14c indicates the axial temperature stress distribution of the side pile under different layout scenarios. The compressive stress is generated in scenarios 1 and 4, especially in scenario 4, where the compressive stress at the top of the pile reached 454 kPa, accounting for 6.8% of the total stress. In the scenario of a non-energy pile, the temperature tensile stress is generated. In scenarios 2 and 3, the tensile stress on the pile top reaches 235 kPa and 238 kPa, respectively, which reduces the compressive stress on the pile top by 3.7%, and the pile body has no tensile stress.

From the above scenario analysis, we understand that the layout and number of energy piles significantly affect the forces distribution of each pile. In the summer condition, the energy pile produces the compressive stress, but the non-energy pile acts as the anchor pile, and the pile body produces the tensile stress. When some piles are not used as energy piles, the temperature compressive stress of energy piles increases significantly. An increased number of non-energy piles can result in temperature compressive stress. The ratio of the induced temperature compressive stress and tensile stress to the strength of concrete ranges from 0.33% to 4.53% and 6.23% to 14.93%, respectively. In the winter conditions, the energy produced with tensile stress, rather than the energy produced with compressive stress, this part of the compressive stress part from pile body temperature to reduce the induced tensile stress; the other part comes from the pile top load generated by the platform, the temperature of the induced tensile stress and compression stress of the concrete strength ranges from 3.27% to 45.33% and 0.78% to 1.49%, respectively.

4. Conclusions

Based on the heat transfer mechanism of energy piles, the differential equation governing the heat transfer of energy piles was proposed. A finite element numerical model was developed and verified, and it can be applied to analyze the unsteady-state temperature distribution and the thermomechanical behavior of the energy pile group. The temperature field evolution of the energy pile under typical intermittent working conditions and heating/cooling load imbalance conditions were analyzed using this model. Moreover, the thermal behavior of the energy pile group was investigated, and the influence of the layout mode on the thermomechanical behavior of the energy piles group was analyzed in six different scenarios. Specific conclusions are as follows:

(1)

Under typical intermittent operation conditions, the temperature of the energy pile and surrounding ground near the heat exchange pipe varies periodically. The closer the distance to the heat exchange pipe, the greater the temperature gradient and the more sensitive to the change in heating load. When the distance from the center of the pile exceeds 1 m, the characteristics of intermittent operation can no longer be reflected;

(2)

For areas with unbalanced cold and heating load, long-term operation of energy piles leads to thermal accumulation, and the maximum temperature of energy piles occurs in the first daily cycle. The cooling and heating load should be leveled by an auxiliary cold source (cooling tower) or auxiliary heat source (solar collector), thus controlling the fluctuation of ground temperature within a reasonable temperature range;

(3)

When the pile group is all energy piles, the DOF at the pile top is the highest, and the pile top load at different positions decreases in the order of corner pile > side pile > center pile. In summer/winter working conditions, the increase/decrease in pile temperature induces axial compression/tensile stress. The axial temperature stress induced at the top of the center pile and the bottom of the corner pile is the largest. In addition, the DOF of each pile varies in the range of 0.985~1, and the constraints between piles are very small;

(4)

When the pile group is partially used as the energy pile, the non-energy pile acts as the “anchor pile”, and it generates the added tensile stress. Decreased number of energy piles generates lower DOF of center piles and induces higher temperature stress. The stress induced by the corner pile does not destroy the pile or reduce its durability. Moreover, due to the constraint of the platform, the bending moment of piles at different positions is different, and the temperature stress increases with the increase in pile spacing.