A numerical simulation of multiple flexible fibers in suspension in Newtonian simple shear flow is presented. The method used is similar to those of previous recent simulation works by Fan et al. [J. Non-Newtonian Fluid Mech. 74 (1998) 113] and Yamane et al, [J. Non-Newtonian Fluid Mech. 54 (1994) 405], however, the method has been modified to allow a small amount of bending and torsion in the fibers. A restoring moment acts to straighten the fibers as they interact in the flow. It is demonstrated that this simulation can be used to extract basic theological information about the suspension including fiber orientations and suspension viscosity. The viscosity of semi-concentrated to concentrated flexible fiber suspensions are shown to increase by a magnitude of the order 7-10% greater than the equivalent rigid fiber suspension tested. This is in qualitative agreement with previous experimental work by Goto et al. [Rheologica Acta 25 (1986) 119] and Blakeney [J. Colloid Interface Sci. 22 (1966) 324]. The implication is that any constitutive relation involving particulate suspensions described by orientation vectors may quantitatively underestimate suspension viscosity, particularly for fibers of large aspect ratio, or low Young's modulus, whereby the tendency to flex is greater [Rheologica Acta 25 (1986) 1 19]. If particulate deformation were accounted for (by whatever means) in the existing constitutive relationship, predictions of bulk suspension parameters such as viscosity should be noticeably improved. A method is developed to modify an existing rigid-fiber viscosity to an equivalent flexible fiber viscosity, hence improving viscosity prediction ability.