Cooling system design in glass bulb pressing operation can greatly affect the productivity and the quality of the final product. The concept of cyclic-averaged steady temperature field is proposed in modeling. Heat transfer in the mold region is considered to be a cyclic-steady, three-dimensional conduction; heat transfer within the glass melt region is treated as a transient, one-dimensional conduction; heat exchange between the cooling system surface and coolant is treated as a steady heat convection. A hybrid model consisting of a three-dimensional boundary element method for the mold region and a finite-difference method with a variable mesh for the melt region is used for numerical simulation. Compared with the experimental data, the numerical model developed here is computationally efficient and sufficiently accurate.