Aiming at simulating the viscoelastic flow problems involving three-dimensional, time-dependent complex flow, a fully three-dimensional (3D) numerical method using an implicit finite volume formulation is developed in the cylindrical coordinate system. We focus our attention on the well known swirling flows in a confined cylinder with a rotating bottom lid, and through careful comparisons with the available experimental results, the suitability of different constitutive models used to describe viscoelastic fluids is evaluated. Three difficulties arising from the coordinate transformation in a structured mesh system are removed by some simple, but effective measures. With this method, the velocity fields measured in experiments are accurately predicted, and the development of vortex breakdown in terms of Reynolds number (Re) and the height to radius ratio ( H/ R), as well as the time, observed in experimental visualisations are exactly reproduced for Newtonian flows. The influence of elasticity on the vertex breakdown processes is investigated for the low viscosity dilute polyacrylamide (PAA) solutions described by the upper converted Maxwell (UCM) model. It is confirmed numerically that the existence domain of vortex breakdown occurs at significantly greater Re than that for Newtonian solvent of similar shear viscosity. It seems that the centrifugal force leading to vortex breakdown in Newtonian fluids is balanced by the normal stresses in elastic fluids. As a result, vortex breakdown would be significantly delayed in terms of Re or even suppressed with increasing elasticities. Numerical simulations of the steady swirling flow involving viscoelastic fluids with a high constant viscosity (low Re) shows that the flow field could partially or entirely reversed depending on the level of fluid elasticity. Our numerical results demonstrate that at a low value of Re, for the fluids with a medium high elasticity characterized by the elastic number El = O(1), an elastic ‘ring’ vortex is developing and growing with increasing elasticities at the edge of the rotating lid, and trying to push the Newtonian vortex to the radial centre of the rotating lid, and finally it occupies the whole flow domain at a certain level of elasticity. It is also found that a small weak edge ring vortex appears at the edge of the rigid cover and counter-rotates with the main vortex when the aspect ratio H/ R ≥ 1, and its size tends to increase with increasing H/ R. The double-cell flow structure for strong shear-thinning viscoelastic fluids observed in visualisation has been successfully predicted numerically by using the simplified Phen-Thien–Tanner (SPTT) model. With steady-state calculations, at a given inertia level, when the flow field is dominated by elasticity, we found that there is a small counter-rotating central ring vortex appears around the centre of the rotating lid, which was usually reported as a sign of instability of the flow field in experimental work. In the second part of this work, we will focus on the numerical inspection of the highly unsteady spiral vortex flows of viscoelastic fluids observed in experiments by using time-dependent calculations, thus, checking the stability criterion of McKinley et al. [G.H. Mckinley, P. Pakdeland, A. Öztekin, J. Non-Newtonian Fluid Mech. 67 (1996) 19–47], and evaluating the suitability of different constitutive equations.