This work presents a finite volume method used to solve the relevant coupled equations for general electro-osmotic flows of viscoelastic fluids, using the upper-convected Maxwell and the simplified Phan-Thien-Tanner models. Three different implementations of the electro-osmosis physical models were carried out, which depend on the required level of approximation. In the first implementation, the Poisson-Nernst-Planck equations were incorporated into the code and the electric charge distribution required to quantify the electric field forcing of the momentum equation is calculated from these fundamental equations. The second implementation is an approximation in which a stable Boltzmann distribution of ions is assumed to occur in the electric double layer, represented by the Poisson-Boltzmann equations. Finally, the so-called Debye-Hikkel approximation was also implemented in the Poisson-Boltzmann-Debye-Huckel model, which is valid for cases with a stable Boltzmann distribution of ions characterized by a small ratio of electrical to thermal energies. Numerical simulations were undertaken in a cross-slot geometry, to investigate the possible appearance of purely electro-elastic flow instabilities, by considering the effect of the electric field. For pure electro-osmotic flow, i.e., in the absence of an imposed pressure gradient, we were able to capture the onset of an unsteady asymmetric flow above a critical Weissenberg number, which is lower than the corresponding value for pressure-driven viscoelastic flow. Additionally, we derive analytically the fully-developed electro-osmosis driven channel flow of polymer solutions described by the FENE-P and PTT models with a Newtonian solvent and assuming a thin electric double layer. (c) 2012 Elsevier B.V. All rights reserved.