A two-dimensional finite element model was developed to investigate thermoelastic instability in multi-layered friction discs with finite thickness, considering deformation modes of steel core. The model was used to simulate four unstable modes that can occur during the engagement process, and the Fourier reduction was applied to calculate the change in critical speed under these modes. Additionally, the influence of thermal physical parameters, including the elastic modulus, thermal expansion coefficient, Poisson's ratio and thermal conductivity of the friction pair, on thermoelastic instability was examined. The findings indicate that the critical speed of the friction pair is lower under the symmetric (friction disc)-antisymmetric (steel disc) mode compared to the other three modes. Consequently, the symmetric-antisymmetric mode is the first to be excited and serves as the dominant mode during thermoelastic instability. Moreover, there exists a specific wave number at which the system exhibits the lowest critical speed and poorest stability. Enhancing the thermal conductivity of the friction disc and steel disc, as well as reducing the thermal expansion coefficient of the steel disc and the elastic modulus and Poisson ratio of both discs, can improve the thermoelastic stability of the friction pair. Notably, the thermal expansion coefficient of the friction disc has minimal impact on thermoelastic instability. These results provide a theoretical foundation for exploring the relationship between the thermal failure of friction pairs and rotational speed, as well as optimizing overall performance design.