We present an extension of multidimensional upwind residual distribution schemes to viscous flows. Following (Ricchiuto et al. , J.Comp.Appl.Math. 2007), we consider the consistent coupling of a residual distribution (RD) discretization of the advection operator with a Galerkin approximation for the second order derivatives. Consistency is intended in the sense of uniform accuracy with respect to variations of the mesh size or, equivalently, for the advection diffusion equation, of the Peclet number. Starting from the scalar formulation given in (Ricchiuto et al. , J.Comp.Appl.Math. 2007), we perform an accuracy and stability analysis to justify and extend the approach to the time-dependent case. The theoretical predictions are cofirmed by numerical grid convergence studies. The schemes are formally extended to the system of laminar Navier-Stokes equations, and compared to more classical finite volume discretizations on the solution of standard test problems.