A higher order Galerkin finite element scheme for simulation of two‐dimensional viscoelastic fluid flows has been developed. The numerical scheme used in this study gives rise to a stable discretization of the continuum problem as well as providing an exponential convergence rate toward the exact solution. Hence, with this method, spurious oscillatory modes are effectively eliminated by increasing the order of the interpolant within each subdomain. In our calculations, an upper limit for the Weissenberg number due to numerical instability was not encountered. However, the memory requirements of the discretization grow quadratically with the polynomial order. Consequently, the maximum attainableWeis determined by the availability of computational resources. The algorithm was tested for flow of upper convected Maxwell and Oldroyd‐B fluids in the undulating tube problem and it was subsequently applied to viscoelastic flow past square cylindrical arrangements. The results obtained show no increase in the flow resistance with increasing elasticity. These findings are in contrast to the experimental data reported in the literature. However, recent experimental investigations have indicated that this dramatic increase in flow resistance is due to a purely elastic instability. Hence, numerical simulation based on two‐dimensional steady flow of viscoelastic fluids cannot be expected to capture this phenomenon.