Analysis of Bearing Characteristics of Energy Pile Group Based on Exponential Model

: Energy piles are an emerging energy technology for both structural and thermal purposes. To support structure load, piles are always used in groups with raft; however, the cost and complexity of ﬁeld tests and numerical modelling limits the research on the bearing characteristics of energy pile group. In this paper, exponential model was applied to simulate the thermo-mechanical soil-pile interaction of energy pile group. Axial load transfer ( τ - z ) analysis was ﬁrst performed to calculate the shear stress distribution in the soil, then matrix displacement method was introduced to determine the thermo-mechanical response of energy pile group. The validity of the analytical model was tested against ﬁeld tests and numerical results. A case study was further performed to analyze the inﬂuence of thermal cycles and arrangement of thermally active piles on the bearing response of the whole pile group. Test results show that with the thermally activated pile in pile group, (1) differential settlement increases with thermal cycle numbers; (2) the axial force of thermally active pile increases during heating process and decreases during cooling process, and this trend varies for the surrounding nonthermal piles; (3) induced load on thermal pile increases with thermal cycles, but decreases for nonthermal piles. The proposed analytical model is expected to serve as a simple and convenient alternative for the preliminary analysis on the bearing characteristics of energy group pile.


Introduction
Energy pile is a novel and environmentally friendly foundation type that achieves the dual role of supporting superstructure and transferring heat between structure and ground. Due to increasing global demand for clean and green energy, research on the thermomechanical behavior of energy piles is gaining increasing popularity [1,2]. Most of the work focuses on the thermal performance and optimal layout of energy pile, while research on its function as structure support has been less reported [3]. Actually, energy piles are under the coupling influence of mechanical and thermal loads. It has been widely reported that the continuous cooling and heating process would impose extra stresses and strains on the pile structure and surrounding soil, thus undermining the long-term safety of the structure [4]. A better understanding of the influence of thermal cycles on the mechanical behavior of energy pile will promote the development and optimization of design guidelines for this new foundation type. For most civil structures, the foundation type always involves groups of piles connected by a raft, rather than a single pile. The thermomechanical behavior of energy pile group is very complex due to the interaction effect between energy piles, nonthermal piles, pile raft and surrounding soil, which adds experimental cost and modelling difficulty to the study of the pile group system [5,6].
Experimental studies, especially field full-scale tests, are by far the most reliable and direct way to observe and analyze the thermomechanical behavior of energy pile groups. Based on field experiments, Mimouni and Laloui [7] revealed that heating the entire pile group could reduce the differential displacements compared to the energy pile group with non-active piles. Peng et al. [8] conducted thermal cycle tests on the floating model energy pile group and single energy pile. The thermomechanical response, including axial pile stress and pile displacement, was quite different between the single energy pile and the pile group. This was explained by the "group effect", which weakened the lateral confining pressure in the pile group. Fang et al. [9] performed field tests on a 2 × 2 pile-raft foundation. Three thermal loading schemes were performed to consider the influence of location arrangement and number of energy piles on the thermomechanical behavior of the pile-raft system. Test results indicated that the heating process would induce thermal expansion for both the energy pile and its surrounding non-active pile, and this influence differed for adjacent piles and diagonal pile. Also, the nonuniform thermal loading scheme would induce differential settlement of the pile raft foundation and the differential settlement was less obvious with increasing active piles in the pile group.
Due to the expensive and time-consuming experimental setup, valid test data are relatively scarce, and numerical modelling provides a convenient and promising alternative to investigate the soil-pile interaction under thermal and mechanical load [10]. Jeong et al. [11] conducted numerical parametrical analysis on the design parameters of energy pile group, including pile spacing, pile location, soil type, and bearing condition. The pile cap seemed to limit the vertical pile displacement during the cooling process. Based on three-dimensional thermomechanical numerical model, Salciarini et al. [12] revealed that the maximum axial load difference occurred at around one month after thermal cycles. Di Donna and Laloui [13] adopted a fully coupled thermo-hydro-mechanical Lagrangianupdated formulation to simulate the thermomechanical behavior of single energy pile and pile group. Though the thermal-induced stress and strain were in the acceptable range, it was recommended to consider the influence of thermal cycles in the design of energy piles. Saggu and Chakraborty [14] analyzed the influence of pile spacing and arrangement of active piles on the thermomechanical response of energy pile group. The cyclic thermal load would induce higher axial displacement and stress compared to the condition with only mechanical load. For numerical case with pile spacing of 2.5 m, a 47% increase in displacement was observed for thermal piles compared to piles under mechanical load only and a recommendation of at least 4 m spacing was finally provided.
The above research provides valuable information about the thermal influence on the mechanical behavior of energy pile group. However, numerical techniques still have shortcomings in simulating complex and intricate models. As a result, a more convenient and simplified approach is necessary to facilitate the practical application. Knellwolf et al. [15] proposed a calculation method of the temperature effect on the energy pile behavior. Onedimensional finite difference model was built to consider the variation of soil properties at different layers. The load-transfer curve proposed by Frank and Zhao [16] was used to simulate the pile-soil interaction, and the thermal loading behavior was determined by the selection of null point. Validity of the proposed method was tested against field test data. Further analysis indicated that the heating process would produce extra compression of the pile, which in turn mobilized more shear stress, while the cooling process led to decreased shear stress and the development of tensile stress in the pile. Chen and McCartney [17] extended the traditional load transfer (τ-z) method by incorporating the thermo-elastic deformation of the foundation the null point criterion. The thermomechanical τ-z method was used to evaluate the influence of different representative boundary conditions and achieved good results. To avoid the selection of null point under the influence of mechanical load and temperature, Fei et al. [18] applied a modified hyperbolic function to simulate the pile-soil interaction and the thermal effect was characterized by thermal load. Though satisfactory results were achieved, the hyperbolic model would produce hysteresis loop during the loading-unloading process; as a result, a reduction factor was needed to simulate the progressive softening behavior of soil shear strain. The above models lay a solid foundation for the analytical analysis of energy group pile. However, further modifications are still necessary to construct models that simplify the selection of input parameter while reflecting the progressive evolvement of soil property.
The objective of this paper is to provide a simple and convenient analytical approach to evaluate the thermal effect on the bearing characteristics of energy group pile. Exponential model is introduced to simulate the nonlinear behavior at the pile-soil interface, and the thermomechanical behavior is simplified as mechanical and thermal load. The validity of the analytical model is checked against reported pile group data. On this basis, the influence of thermal cycle numbers and the arrangement of thermally active piles on the bearing performance, in terms of axial force and settlement, of the energy pile group is evaluated.

