A novel theory which employs the K-BKZ viscoelastic integral constitutive model was developed in this study for simulating the non-isothermal injection molding filling process and the frozen stresses of a three-dimensional thin part. The simulation of viscoelastic model in such a problem has not yet been discussed in most previous works as a result of the complexity of this problem. A quasi-steady approximation concept was developed in the present investigation for solving the non-isothermal filling process via K-BKZ viscoelastic integral constitutive model. Additionally, the numerical method of flow field was based on the control volume finite element method. The flow field of generalized Newtonian fluid was used as the initial guess of flow kinematics. The quasi-steady state approximation was introduced for each element to calculate the flow kinematics and stress profile from the K-BKZ integral constitutive model. The finite difference method with streamline upwind characteristic was adopted here for calculating the temperature field of process. When the cavity is fully filled, the subsequent non-isothermal stresses would relax after cessation of flow. Thereby, the frozen stresses (or frozen birefringence) could be obtained.