This study investigates the interaction between energy piles and layered saturated soils, considering the consolidation induced by the thermal loads and mechanical loads. Initially, the coupled thermo-hydromechanical solution of layered media is obtained by utilizing the boundary element method (BEM) and the transformed differential quadrature method. Subsequently, the energy piles are discretized and modelled by the finite element method (FEM), and the solving equation for piles is established. To reflect the interaction between piles and soils, a coupled BEM-FEM matrix equation is formulated and solved by incorporating displacement coordination conditions and force equilibrium conditions. This approach facilitates the analysis of the temporal evolution of displacements and temperatures of piles and surrounding soils. The proposed methodology is validated through comparisons with monitoring data of field tests and results from simulations. Ultimately, the key factors, including the temperature increments, mechanical loads, length-diameter aspect ratio are examined through examples.