Pile Shaft Response
Settlement calculation for pile group is more complex compared to that for single pile due to the interaction of pile-soil-pile-raft. To simplify the calculation, settlement of pile group is always divided into two parts, settlement of single pile and induced settlement caused by group effect.
Under vertical load, a certain range of soil around the pile would enter the plastic stage, while soil outside this range is still in an elastic state. It is acknowledged that the nonlinear load-settlement response of single pile depends on the relative displacement at the pile-soil interface. As a result, the pile displacement u at certain depth is decomposed into the relative displacement s at the pile-soil interface and the displacement of the surrounding soil w, as The shear stress and relative displacement at pile shaft are calculated based on the exponential load transfer model [19,20], as where τ is the shear stress at the pile shaft, γ is the relative displacement at the pile-soil interface, a and b are model parameters, which are defined as where τ ult is the ultimate shear stress at the pile-soil interface, G max is the initial shear stiffness at the pile-soil interface. The elastic displacement of the surrounding soil w is calculated according to Randolph and Wroth [21] as where r is the pile radius, r m is the influence radius around the pile, G is the shear modulus of the soil. Combining Equations (1), (2) and (5) produces the relationship between incremental displacement and shear stress as Equation (6) is normally used for monotonic loading. As for energy piles under thermal loading, heating and cooling cycles would lead to the continuous loading and reloading behavior of the pile-soil system. To simulate the cyclic behavior, a load transfer model considering cyclic load effect is constructed according to Masing criterion [18].
At unloading stage, At reloading stage, where τ 0 and γ 0 are the shear stress and relative displacement at the transition point between loading and unloading stage, a 1 , b 1 and a 1 , b 1 are model parameters, which are defined as By differentiating Equations (7) and (8), the relationship between incremental displacement and shear stress under cyclic loading, as shown in Figure 1, is described mathematically as: At loading stage, At reloading stage, Equation (6) is normally used for monotonic loading. As for energy piles under thermal loading, heating and cooling cycles would lead to the continuous loading and reloading behavior of the pile-soil system. To simulate the cyclic behavior, a load transfer model considering cyclic load effect is constructed according to Masing criterion [18].
At unloading stage, By differentiating Equations (7) and (8), the relationship between incremental displacement and shear stress under cyclic loading, as shown in Figure 1, is described mathematically as: At loading stage, At reloading stage,

Pile Base Response
The bearing stress and relative displacement at pile base are calculated based on the exponential load transfer model as where p is the base resistance and γ b is the relative displacement at the pile base, a 2 and b 2 are model parameters, which are defined as where τ ult is the ultimate shear stress at the pile base; G max is the initial shear stiffness at the pile base. And the relationship between incremental displacement and shear stress under cyclic loading is At unloading stage, At reloading stage,

