A numerical scheme is presented for simulating electrokinetic microfluidics in systems with arbitrary morphology. This scheme is based on a numerical solution of the coupled Poisson, Nernst-Planck, and Navier-Stokes equations. While traditional finite-difference methods were used to resolve the first two problems, the lattice Boltzmann method was applied to the latter. The developed numerical approach was used for the simulation of electroosmotic flow through a simple cubic array of hard ( impermeable, nonconducting) micro-sized spheres. Volumetric electroosmotic flow was studied for dependence on electrical field strength, zeta-potential at the solid-liquid interface, electrical double layer interaction, and numerical grid resolution. Colloid stability and electrokinetics in microfluidic devices with particulate or monolithic fixed-bed elements represent two potential applications of this work.