Boundary element method (BEM) is widely used in engineering analysis, especially in solving the transient heat conduction problem because of the advantage that only boundary of the problem needs to be discretized into elements. The general procedure of solving the variable-coefficient transient heat conduction problem by using the BEM is as follows. First, the governing differential equations are transformed into the boundary-domain integral equations by adopting the basic solution of the linear and homogeneous heat conduction problemGreen function. Second, domain integrals in the integral equation are converted into boundary integrals by the radial integral method or the dual reciprocity method. Finally, the time difference propulsion technology is used to solve the discrete time differential equations. A large number of practical examples verify the correctness and validity of the BEM in solving the variable coefficient of transient heat conduction problem. However, two deficiencies are encountered when the system of time differential equations is solved with the time difference method, i.e., one is the stability of the algorithm, which is closely related to the time step size, and the other is time-consuming when the freedom degree of the problem is large and all specified time steps are considered, because a system of linear equations needs to be solved in each time step. Therefore, in this paper we present a reduced order model analysis method of solving the variable-coefficient transient heat conduction problem based on BEM by using the model reduction method of proper orthogonal decomposition (POD). For variable-coefficient transient heat conduction problems, the discrete integral equations which are suitable for order reduction operation are deduced by using the BEM, the reduced order model is established by using the model reduction method of POD, and a lowdimensional approximate description of the transient heat conduction problem under time-varying boundary condition is obtained by projection of the initial discrete integral equations on some few dominant POD modes obtained from the problem under constant boundary conditions. First, for a variable coefficient transient heat conduction problem, boundary-domain integral equations are established and the domain integrals are transformed into boundary integrals by using the radial integration method. Second, the time differential equations with discrete format which is suitable for order reduction operation are obtained by reorganizing the integral equations. Third, the POD modes are developed by calculating the eigenvectors of an autocorrelation matrix composed of snapshots which are clustered by the given results obtained from experiments, BEM or other numerical methods for transient heat transfer problem with constant boundary conditions. Finally, the reduced order model is established and solved by projecting the time differential equations on reduced POD modes. Examples show that the method developed in this paper is correct and effective. It is shown that 1) the low order POD modes determined under constant boundary conditions can be used to accurately analyze the temperature field of transient heat conduction problems with the same geometric domain but a variety of smooth and time-varying boundary conditions; 2) the establishment of low order model solves the problem of heavy workload encountered in BEM where a set of large linear equations will be formed and solved in each time step when using the time difference method to solve the large time differential equations.