Pile-Pile Interaction
For pile i in a pile group, the relative displacement at the pile shaft could be decomposed into the nonlinear displacement induced by the shear stress at the pile shaft and the induced elastic displacement caused by the m piles in the influence range around pile i, as where δu iiz is the displacement caused by the shear stress at pile i, which could be calculated by Equation (19). According to Lee and Xiao [22], δu iiz can be determined as where S ij is the distance between pile i and pile j; when S ij > r m , the influence of surrounding pile could be ignored. The coefficient r m is influenced by the pile length and soil layer, and its value tends to increase with depth, which could be calculated as where δ is the soil shear modulus, µ s is poisson ratio, L is pile length. In pile group, assuming the total pile number is n, then the displacement of single pile i is According to Randolph and Wroth [21], the displacement increment at the pile base caused by pile interaction could be seemed as elastic. On this basis, the base displacement increment of pile i in a pile group could be calculated as

Solution of the Analytical Model
In general, the energy pile is first subjected to the load of superstructure. When the structure needs to be cooled or heated, the energy pile is then put into operation and begins to suffer the cyclic heating and cooling process [23]. On this basis, the thermomechanical process is simplified as mechanical loading and thermal loading.

Governing Equation for Single Pile
According to the equilibrium state of pile element, the following equation is obtained where N is the pile axial force, z is depth.
Considering the thermal effect, the axial strain of the pile can be calculated as where ε is axial strain, E p is Young's modulus of the pile, A is the cross-section area of the pile, a T is the thermal expansion coefficient of the pile, and ∆T is the temperature variation. By differentiating Equation (30) and combining Equation (29), one obtains In Equation (31), the basic assumption is that the temperature increment keeps the same along the pile length. The pile response under mechanical loading is first solved, the solution then serves as the initial condition at thermal loading stage.

Solution of Governing Equation under Mechanical Load
For a single pile with length of L, it is discretized into N sections, with each section a length of ∆h and a total of M = N + 1 nodes. By combing Equations (25) and (31), the solution based on 2nd order difference formula is To meet the solution accuracy for the pile top and bottom nodes, two fictitious nodes are added at position 0 and M + 1, and the formula could be solved considering the boundary condition.
Due to the raft constraint at pile head, the displacement values at pile top could be assumed as equal as long as the pile raft is stiff enough, and the external load could be assumed as the sum of the load on each pile. According to this assumption, the following equations are obtained δu 1,1 = δu 2,1 = · · · = δu n,1 (34) Combining the above equations with Equation (30), Equation (35) could be present in the form of Equation (36) based on the 1st order difference formula, that is Similarly, the iterative relationship of node M at pile base is derived as Combining Equations (33), (34), (36) and (37) in a matrix form, the governing equation for pile group is obtained as The stiffness matrix [K i,j ] is a M + 2 order square matrix, and the non-zero elements are Elements at the first and M + 2 rows are fictitious nodes, which can be determined as

Solution of Governing Equation under Thermal Load
For thermal loading simulation, the solution under mechanical load is used as the initial condition, and the governing equation keeps the same form. For energy pile group with thermally-active piles and nonactive piles, piles at different locations would experience different load condition, which would lead to different pile response. According to the research of Fei et al. [18], the bottom of the pile raft is assumed as plane. The center point of pile 1 is set as the origin coordinate, the direction of x axis is along the connecting line of pile 1 and pile 2, the direction of y axis is along the connecting line of pile 1 and pile 3, as shown in Figure 2. with thermally-active piles and nonactive piles, piles at different locations would experience different load condition, which would lead to different pile response. According to the research of Fei et al. [18], the bottom of the pile raft is assumed as plane. The center point of pile 1 is set as the origin coordinate, the direction of x axis is along the connecting line of pile 1 and pile 2, the direction of y axis is along the connecting line of pile 1 and pile 3, as shown in Figure 2. According to the assumption of stiff deformation, displacement at the top of pile i could be calculated based on the top displacement of pile 1, 2, 3 as where xi and yi are the center coordinates at pile top. As the boundary condition at pile top changes during thermal loading, the first row of the stiffness matrix [Ki,j] is cleared and the value is reassigned according to Equation (42). According to the assumption of stiff deformation, displacement at the top of pile i could be calculated based on the top displacement of pile 1, 2, 3 as where x i and y i are the center coordinates at pile top. As the boundary condition at pile top changes during thermal loading, the first row of the stiffness matrix [K i,j ] is cleared and the value is reassigned according to Equation (42).
In Equation (42), there are n − 3 equations. According to the assumption that the mechanical load keeps stable during thermal loading, Equation (41) is used to solve the governing equation of Equation (38).
Based on Equation (43), the left three rows of stiffness matrix [K i,j ] could be solved as Energies 2021, 14, 6881 9 of 16 As for the boundary condition at pile base, the effect of thermal loading is introduced into Equation (37) to produce The result of stiffness matrix is uninfluenced by the modification of boundary condition. However, with the introduction of thermal effect, there will be extra term in Equations (43) and (45), which may influence the load matrix in the governing equation set. To solve this problem, the load matrix at mechanical load stage is modified as To meet the calculation accuracy, the midpoint increment method [18] is used to solve Equation (38) to obtain the node displacement at different section of the pile. Then Equations (30) and (31) are used to calculate the shaft resistance and pile axial force.

