In this paper we numerically investigate the motion of viscoelastic liquids passing through two-dimensional periodic arrays of cylindrical particles using the finite element method. The viscoelastic liquid is modeled by the Chilcott-Rallison version of the finitely extensible, nonlinear elastic (FENE) dumbbell model. The permeability and the viscoelastic stress distribution are studied as functions of the dimensionless relaxation time De and the dimensionless wave number kD, where k= 2φ-/λ is the wave number, λ is the distance between the particles in the flow direction and D is the cylinder diameter. The porosity and D are held fixed. Our simulations show that for a fixed value of De the viscoelastic permeability increases with kD, but, as is the case for Newtonian fluid [Alcocer, Kumar, and Singh, Phys. Rev. E 59, 711 (1999)], this increase is not monotonic. The permeability decreases between kD≈5.0, where it is locally maximum, and kD≈7.7, where it is locally minimum. The difference between the locally maximum and minimum values of permeability increases with increasing De. When De= 0(1) the locally minimum value of the permeability is ∼40% smaller than the value at the local maximum. This implies that a substantial change in permeability can be achieved by changing the distance between the particles in the flow direction while keeping De, D, and porosity fixed. [ABSTRACT FROM AUTHOR]