Finite element analyses of energy piles using different constitutive models

. Energy piles are foundation elements having the double scope of transferring structural loads from the structure to the ground and of exchanging heat with the surrounding soil. It follows that pile state of stress and settlement are altered by the time-dependent temperature change in both pile and soil. This work is aimed at investigating the effect of thermal cycles on the behaviour of a single energy pile. To this end, fully coupled thermo-hydro-mechanical analyses have been carried out using the Finite Element code ABAQUS. The single pile is installed in a normally consolidated clay behaving according to different constitutive models involving Mohr-Coulomb, Modified Cam Clay and Hypoplastic. The latter is employed with and without the thermal formulation capable of accounting for the thermal collapse of NC clays during heating. A single free-head pile is considered and the results are presented in terms of pile axial force and settlement developed cycle by cycle.


Introduction
Results of laboratory tests on soil specimens have proven that normally consolidated and slightly overconsolidated clays experience volumetric contraction when heated in drained conditions [1][2][3][4][5][6][7], and this certainly affects the behaviour of a pile embedded in this soil. The thermallyinduced collapse stabilizes after few thermal cycles when new contacts between the particles are formed. Centrifuge and small scale tests showed that settlements are accumulated cycle after cycle [8,9].
To investigate the effects of the thermal cycles on the global performance of a single energy pile embedded in NC clay and subjected to mechanical and cyclic thermal loads, fully coupled thermo-hydro-mechanical analyses have been carried out.
In the FE analyses, the soil behaviour is modelled employing different constitutive models, i.e. the Mohr-Coulomb (MC), the Modified Cam-Clay (MCC), and the Hypoplastic model with and without the thermal formulation (Hypo and Hypo-T, respectively) The classical elasto-plastic models are developed in the framework of the perfect or hardening plasticity theory and are not capable of properly describing the accumulation of irreversible strain and pore water pressures caused by the cyclic thermal loading. In recent years, in order to improve the predictive capabilities of the conventional constitutive models, new approaches have been developed, among which the Hypoplastic theory with internal state variables formulated as an extension of the Hypoelastic theory [10,11]. Hypoplastic model may also include the thermal formulation capable of reproducing the thermal collapse of the NC clay [12,13].

Finite element analyses
The fully coupled thermo-hydro-mechanical analyses have been performed via the FE code ABAQUS. The single energy pile has a slenderness ratio L/d = 50, where L and d are the pile length and diameter, respectively. Table 1 reports a summary of the geometrical, mechanical and thermal data of the pile.
The axial symmetry of the problem allows to consider a 2D model having dimensions of 50x50 m. As an outcome of a sensitivity study, the elements used have a variable geometry. In particular, a structured mesh composed of 4-node axisymmetric quadrilateral elements strongly refined at the interface (where the dimensions are the same as for the elements of the pile, i.e. base of 0.05 m and a height of 0.125 m) and in the vicinity of the pile tip, has been used (Fig. 1). The pile is a non-porous elastic medium; the elements employed are of the type CAX4T with displacement and temperature degrees of freedom. The soil is modelled with elements of the type CAX4PT which also include the pore pressure degree of freedom. The total number of elements is 1000 and 27956 for pile and soil, respectively. At the interface perfect contact between pile and soil is modelled with the use of tie constraints. This assumption is justified considering that for a pile in clay the degradation phenomenon due to the thermal cycles is of minor concern and the increase of temperature improves the interface strength [14]. The groundwater table is located at the ground surface; the soil permeability, kw, is set to a value of 10 -10 m/s while the saturated soil unit weight, γsat, is equal to 20 kN/m 3 . The thermal properties used for the thermal characterization of the soil are reported in Table 1. It is worth noting that, in order to reduce the computational efforts, the thermal expansion coefficient of the water, αw, is taken constant throughout the analyses and equal to the value at the initial temperature. The first step of the analysis involves a geostatic initialization. The loading phase starts with a mechanical step in which the mechanical load is applied at pile head in drained conditions, thus resulting in no excess pore water pressures at the end of the loading stage. The load level, Q/Rt, is equal to 30%, Q being the applied load and Rt = 1490 kN the bearing capacity of the pile evaluated with the β method, considering the pile as wished in place. Then, the temperature is varied according to the path shown in figure 2 for a total of 5 thermal cycles (i.e. 5 years of operation). A 1 month of recovery (i.e. a rest phase) between the heating and the cooling is also simulated, in which the temperature is not imposed as a boundary condition to the pile, whose temperature thereby varies as a result of the thermal equilibrium. Since the scope of the present work is the evaluation of the effect of the thermal cycles on the pile behaviour, the circulation of the heat-carrier fluid has not been modelled and the temperature variation has been uniformly applied on all the elements of the pile. Past works have shown that the choice of a uniform temperature applied to each pile cross section and along its length does not significantly affect the results [15,16].