Case with Mechanical Loading
The first case is from O'Neill et al. [24], where a pile group, comprising 9 closed-ended piles in 3 × 3 setup, was loaded to failure in stiff over-consolidated clay. The pile group configuration is presented in Figure 3. The pile had a diameter (d) of 274 mm, a wall thickness (t) of 9.3 mm, an embedment depth (H) of 13.1 m and a center-to-center spacing (s) of 822 mm. The side length (B) of the square raft was 2.7 m. Castelli and Maugeri [25] evaluated the field soil properties based on the back analysis of the data from O'Neill et al. [24]. According to Castelli and Maugeri [25], the shear modulus for soil (E s ) close to the ground surface was 47.9 MPa, which increased linearly with depth to 151 MPa at the pile base. The unit shaft resistance (f s ) increased linearly from 19 MPa to 93 MPa. The ultimate unit base resistance (q b ) was 2.15 MPa and Poisson's ratio (µ s ) was set as 0.5. The soil deformation modulus (E) was set as 1.95 × 10 5 kPa and the Young's modulus (G) for the pile was 2.1 × 10 8 kPa.
Based on the above parameters, the load response of the pile group is calculated based on the proposed method. As shown in Figure 4, the trend of the calculated Q-S curve is basically consistent with the test data. At initial loading stage, pile group settlement increases linearly with load. At a load level around 5600 kN, a sudden drop occurs. This may indicate that the pile reaches its ultimate state, pile displacement increases rapidly, which leads to the failure of the pile group. Based on the Q-S curve, the ultimate capacity of the pile group is determined as 5600 kN.
To further analyze the load response of pile group foundation, load variations for different piles against the total load is presented in Figure 5. The calculated load shared by different piles shows good agreement with the measured value. Also, the difference of the load value at different locations is minimal; this may be due to the constraint effect of the raft. However, a difference still exists. It seems that piles located at the diagonal corner bear the largest load and pile at the center supports the minimum weight. When the total load is 809 kN, the loads supported by the center pile are 10.4% (calculated) and 10.2% (measured) smaller than that by the diagonal pile. When the total load is 5324 kN, these values are 8.6% (calculated) and 6.2% (measured), respectively. Based on the above analysis, it is believed that the proposed method could be used to simulate the load response of the pile group.
group configuration is presented in Figure 3. The pile had a diameter (d) of 274 mm, a wall thickness (t) of 9.3 mm, an embedment depth (H) of 13.1 m and a center-to-center spacing (s) of 822 mm. The side length (B) of the square raft was 2.7 m. Castelli and Maugeri [25] evaluated the field soil properties based on the back analysis of the data from O'Neill et al. [24]. According to Castelli and Maugeri [25], the shear modulus for soil (Es) close to the ground surface was 47.9 MPa, which increased linearly with depth to 151 MPa at the pile base. The unit shaft resistance (fs) increased linearly from 19 MPa to 93 MPa. The ultimate unit base resistance (qb) was 2.15 MPa and Poisson's ratio (μs′) was set as 0.5. The soil deformation modulus (E) was set as 1.95 × 10 5 kPa and the Young's modulus (G) for the pile was 2.1 × 10 8 kPa. Based on the above parameters, the load response of the pile group is calculated based on the proposed method. As shown in Figure 4, the trend of the calculated Q-S curve is basically consistent with the test data. At initial loading stage, pile group settlement increases linearly with load. At a load level around 5600 kN, a sudden drop occurs. This may indicate that the pile reaches its ultimate state, pile displacement increases rapidly, which leads to the failure of the pile group. Based on the Q-S curve, the ultimate capacity of the pile group is determined as 5600 kN.  To further analyze the load response of pile group foundation, load variations for different piles against the total load is presented in Figure 5. The calculated load shared by different piles shows good agreement with the measured value. Also, the difference of the load value at different locations is minimal; this may be due to the constraint effect of the raft. However, a difference still exists. It seems that piles located at the diagonal corner bear the largest load and pile at the center supports the minimum weight. When the total load is 809 kN, the loads supported by the center pile are 10.4% (calculated) and 10.2% (measured) smaller than that by the diagonal pile. When the total load is 5324 kN, these values are 8.6% (calculated) and 6.2% (measured), respectively. Based on the above analysis, it is believed that the proposed method could be used to simulate the load response of the pile group.  To further analyze the load response of pile group foundation, load variations for different piles against the total load is presented in Figure 5. The calculated load shared by different piles shows good agreement with the measured value. Also, the difference of the load value at different locations is minimal; this may be due to the constraint effect of the raft. However, a difference still exists. It seems that piles located at the diagonal corner bear the largest load and pile at the center supports the minimum weight. When the total load is 809 kN, the loads supported by the center pile are 10.4% (calculated) and 10.2% (measured) smaller than that by the diagonal pile. When the total load is 5324 kN, these values are 8.6% (calculated) and 6.2% (measured), respectively. Based on the above analysis, it is believed that the proposed method could be used to simulate the load response of the pile group.

