Abstract A numerical method is constructed for two-dimensional Navier–Stokes flows in a circular domain. Adaptive computational grids are employed to resolve the boundary layer and the boundary vorticity values are obtained in high order of accuracy. Numerical results show a quite good agreement with the asymptotic expansion calculation (Kim, SIAM J. Appl. Math. 58(5) (1998) 1394–1413). Extension to general shape domains is briefly addressed.