This paper deals with the computational modeling and numerical simulation of the material flow around the probe tool in a Friction Stir Welding (FSW) process. Within the paradigmatic framework of the multiscale stabilization methods, suitable subgrid scale stabilized coupled thermomechanical formulations have been developed using an Eulerian description. Norton-Hoff and Sheppard-Wright thermo-rigid-viscoplastic constitutive material models have been considered. Constitutive equations for the subgrid scale models have been proposed and an approximation of the subgrid scale variables has been given. In particular Algebraic Subgrid Scale (ASGS) and Orthogonal Subgrid Scale (OSGS) methods for P1/P1/P1 linear elements have been considered. Furthermore, it has been shown that well known classical stabilized formulations, such as the Galerkin Least-Squares (GLS) or Streamline Upwind/Petrov-Galerkin (SUPG) methods, can be recovered as particular cases of the multiscale stabilization framework considered. Within the framework of a product formula algorithm, the resulting algebraic system of equations has been solved using a staggered procedure, in which a mechanical problem, defined by the plastic strain rate incompressibility equation and the quasi-static linear momentum balance equation, is solved at a constant temperature and a thermal problem, defined by the energy balance equation, is solved keeping constant the mechanical variables. The computational model has been implemented in the in-house developed FE software COMET. An assessment of the influence of the thermal deformation in the formulation has been carried out. Results obtained show that the influence of the thermal deformation is very small and can be neglected, getting a fully incompressible formulation. Finally, the computational model implemented in COMET has been validated through a number of examples, including a 3D numerical simulation of a FSW process. Numerical results obtained have been compared with experimental results available in the literature. A good agreement on the temperature distribution has been obtained and predicted peak temperatures compare well, both in value and position, with the experimental results available.