A numerical simulation method has been developed to predict the rheological behavior of fiber suspensions in shear flow by using the particle simulation method (PSM), in which a fiber is modeled by arrays of spheres. The hydrodynamic interaction among fibers is considered by decomposing into intra- and inter- fiber ones. In the former, the many-body problem is solved by calculating the mobility matrix for each fiber to obtain the hydrodynamic force and torque exerted on each sphere. In the latter, only the near-field lubrication force is considered between spheres belonging to different fibers. The methodology was applied to predict the microstructure and the rheological properties of suspensions of rigid fibers with aspect rations 5 and 10 at several cancentrations. The overshoot of suspension viscosity due to the transient change from random to planar orientation in the microstructure was observed in semi-dilute to concentrated regimes.