A method is described for obtaining viscoelastic flow solutions, based upon time-integral constitutive equations, using a general purpose CFD package. The method is general enough to be applied to any available software that has rudimentary input and output facilities and can solve a Stoke＇s flow problem. From this basis, flexibility of choice of constitutive equation and computational techniques is available. The method is presented in a form appropriate for solving both planar and axisymmetric flows. Delaunay triangulation is used to reconstruct a mesh for external code, and stress computation procedures are performed on this mesh. The method has only two particular requirements for the CFD package used - it must be able to output nodal values (of position and velocity) to file and it must be able to read body-forces from a file. Two methods of velocity adjustment were compared: an incremental method and a method whereby the viscoelastic stresses were incorporated directly as body-forces. Results and convergence from the two methods were found to be essentially identical, hence the direct body-force method (which is considerably easier to implement) is described in detail. The method is applied to a well-known flow problem of LDPE melt through a 4 : 1 abrupt contraction axisymmetric die. Convergence was obtained up to nearly the highest value of apparent shear rate for which published simulation results are available. Quantitative results for vortex strength, vortex opening angle and Couette correction are presented which are compared with earlier work on the problem using other methods. Agreement is generally good, giving confidence in the method. Simulations of planar flows of the same melt are performed: a decreasing corner vortex was observed. This phenomenon has been observed experimentally for flows of other substances, but is not expected for flows of LDPE melt. A parametric study of a critical strain hardening parameter is conducted to help explain the cause of the results.