Abstract: The prediction of hypersonic aerodynamic environment is highly sensitive to numerical algorithms and computational grids. Due to the geometry complexity of modern lifting-body hypersonic vehicles, the time cost of generating high quality structured grids increases horribly, and it is difficult to meet the need of engineering applications. In order to shorten the task cycle, it is necessary to develop numerical methods for accurate aero-heating prediction on hybrid unstructured grids, due to the capability for complex configurations. Gradient reconstruction method is one of the important factors influencing the accuracy of heat flux prediction on unstructured grids. In this work, a multi-dimensional gradient reconstruction was introduced into the second-order cell-centered unstructured finite-volume solver developed by the authors, and compared with the widely used Green-Gauss methods and least squares methods in aerodynamic force simulations. The benchmark cases of hypersonic flows over a 2D cylinder, were tested to validate the multi-dimensional approach. The numerical results demonstrate that the multi-dimensional approach can improve the accuracy of aero-heating prediction and the convergence performance. Finally, the multidimensional gradient reconstruction method is applied to predict the heat flux over a 3D cylinder and the double ellipsoid configuration on the general hybrid grids. The numerical results agree well with experimental data, which demonstrates the potential capability for complex engineering applications.