A numerical method is proposed to approach the Approximate Inertial Manifolds (AIMs) in unsteady incompressible Navier-Stokes equations, using multilevel finite element method with hierarchical basis functions. Following AIMS, the unknown variables, velocity and pressure in the governing equations, are divided into two components, namely low modes and high modes. Then, the couplings between low modes and high modes, which are not accounted by standard Galerkin method, are considered by AIMs, to improve the accuracy of the numerical results. Further, the multilevel finite element method with hierarchical basis functions is introduced to approach low modes and high modes in an efficient way. As an example, the flow around airfoil NACA0012 at different angles of attack has been simulated by the method presented, and the comparisons show that there is a good agreement between the present method and experimental results. In particular, the proposed method takes less computing time than the traditional method. As a conclusion, the present method is efficient in numerical analysis of fluid dynamics, especially in computing time.