The challenges of computational cost and robustness are critical obstacles in topology optimization methods, particularly for the iterative process of optimizing large-scale multiphysical structures. This study proposes an efficient and robust topology optimization method for minimizing the thermoelastic damping of large-scale microresonators. An evolutionary structural optimization method is adopted to passively determine the search direction of optimizing large-scale thermoelastic structures. To efficiently reduce the computational cost of the iterative process of an optimizing process, a model reduction method is developed based on the projection-based model reduction method whose reduced basis is generated within the Neumann series subspace. However, the projection-based model reduction method may be unstable when topology modifications are made during an iteration optimization process. To ensure robustness, a modal validation technique is first implemented in the iterative process to stabilize the iteration and narrow down the search domain, and a posterior evaluation of the Neumann series expansion is then developed to retain the convergence of the projection-based model reduction method. Furthermore, the efficiency and accuracy of the proposed topology optimization method are validated through numerical examples. Two large-scale numerical models are also used to demonstrate the advantage of the proposed method. It is found that large-scale thermoelastic structures with a phase-lag heat conduction law can be designed passively, precisely, and efficiently by using the proposed topology optimization method.