The characteristics of fibers suspensions depend on the properties of fibers, the suspending fluid, and fiber-fiber interactions. This thesis demonstrates the development and application of a novel coupled method (lattice Boltzmann and finite element methods) to investigate these relationships. Fibers are modeled as flexible rod particles which are simulated by the finite element method. The fluid flow that causes the fibers to deform is calculated by the lattice Boltzmann method. The method is extended from the two dimensional case to the three dimensional case. Results from simulation show the rigid fiber in simple shear flow produces a good agreement for orientation of a fiber relative to the theoretical study by Jeffery (1922). The flexible fiber exhibits an increase on the rotational period from the rigid fiber due to more deformation shape is revealed during rotation. The simulation technique demonstrates the ability to simulate fiber-fiber interactions to further study of relative viscosity of suspensions in shear flow. Simulation results show that fiber orientation and relative viscosity depend on the fiber characteristics (fiber aspect ratio, fiber flexibility, and volume fraction). The results are verified against known experimental measurements and theoretical results. The broad aim of this research is to better understand the behavior of fibers in fluid flow. It is hoped that future researchers may benefit from the new technique and algorithms developed here.