The major advancement in the unstructured grid capabilities provided a promise for easier and faster grid generation on complex geometries. However, the need for the gradients has limited the promise of unstructured grid technology because the gradients are susceptible to numerical noise and exhibit significant irregularity and oscillations on non-hexahedral unstructured grids. This project addresses a long-standing issue in predicting accurate and smooth solution gradients for purely tetrahedral elements. Today's state of the art is to completely bypass the problem by reverting to structured-type layers within the boundary layer (BL), and also limiting grid adaptation within the BL and possibly around discontinuities (shocks). The noisy and inaccurate reconstructed solution gradients within the computational domain are still inevitable and used in adjoint error estimation, turbulence modeling, transition prediction, etc. These significantly limit the unstructured grid technology that was once promised to make CFD calculations easier than structured grid methods. The novelty of this proposed work lies in: (1) hyperbolic formulation of the Navier-Stokes (NS) equations, where the gradients are solved together with the flow variables, and therefore not reconstructed, and (2) construction of Newton's method with the exact Jacobian using a second-order Residual-Distribution (RD) scheme. The idea of the hyperbolic RD NS solver is unique and truly beyond currently available approaches. First, the NS equations are reformulated to include the gradients as additional variables, so that the gradients are computed directly, not by reconstruction from the flow variables. Second, the reformulated NS system has a unified characterization, i.e., hyperbolic for all terms, whereas the original NS system is of mixed nature, having hyperbolic and parabolic terms. Numerical techniques have been well developed for the hyperbolic terms, but there remain technical difficulties for the parabolic terms, which govern, for example, the heat fluxes and the viscous stresses. Also, the parabolic terms introduce stronger stiffness than hyperbolic terms and thus, impede construction of an efficient solver. Elimination of the parabolic nature by the proposed method, therefore, brings a dramatic simplification in the construction of numerical schemes, which consequently leads to drastic improvements in accuracy as well as a tremendous advantage for the solver efficiency. The RD scheme has a very attractive feature of being second-order accurate on a compact stencil for hyperbolic terms. The feature is lost for the parabolic terms, and it is one of the reasons that the scheme has not been widely employed in practical codes. But it is now possible to keep the attractive feature of the RD scheme for the full NS systems by making all terms hyperbolic. The numerical code under development in the proposed work will have the following distinguished capabilities: (1) accurate gradient predictions on highly-stretched and skewed unstructured grids, (2) allowing highly arbitrary grid generation and adaptation, (3) equal order of accuracy for both solutions and gradients on irregular grids, and (4) compact second-order discretization.