Within the framework of the high-order finite volume (FV) method, a high-order gas kinetic flux solver (GKFS) is developed in this work for simulation of two-dimensional incompressible flows. Generally, in the conventional high-order FV method, the inviscid and viscous fluxes are treated separately. However, different from the conventional high-order FV method, the high-order GKFS evaluates the inviscid and viscous fluxes simultaneously from the local asymptotic solution to the Boltzmann equation, which consists of the equilibrium distribution function and its substantial derivative at the cell interface. By introducing a difference scheme with the high-order accuracy in space to discretize the substantial derivative, a high-order accurate local asymptotic solution to the Boltzmann equation can be obtained. The numerical flux of the Navier–Stokes equations can then be calculated by the moments of the local asymptotic solution. Since this local asymptotic solution is relatively simple, the numerical fluxes of the Navier–Stokes equations can be given explicitly for the high-order GKFS, which is the function of the left and the right states and their first-order derivatives. Numerical results showed that the developed solver can achieve the desired accuracy on both the quadrilateral mesh and the triangular mesh and its efficiency is higher than the second-order counterpart when achieving comparable accuracy of solution.