A numerical model of the ion beam in field emission electric propulsion thrusters is presented. The objective of the study is to develop a computational tool capable of predicting the effects of thruster configuration parameters (geometry, voltage distribution between the electrodes) on thruster performance (beam divergence, ion exhaust velocity). A two-dimensional, steady model has been adopted. A simple emission model has been assumed, neglecting the detailed, microscopic description of the physical emission mechanism. The electric field is evaluated by a finite element scheme; boundary conditions are provided by an auxiliary technique based on the placement of fictitious charges in proximity to the electrode surfaces. The beam is simulated through a discrete number of charged particles; the effects of the associated space charge density distribution on the electric field is considered. The results of the simulations of typical thruster configurations are presented and discussed. The need to take space charge effects into account to attain a sufficiently accurate solution is indicated: the beam charge density modifies the electric potential distribution generated by the electrodes voltage difference, causing a decrease in the local electric field at emitter tip and an increase in beam divergence. Neutralization at a finite distance must be taken into account to fully appreciate the effects of different voltage distributions between the electrodes: a higher accelerator electrode voltage determines a lower ion exhaust velocity and a larger beam divergence. The capability of an additional, decelerating neutral electrode to focus the ion beam is also investigated, showing that only little improvement in beam containment is attained.