A general description of a thermally coupled fluid flow is given by the incompressible Navier-Stokes equations coupled with the heat equation using Boussinesq approximation, whose mathematical structure is much well understood. A variational multiscale finite element approximation has been considered for the formulation of incompressible Navier-Stokes equation and heat equation. The complexity of these problems makes their numerical solution very difficult as the standard finite element method is unstable. In the incompressible Navier Stokes equations, two well known sources of numerical instabilities are the incompressibility constraint and the presence of the convective term. Many stabilization techniques used nowadays are based on scale separation, splitting the unknown into a coarse part induced by the discretization of the domain and a fine subgrid part. The modeling of the subgrid scale and its influence leads to a modified coarse scale problem providing stability. In convection-diffusion problem once global instabilities have been overcome by a stabilization method, there are still local oscillations near layers due to the lack of monotonicity of the method. Shock capturing techniques are often employed to deal with them. Proper choice of stabilization and shock capturing techniques can eliminate the local instabilities near layers of convection-diffusion equation. A very important issue of the formulation presented in this thermally coupled incompressible flow is the possibility to model turbulent flows. Some terms involving the velocity subgrid scale arise from the convective term in the Navier-Stokes equations which can be understood as the contribution from the Reynolds tensor of a LES approach and the contribution from the cross stress tensor. This opens the door of modeling thermal turbulence using LES automatically inherited by the formulation used in this work. Different classical benchmark problems are numerically solved in this thesis work for the convection-diffusion equation to show the capabilities of different combination of stabilization and shock capturing methods. In the case of thermally coupled incompressible flows some numerical and industrial examples are exhibited to check the performance of the different combination of stabilization and shock capturing methods and to compare them. The objective is to conclude which method works better to approximate the exact solution and eliminate instabilities and local oscillations.