Case with Thermomechanical Loading
The second case is from Jeong et al. [11], where a series of 3D numerical models were constructed to analyze the influence of thermal loading, pile spacing, active pile location, and pile end restraint on the bearing performance of energy pile group. The selected case is based on a 3 × 3 pile group, with 3 thermally-active piles and 6 non-active piles. The location of the energy pile is shown in Figure 6. The pile was made of pre-tensioned spun high strength concrete, with a diameter (d) of 500 mm and pile length (L) of 20 m. The side length (B) of the square raft was 6.5 m, with thickness (t) of 1 m. Pile heads were rigidly connected to the raft. The Young's modulus (G) of the numerical soil was 50 kPa, Poisson's ratio (μs') was 0.3, friction angle (φ) was 35°. The pile shaft stress (fs) was 120 kPa and the ultimate base stress (qb) was 1150 kPa; the soil shear modulus (Es) was 13.6 MPa. The val-

Case with Thermomechanical Loading
The second case is from Jeong et al. [11], where a series of 3D numerical models were constructed to analyze the influence of thermal loading, pile spacing, active pile location, and pile end restraint on the bearing performance of energy pile group. The selected case is based on a 3 × 3 pile group, with 3 thermally-active piles and 6 non-active piles. The location of the energy pile is shown in Figure 6. The pile was made of pre-tensioned spun high strength concrete, with a diameter (d) of 500 mm and pile length (L) of 20 m. The side length (B) of the square raft was 6.5 m, with thickness (t) of 1 m. Pile heads were rigidly connected to the raft. The Young's modulus (G) of the numerical soil was 50 kPa, Poisson's ratio (µ s ) was 0.3, friction angle (ϕ) was 35 • . The pile shaft stress (f s ) was 120 kPa and the ultimate base stress (q b ) was 1150 kPa; the soil shear modulus (E s ) was 13.6 MPa. The validation scenario is that the temperature is cooled down by 9 • C and the external load is 9 MN. The focus is on the influence of pile spacing on the development of axial force and settlement of pile group.  The influence of pile spacing on the distribution of axial force of energy piles is illustrated in Figure 7. At pile spacing of 3d, the diagonal pile supports the largest load while the center pile bears the minimum load. And the axial force deceases with increasing distance to the pile top. When pile spacing increases to 5d, the axial force of the center pile shows a noticeable increase while the axial forces for piles at the side and diagonal corner are basically at the same level. This may be explained by the weaker "group effect" at larger pile spacing [26,27]. When pile spacing is small, the center pile would be dragged by the surrounding pile. This dragging effect would limit the relative displacement among piles, which further influences the development of shaft resistance of the center pile. As pile spacing increases, the bearing capacity of center pile increases and bearing performance of the pile group is improved. Another finding is that at pile spacing of 5d, the calculated axial force is obviously smaller that the numerical results, especially at the lower part the pile. It seems that the analytical model produces better prediction of axial force at small pile spacing. This may be the result of simplification of the raft-soil interaction in the proposed method. When pile spacing is small, the pile group works as a whole and the influence of pile raft is less obvious. When pile spacing increases, the influence of raft-soil interaction needs to be considered. On this basis, the proposed model is recommended for pile raft foundation with small pile spacing. The influence of thermal load and pile spacing on pile top settlement is presented in Figure 8. Since piles are assumed to be rigidly connected to the raft, pile top displacement The influence of pile spacing on the distribution of axial force of energy piles is illustrated in Figure 7. At pile spacing of 3d, the diagonal pile supports the largest load while the center pile bears the minimum load. And the axial force deceases with increasing distance to the pile top. When pile spacing increases to 5d, the axial force of the center pile shows a noticeable increase while the axial forces for piles at the side and diagonal corner are basically at the same level. This may be explained by the weaker "group effect" at larger pile spacing [26,27]. When pile spacing is small, the center pile would be dragged by the surrounding pile. This dragging effect would limit the relative displacement among piles, which further influences the development of shaft resistance of the center pile. As pile spacing increases, the bearing capacity of center pile increases and bearing performance of the pile group is improved. Another finding is that at pile spacing of 5d, the calculated axial force is obviously smaller that the numerical results, especially at the lower part the pile. It seems that the analytical model produces better prediction of axial force at small pile spacing. This may be the result of simplification of the raft-soil interaction in the proposed method. When pile spacing is small, the pile group works as a whole and the influence of pile raft is less obvious. When pile spacing increases, the influence of raft-soil interaction needs to be considered. On this basis, the proposed model is recommended for pile raft foundation with small pile spacing.
The influence of thermal load and pile spacing on pile top settlement is presented in Figure 8. Since piles are assumed to be rigidly connected to the raft, pile top displacement is the same under sole mechanical load. When piles are under thermal load, the cooling process would lead to contraction of the pile body, thus increasing pile deformation. It can be seen that the largest settlement occurs on the diagonal pile, this is in agreement with the numerical result. When pile spacing decreases to 3d, pile settlement increases. This is the result of group effect that limits the full development of soil resistance, especially for soil around the center pile. calculated axial force is obviously smaller that the numerical results, especially at the lower part the pile. It seems that the analytical model produces better prediction of axial force at small pile spacing. This may be the result of simplification of the raft-soil interaction in the proposed method. When pile spacing is small, the pile group works as a whole and the influence of pile raft is less obvious. When pile spacing increases, the influence of raft-soil interaction needs to be considered. On this basis, the proposed model is recommended for pile raft foundation with small pile spacing. The influence of thermal load and pile spacing on pile top settlement is presented in Figure 8. Since piles are assumed to be rigidly connected to the raft, pile top displacement is the same under sole mechanical load. When piles are under thermal load, the cooling process would lead to contraction of the pile body, thus increasing pile deformation. It can be seen that the largest settlement occurs on the diagonal pile, this is in agreement with the numerical result. When pile spacing decreases to 3d, pile settlement increases. This is the result of group effect that limits the full development of soil resistance, especially for soil around the center pile.

Simulation Background
In practice, there will be thermally active piles and non-active piles in an energy pile group. The influence of thermal cycles on the loading and deformation characteristics of active piles has been studied based on the above analysis. However, the influence of thermal loading on the bearing behavior of non-active piles and the whole pile group is not clear.
In this section, the proposed analytical approach is further applied to analyze the influence of thermal activity on the bearing behavior of energy pile group based on the field test parameters from O'Neill et al. [24]. The geometrical and mechanical parameters of the pile raft system and the soil foundation keep the same as in Section 4.1, while the pile sequence number is modified according to Figure 2. The main difference is that pile 1, which locates at the diagonal corner, serves as a thermally active pile. Pile 1 suffers 15 cycles of heating and cooling and the total load on the pile raft is 3000 N, which is 52.9% of the ultimate bearing capacity of the pile group. Detailed information of the loading scheme is provided in Figure 9. The focus of this simulation is to reveal the thermal influence on the development of axial force and settlement of non-active piles.

Simulation Background
In practice, there will be thermally active piles and non-active piles in an energy pile group. The influence of thermal cycles on the loading and deformation characteristics of active piles has been studied based on the above analysis. However, the influence of thermal loading on the bearing behavior of non-active piles and the whole pile group is not clear.
In this section, the proposed analytical approach is further applied to analyze the influence of thermal activity on the bearing behavior of energy pile group based on the field test parameters from O'Neill et al. [24]. The geometrical and mechanical parameters of the pile raft system and the soil foundation keep the same as in Section 4.1, while the pile sequence number is modified according to Figure 2. The main difference is that pile 1, which locates at the diagonal corner, serves as a thermally active pile. Pile 1 suffers 15 cycles of heating and cooling and the total load on the pile raft is 3000 N, which is 52.9% of the ultimate bearing capacity of the pile group. Detailed information of the loading scheme is provided in Figure 9. The focus of this simulation is to reveal the thermal influence on the development of axial force and settlement of non-active piles.
of the pile raft system and the soil foundation keep the same as in Section 4.1, while the pile sequence number is modified according to Figure 2. The main difference is that pile 1, which locates at the diagonal corner, serves as a thermally active pile. Pile 1 suffers 15 cycles of heating and cooling and the total load on the pile raft is 3000 N, which is 52.9% of the ultimate bearing capacity of the pile group. Detailed information of the loading scheme is provided in Figure 9. The focus of this simulation is to reveal the thermal influence on the development of axial force and settlement of non-active piles. Figure 9. Loading scheme. Figure 9. Loading scheme. Figure 10 shows the influence of thermal cycles on the development of pile top settlement. For pile 1, the active pile, a 61% increase of settlement compared to that under only mechanical load can be observed after 15 thermal cycles. During the continuous heating and cooling process, pile body experiences repeated expansion and contraction. This cyclic deformation would lead to the cyclic shearing of surrounding soil, which further induces degradation of shaft resistance and accumulation of plastic deformation. As for pile 2, a non-active pile surrounding pile 1, a 30% increase of settlement can also be observed. The settlement of pile 2 may be due to the constraint effect of pile raft, as piles are all rigidly connected to the pile raft and the settlement of pile 1 would drag pile 2 downwards. As a result, the increase of settlement is not as remarkable as for pile 1. As for pile 7 and 9, which locate at the opposite side of the raft, a reverse trend of settlement development occurs. A notable lift of pile position, around 58% of the initial settlement, occurs for pile 9, which locates at opposite diagonal corner to pile 1. It seems that pile 7 and pile 9 are pulled out by the pile raft, which indicates the rotation movement of pile group along the connecting line of pile 4 and pile 6. The differential settlement caused by thermal loading would pose a great threat to the long-term stability of the superstructure, thus it is necessary to select the location of active piles properly in energy pile group.  Figure 10 shows the influence of thermal cycles on the development of pile top settlement. For pile 1, the active pile, a 61% increase of settlement compared to that under only mechanical load can be observed after 15 thermal cycles. During the continuous heating and cooling process, pile body experiences repeated expansion and contraction. This cyclic deformation would lead to the cyclic shearing of surrounding soil, which further induces degradation of shaft resistance and accumulation of plastic deformation. As for pile 2, a non-active pile surrounding pile 1, a 30% increase of settlement can also be observed. The settlement of pile 2 may be due to the constraint effect of pile raft, as piles are all rigidly connected to the pile raft and the settlement of pile 1 would drag pile 2 downwards. As a result, the increase of settlement is not as remarkable as for pile 1. As for pile 7 and 9, which locate at the opposite side of the raft, a reverse trend of settlement development occurs. A notable lift of pile position, around 58% of the initial settlement, occurs for pile 9, which locates at opposite diagonal corner to pile 1. It seems that pile 7 and pile 9 are pulled out by the pile raft, which indicates the rotation movement of pile group along the connecting line of pile 4 and pile 6. The differential settlement caused by thermal loading would pose a great threat to the long-term stability of the superstructure, thus it is necessary to select the location of active piles properly in energy pile group. The variation of axial force on pile 1 after the first heating is presented in Figure 11. After mechanical loading, the axial force decreases along the pile length. At pile bottom, the base resistance is 77 kN, which is about 22.5% of the total load on the pile head. The distribution of axial force follows the same trend for other piles, as all piles are assumed to be rigidly connected to the pile raft and external load is uniformly applied to the pile raft. After the first heating, the axial force of pile 1 shows an obvious increase, while the axial force of pile 2 decreases slightly, especially at the upper part. It is believed that heat- The variation of axial force on pile 1 after the first heating is presented in Figure 11. After mechanical loading, the axial force decreases along the pile length. At pile bottom, the base resistance is 77 kN, which is about 22.5% of the total load on the pile head. The distribution of axial force follows the same trend for other piles, as all piles are assumed to be rigidly connected to the pile raft and external load is uniformly applied to the pile raft. After the first heating, the axial force of pile 1 shows an obvious increase, while the axial force of pile 2 decreases slightly, especially at the upper part. It is believed that heating pile 1 would lead to expansion of the pile body, which further increases the axial force due to the constraint of pile raft. The upward movement of pile 1 would pull the surrounding non-active pile 2 upwards, as a result, the axial force of pile 2 decreases slightly. For pile 4, which locates on the rotation axis of pile group (as revealed from Figure 10), the influence of heating on the distribution of axial force is minimal. As for pile 9, the rotation of pile raft would lead to the compression of pile body. As a result, the axial force of pile 9 increases, though with smaller magnitude than that of pile 1. To better understand the influence of thermal cycles on the loading performance, variations of the pile head load on pile 1 and pile 2 at different thermal stages are presented in Figure 12. For pile 1, heating will lead to the expansion of pile body, which further induces the increase of pile head load due to the constraint of pile raft, while the following cooling process will cause pile contraction. After the first heating and cooling cycle, the pile head load is smaller than that after mechanical loading. With continuous heating and cooling, load on pile 1 keeps decreasing. And the load variations due to thermal cycles gradually stabilizes at around 4~5 cycles. For pile 2, the load variance shows a reverse trend. Heating on pile 1 will lead to load decrease on pile 2. After each thermal cycle, there will be a slight increase of the pile head load, and load variance also decreases with increasing thermal cycles. After 15 heating and cooling cycles, pile head loads on pile 1 and pile 2 increase by 20% and 8.5%, compared to that under mechanical loading. The continuous heating and cooling process will cause the repeated expansion and To better understand the influence of thermal cycles on the loading performance, variations of the pile head load on pile 1 and pile 2 at different thermal stages are presented in Figure 12. For pile 1, heating will lead to the expansion of pile body, which further induces the increase of pile head load due to the constraint of pile raft, while the following cooling process will cause pile contraction. After the first heating and cooling cycle, the pile head load is smaller than that after mechanical loading. With continuous heating and cooling, load on pile 1 keeps decreasing. And the load variations due to thermal cycles gradually stabilizes at around 4~5 cycles. For pile 2, the load variance shows a reverse trend. Heating on pile 1 will lead to load decrease on pile 2. After each thermal cycle, there will be a slight increase of the pile head load, and load variance also decreases with increasing thermal cycles. After 15 heating and cooling cycles, pile head loads on pile 1 and pile 2 increase by 20% and 8.5%, compared to that under mechanical loading.

Simulation Results
The continuous heating and cooling process will cause the repeated expansion and contraction of the energy pile. Due to the constraint of pile raft, differential settlement may occur for piles at different location, which would further lead to the stress redistribution of each pile and influence the long-term bearing performance of energy pile group. According to the above analysis, the effect of thermal loading on the pile group behavior is very complex and depends on pile spacing, active pile number, active pile locations, etc. The proposed analytical approach is expected to provide a preliminary evaluation of the energy group design before sophisticated numerical modelling and experimental research. cycle, the pile head load is smaller than that after mechanical loading. With continuous heating and cooling, load on pile 1 keeps decreasing. And the load variations due to thermal cycles gradually stabilizes at around 4~5 cycles. For pile 2, the load variance shows a reverse trend. Heating on pile 1 will lead to load decrease on pile 2. After each thermal cycle, there will be a slight increase of the pile head load, and load variance also decreases with increasing thermal cycles. After 15 heating and cooling cycles, pile head loads on pile 1 and pile 2 increase by 20% and 8.5%, compared to that under mechanical loading. The continuous heating and cooling process will cause the repeated expansion and contraction of the energy pile. Due to the constraint of pile raft, differential settlement may occur for piles at different location, which would further lead to the stress redistribution of each pile and influence the long-term bearing performance of energy pile group. According to the above analysis, the effect of thermal loading on the pile group behavior is very complex and depends on pile spacing, active pile number, active pile locations, etc. The proposed analytical approach is expected to provide a preliminary evaluation of the energy group design before sophisticated numerical modelling and experimental research.

Conclusions
(1) In this paper, a simple and convenient analytical approach is proposed to evaluate the thermal effect on the bearing characteristics of energy group pile. The nonlinear pilesoil interaction is simulated based on an exponential model and the axial load transfer (τ-z) method is used to calculated the distribution of shear stress. The thermomechanical behavior is simplified as mechanical and thermal load.
(2) The proposed method shows good prediction of axial force distribution for pile raft foundations with small pile spacing. When pile spacing increases, the influence of raft-soil interaction needs to be considered.
(3) For energy pile group with thermally-active piles and non-active piles, the heat exchange of the energy pile may cause differential settlement and redistribution of axial force of the pile group. This thermal effect is highly related to pile spacing, energy pile number, and energy pile location.
(4) For energy pile located at the corner of the pile group, thermal cycles would lead to the rotation of the whole structure. During heating, load on the head of energy pile would increase due to the interaction between pile and raft. With continuous heating and cooling, pile head settlement gradually accumulates and load decreases. For non-active piles in the pile group, the settlement and load behavior may differ at different location.  Data Availability Statement: All the data used in this research can be easily accessible by downloading the various documents appropriately cited in the paper.