This study utilizes the discrete element method (DEM) to explore the thermal–mechanical coupling behavior of slate. Initially, a series of triaxial compression tests are conducted on slate specimens with varying orientation angles, temperatures, and confining pressures. Subsequently, the DEM is employed to replicate the mechanical and failure characteristics of the slate under identical laboratory conditions. The DEM constitutive models incorporate the parallel bond model for the rock matrix and the smooth-joint model for foliation. A thermal-degradation function for the parallel bond model is proposed, reasonably capturing the temperature's effect on the strength and failure characteristics of the slate. The simulation results exhibit a series of U-shaped curves in peak strength at different orientation angles, mirroring experimental findings. The study's findings indicate that orientation angles significantly influence the mechanical behavior and failure characteristics of slate. Elevated temperatures have minimal effects on strength anisotropy. Furthermore, applying high confining pressure is observed to mitigate the impact of high temperatures on bond weakening. This proposed method provides a useful tool for tunnel designs and geothermal energy development.