The Galerkin/least-square hp finite element method is employed to solve the problem of steady flows of a sphere falling in a tube filled with viscoelastic fluids, the diameter ratio being 2.0. A global a posteriori error estimation method is proposed and applied to demonstrate the convergence and error distributions of the solutions. By adopting appropriate high-order interpolation distributions the computations have discovered three boundary layers of the extra stress components near the rear stagnation point when the solutions approach to the limiting points for the upper convected Maxwell (UCM) and Oldroyd-B fluids, meanwhile being able to satisfy the positive definiteness of the Finger stress tensor. These boundary layers, characterized by rapidly increasing stress gradients and rapidly decreasing thickness, pose a great challenge for accurate numerical predictions and compose the major difficulty for reaching considerably higher Deborah numbers. The limiting Deborah number for the Oldroyd-B fluid with the viscosity ratio of 0.5 is about 1.3, which agrees with Chauviere and Owens' prediction and is considerably lower than the limiting Deborah number of 2.2 for the UCM fluid. The computations for a finitely extendable nonlinear elastic spring (FENE) fluid with low extensibility of the molecular chain can easily exceed the limiting Deborah number of the UCM fluid, but with the high extensibility the solutions exhibit the same characteristics of the stress boundary layers as in the UCM case. This indicates that the near-singular behaviors of the UCM or Oldroyd-B solutions near the limiting points are rooted in the infinitely extensible Gaussian chain in these two models. (C) 2003 Elsevier Science B.V. All rights reserved.