Calibration of the constitutive models
The Hypoplastic model is an advanced incrementally non-linear constitutive model developed on the basis of the Critical State Soil Mechanics. In the present analyses, the basic version enhanced for the prediction of the small strain stiffness (i.e. with the intergranular strain concept [17]) has been used. In particular, a robust calibration of the mechanical parameters is available in Mašín [11] for London clay ( Table 2). It has to be noted that since London clay is overconsolidated, all the relevant parameters have been adjusted in order to simulate the normally consolidated state (for instance Ag and ng have been selected following Viggiani and Atkinson [18]). Regarding the MC and the MCC, the mechanical parameters have been chosen using the Hypoplastic model as a benchmark. In this way the results of the FE * φ is the friction angle; λ*, κ*, N* are the parameters of the Critical State Soil Mechanics defined in the ln(1+e) -lnp' plane; ν is the Poisson's ratio; Ag and ng control the magnitude of the small strain stiffness G0: G0 = Ag (p') ng ; R is the size of the elastic range; βr regulates the rate of evolution of the intergranular strain tensor; χ controls the interpolation between the reversible elastic and non-linear hypoplastic responses; mrat represents the ratio G90/G0, being G90 the stiffness after a 90° change in the strain path; ϑ regulates the accumulation of strains or stresses in cyclic loading paths.
analyses in terms of axial force and settlements can be compared. In particular, for the MC model the soil stiffness used has been selected in order to obtain the same settlement as for the Hypoplastic model at a load level of Q/Rt ∼ 40% resulting in G0,MC = G0/2.1; for the MCC the parameters have been derived from the Hypoplastic ones considering that the Critical State theory is incorporated in both models with the difference that the MCC is developed with reference to the (1+e)lnp' plane while the Hypoplastic is developed with reference to the ln(1+e)-lnp' plane. It is worth noting that, since the elastic part is only defined by the bulk stiffness and the Poisson's ratio, the MCC model cannot reproduce the small-strain response of the Hypoplastic model. A summary of the selected parameters for the classical models is reported in Table 3. Concerning the thermal part of the thermo-Hypoplastic model, to account for the accumulation of non-recoverable volumetric contraction due to thermal collapse of the NC clay, Mašín and Khalili [12] introduced the thermal stabilization line, TSL, a straight line in the ln(1+e) -lnp' plane parallel to the normal consolidation line, NCL. According to the model, the soil attains a stable structure (i.e. zero volumetric contraction) when its status lies on the TSL. Since no tests were available for the thermal characterization of London clay, the thermal parameters were estimated following Ma et al. [13] (Table 4).  While the mathematical formulation of the MC and the MCC models is already available in the ABAQUS library, an external subroutine (UMAT) with the Hypoplastic model written in Fortran language has been used (accessible from www.soilmodels.com). The thermal part has been implemented by modifying the original algorithm.

Results and discussion
The temperature variation imposed to the pile ( fig. 2) leads to a thermal field that is independent of the constitutive model chosen to simulate the soil behaviour. Therefore, in all the cases the pile is subjected to expansion at increasing temperature and to contraction when the temperature decreases.
After the application of the mechanical load the pile is heated up to 29°C. Since the maximum temperature is applied in one day, only the soil in the vicinity of the pile is affected by the heat conduction while the remaining part is undisturbed. In this situation, the pile free expansion is partly prevented by the surrounding soil for two main reasons: first, the temperature distribution in the soil remains practically constant and equal to the initial temperature; further, the thermal expansion coefficient of the pile and the soil are different. As a consequence, a thermally-induced axial compression (positive value in the following figures) develops along the pile shaft and at the base, that sums up with the mechanical axial load. The pile head moves upward recovering part of the settlement (downward movements are taken as positive) experienced after the application of the mechanical load. During the 5 months in which heating is kept constant, owing to heat conduction, the temperature increase involves a wider portion of the soil leading to a decrease of the thermally-induced axial force, whereas the pile head moves upward; as it will be shown later, except for the Hypoplastic model with the thermal formulation, this is due to the thermal expansion of the soil. In the following step, no temperature condition is imposed to the pile and a thermal equilibrium is reached according to the boundary conditions of the model; the temperature of both pile and soil reduces determining the reduction of the axial force and the pile settlement. When the temperature is reversed, the upper part of the pile moves downward while the lower part goes upward; the soil opposes to pile shortening and the generated tensile stress determines the reduction of the axial force. As in the continuous heating phase, during the continuous cooling the temperature decrease involves further portion of the soil and, therefore, the tensile thermally-induced axial force reduces, causing the increase of the total axial force, whereas downward movements take place. A further increase of the axial force and decrease of the head settlement is encountered in the following thermal rest phase, in which the pile temperature, variable with depth, grows up.
In figure 3 the cyclic axial force distribution is shown for all the models. Since the maximum and minimum thermally-induced axial forces are experienced at the end of the heating and cooling phases, respectively, only the results of this two phases are reported. The results in terms of cyclic settlements at pile head are presented in figure 4 with reference to all the steps. It is worth noting that the final configuration in terms of temperature Concerning the MC model the highest thermallyinduced forces are equal to about 9%Rt and -10%Rt for heating and cooling respectively; the maximum overall axial force (32%Rt) is located along the shaft. As a consequence of the different thermal field, a slight difference is found in the axial loads distributions from the first to the second cycle. At the end of the mechanical loading phase the pile displacement is equal to 0.57%d. The maximum and minimum values of the thermally-induced settlement are equal to ± 0.5%d and are attained during the continuous cooling and continuous heating, respectively. It is possible to notice that, after the first cycle, there is no accumulation of irreversible displacements.
Compared to the other models, the calibration of the MCC model resulted in a lower initial stiffness profile and thereby in a larger pile-to-soil stiffness ratio. This leads to a higher value of the axial force along the shaft due to the mechanical load and to a lower value of the thermally-induced axial force. For the same reason, settlements are larger in magnitude, but the global response of the system remains reversible. The highest thermally-induced forces are equal to about 5%Rt and -5%Rt for heating and cooling, respectively; the increase in axial force along the shaft for increasing temperature causes a maximum axial force equal to 31%Rt. Similarly to the MC, the slight difference between the load distributions in the first and second cycles is due to the different thermal field at the beginning of the cycle. At the end of the mechanical loading phase the pile displacement is equal to about 1.4%d. The maximum and minimum thermally-induced settlement develop at the end of the continuous cooling and continuous heating phases, their values being equal to about ± 0.54%d. As expected, no accumulation of irreversible settlements is found.
For the Hypoplastic and the thermo-Hypoplastic models, the largest and lowest thermally-induced forces are observed during heating and cooling of the first cycle and are equal to about 6%Rt and -4%Rt, respectively. The axial force envelope moves rightward cycle after cycle at decreasing rate. Similar results are obtained by Suryatriyastuti et al. [19]. After 5 thermal cycles, the maximum overall force is located along the shaft and is equal to 36%Rt (Hypoplastic) and 34%Rt (thermo-Hypoplastic), respectively. It can be noticed that thermally-induced load is not affected by the activation of the thermal term. Although this term allows to account for the volumetric collapse of the soil as a consequence of the temperature increase, the undrained conditions assumed in the analyses prevent the volumetric collapse and cause further increase in the excess pore pressures. While in models without the thermal formulation both soil and water would expand for a positive temperature variation and the pore pressures change as function of the difference between the thermal expansion coefficient of water and soil, in the case of the thermo-Hypoplastic model the soil would contract, and this contraction is contrasted by the expansion of the water, which must thus adsorb a larger pressure to guarantee the compatibility. The similar volumetric behaviour of the Hypoplastic and thermo-Hypoplastic models is responsible of the similar axial force distribution along the pile.
Concerning the displacements, the activation of the thermal term is responsible for contraction during the continuous heating phase (partially drainage occurs in this phase); therefore, differently from the other constitutive models, the pile experiences additional downward movements. The overall behaviour is characterized by irreversible accumulation of cyclic settlements that, after 5 thermal cycles, are equal to 1.97%d and to 3.75%d for the Hypoplastic and the thermo-Hypoplastic models, respectively. With respect to the mechanical loading phase (0.47%d for both models) the settlement is about 4.2 (Hypoplastic) and 7.9 (thermo-Hypoplastic) times larger. The accumulation slightly decreases cycle after cycle; hence, the highest irreversible displacement occurs during the first cycle, being about 1.75 (Hypoplastic) and 2 (thermo-Hypoplastic) times larger than that experienced at the end of the mechanical loading phase. Similar results are obtained by Wu et al. [9], who performed a series of tests on small scale floating energy piles subjected to mechanical load (around the 40% of the bearing capacity) and to 5 thermal cycles (14°C and -13°C). The piles were installed in normally consolidated saturated clays. After the application of the mechanical load plus 5 thermal cycles, the free pile continues to settle showing a similar trend as in figure 4. The overall displacement due to the temperature variation was equal to about 1.5%d. Other studies on the response of the single pile in a normally consolidated soil also showed similar trend. For instance, for a lightly overconsolidated clay (OCR=1.7), Ng et al. [8] showed that, after performing 5 cycles (ΔT=13°C and ΔT=-10°C) the overall displacement was in the order of 3.8%d, with the settlement due to the application of the mechanical load (about 40% of the pile bearing capacity) equal to 0.67%d. Likewise, the increase of irreversible settlement cycle after cycle as a function of the load level is reported by Vieira and Maranha [20].

Conclusions
The scope of the study presented herein is to investigate the cyclic performance of a single energy pile embedded in NC clay. Fully coupled thermo-hydro-mechanical finite element analyses have been performed employing constitutive models of different complexity. Two elastoplastic models available in the common FE codes along with an advanced constitutive model based on the hypoplasticity theory has been employed. The thermal collapse of the NC clay during heating has been accounted enhancing the advanced model with a thermal term. The single energy pile is first subjected to an axial mechanical load and then to 5 thermal cycles.
The main results of this study are that (i) the use of advanced models like the Hypoplastic model is the only option to qualitatively reproduce the cyclic accumulation of irreversible settlements; (ii) as concerns the axial load distribution, differences in the predictions by different models are found.