We describe a computational method for the numerical simulation of three-dimensional transient flows of polymer solutions that extends the work of Harlen et al. [O.G. Harlen, J.M. Rallison, P. Szabo, A split Lagrangian-Eulerian method for simulating transient viscoelastic flows, J. Non-Newtonian Fluid Mech. 60 (1995) 81-104]. The method uses a Lagrangian computation of the stress together with an Eulerian computation of the velocity field. Adaptive mesh reconnection based on Delaunay tetrahedra is used to ensure well-shaped elements. Additional shape-quality improvement procedures are developed to improve the algorithm. We validate the method for the benchmark problem of a rigid sphere falling in a cylindrical pipe. Inertia is neglected. We compare results for the axisymmetric case with previous work (using a FENE model), and then consider the off-axis non-axisymmetric case. In the latter case, we find that as the sphere falls, it drifts across the pipe, a phenomenon previously observed in experiments but not fully explained. The physical mechanisms that cause the time-dependent drift are identified, and a simple model based on the normal stresses in the fluid is shown to predict the magnitude of the drift velocity. We also consider a second benchmark problem involving a constriction in an axisymmetric pipe. Numerical difficulties associated with ill-shaped elements near the concave boundary arise for higher Weissenberg numbers. The merits and drawbacks of the new numerical method, and its applicability to various flow geometries are discussed. (C) 2010 Elsevier B.V. All rights reserved.