An effective approach is presented for the numerical solution of the equations governing steady laminar and turbulent flow, heat and mass transfer at low Mach number. The approach adopted combines a compact and accurate discretization using the residual-distribution (RD) approach with a fully-coupled implicit solution procedure. The system RD approach adopted employs genuinely-multidimensional upwinding to achieve accurate and stable discrete equations on a highly compact computational stencil. This combines very naturally with the fully-implicit coupled solution procedure for which the number of non-linear iterations required is essentially independent of the grid size. This contrasts with other widely used segregated approaches in which the pressure-velocity system is discretized and solved as a set of scalar equations. A further key distinction from many other implicit methods is that the compact nature of the discretization allows the full convection and diffusion terms in all equations to be treated implicitly without any form of deferred correction. The code developed solves the 2D, axisymmetric and 3D Navier-Stokes equations in incompressible or weakly-compressible form on unstructured grids of triangles or tetrahedra. The RD form of the system Lax-Wendroff scheme is applied to the convection and pressure terms while the viscous terms are treated using the Galerkin finite-element method. The system scheme for convection provides natural stabilization, allowing a collocated variable arrangement to be used. Convection terms in scalar equations are treated using scalar multidimensional RD schemes such as the N and PSI schemes, both of which are positive. The discrete equation system is solved using a fully-coupled implicit approach based on Picard/Newton linearization with the linear system solved using standard Krylov subspace methods (e.g. GMRES or BiCGStab) with ILU(0) preconditioning. A number of practical issues relating to the solution procedure have been investigated including parallelization and equation segregation. A domain decomposition method has been developed for the treatment of large problems using a multiplicative-Schwarz approach with arbitrary inter-block overlap. The approach has been extensively validated on a number of 2D, axisymmetric and 3D test cases, with and without heat and mass transfer. This has included direct comparisons with commercial unstructured flow solvers with very promising results - showing equivalent levels of accuracy but reduced computational times and less sensitivity to grid quality. Examples of the validation and application of the method are presented