Flow induced residual stress is the major reason for stress cracking and warping of plastic products, the study on which is significant to overcome the flaws of products. In this paper, the energy equation based on Extended Pom-Pom constitutive relationship is deduced. A non-isothermal viscoelastic-Newtonian two-phase fluid model for mold filling process of viscoelastic materials is set up. The conservative interface capturing technique and the flow field solving method are coupled to perform a dynamic simulation. The distribution of the frozen skin layer and the shear rates are given. The flow induced residual stress is predicted and analyzed. The numerical results show that the thickness of the frozen skin layer is dependent on the injection velocity and a higher injection velocity corresponds to a thin frozen skin layer. Near the walls of the product, the shear rate and the residual stress are almost zero. At the position of subsurface, the shear rate and the residual stress reach their largest values. At the positions far away from the walls of the product, the shear rate and the residual stress are small.