The static structure and conformational properties of salt-free polyelectrolyte solutions are studied using many-chain Monte Carlo simulations and integral equations. The polymer molecules are modeled as freely jointed tangent sphere chains with a hard sphere plus screened Coulomb potential between the beads. Monte Carlo simulations are for the most part in agreement with previous molecular dynamics simulations of a related model of polymers but the present work facilitates an unambiguous comparison between theory and simulation. An integral equation theory is implemented where the properties of the polymer are obtained from a single chain simulation where the beads interact via the bare potential plus a self-consistently determined solvation potential. The theory overestimates the degree of liquidlike order and underestimates the size of the chains (when compared to many chain simulations) especially for long chains and high densities. A more approximate theory that employs a thread model for the polymers is actually more accurate.