Modern civil jet engines arrange components on spools with different rotational speeds in order to improve compressor stall margin, overall engine performance, etc. The unsteady interactions among these components can be significant and should be considered at an early design stage if possible. Unsteady Reynolds-averaged Navier–Stokes (URANS) is a common approach to simulate these unsteady effects, but the disparity in time scales in a multispool simulation can lead to expensive URANS simulations. Harmonic methods are effective and efficient approaches to simulate unsteady interactions among turbomachinery components, but their applications to multispool simulations remain a challenge. The objective of this paper is to address this challenge. This paper extends the Favre-averaged non-linear harmonic method to simulate multispool turbomachinery components using a unified bladerow interface which transfers disturbances through bladerows with arbitrary blade counts at any rotational speed. The regularization of non-reflective boundary condition is described for certain circumferential wave number of the zero-frequency mode. The capability of the proposed approach is demonstrated by simulating the transfer of hot streaks through full 3D high- and intermediate-pressure turbines in a three-shaft engine. The temperature distributions from the harmonic method show good agreement with direct unsteady simulation in terms of the mean flow and the instantaneous flow. The radial migration of the hot streaks towards the hub are captured very well by the proposed harmonic method. The required wall-clock time of the harmonic method is roughly 240 times smaller than the whole annulus URANS simulation. This demonstrates that the proposed method can be an efficient design tool to trace hot streaks in multispool turbines at the early design stage.