The aim of this article is to describe the development of a mathematical model for the simulation of the flow of polymer melts inside internal mixers. The rheological behavior of the polymeric fluid is assumed to be described by the Carreau equation. The flow regime is considered to be non-isothermal. The set of the governing equations are solved using the finite element method for both steady-state and transient conditions. In the steady state case, the flow equations are solved by the penalty method using the standard Galerkin technique. Petrov-Galerkin schemes based on both consistent and inconsistent streamline upwinding methods are employed to solve the energy equation and the obtained results are compared. Transient velocity, pressure, and stress fields are modeled using implicit theta method. Ln addition to implicit theta method, two versions of the Taylor-Galerkin approach are used to solve the transient energy equation. Slip-stick on the solid walls, encountered in the flow of viscous fluids, is incorporated in the model by the use of Navier＇s slip conditions. We describe two new methods for the inclusion of this condition in the working equation. Our simulations yield the velocity field, distribution of pressure, stress and temperature in the steady state and the variations of these parameters with respect to time under transient conditions. As an example of the applicability of the developed model, a typical mixing problem which involves convection of carbon black with flowing rubber in a domain representative of the section under the blade of a tangential rotor mixer is simulated. Concentration profiles of carbon black in the rubber matrix in this case is obtained by the solution of carbon mass continuity equation in conjunction with the flow model. This solution gives the distribution of filler volume fraction at different mixing times in the mixer. Comparison of the obtained results with the available experimental data gives some indication of the validity of the model.