A numerical simulation method has been developed to predict the motion of fiber suspensions in various flows by using the particle simulation method (PSM), in which a fiber is modeled by an array of spheres. The hydrodynamic interaction among fibers is considered by decomposing into intra-fiber and inter-fiber interactions. For the intra-fiber case, 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. For the inter-fiber case, only the near-field lubrication force is considered beta een spheres belonging to different fibers. The methodology was applied to predicting the microstructure and rheological properties of rigid fiber suspensions in simple shear now and fiber orientation in two-dimensional diverging now by combining this method with the finite element method (FEM). In the former, the overshoot of suspension viscosity, due to the transient change of the microstructure, was observed in semi-dilute to concentrated systems. In the latter, the calculated fiber orientations agreed with expectations of theory.