  We present an approach for numerically simulating the dynamics of flexible
fibers in a three-dimensional shear flow using a scalable immersed boundary
(IB) algorithm based on Guermond and Minev's pseudo-compressible fluid solver.
The fibers are treated as one-dimensional Kirchhoff rods that resist
stretching, bending, and twisting, within the generalized IB framework. We
perform a careful numerical comparison against experiments on single fibers
performed by S. G. Mason and co-workers, who categorized the fiber dynamics
into several distinct orbit classes. We show that the orbit class may be
determined using a single dimensionless parameter for low Reynolds flows.
Lastly, we simulate dilute suspensions containing up to hundreds of fibers
using a distributed- memory computer cluster. These simulations serve as a
stepping stone for studying more complex suspension dynamics including
non-dilute suspensions and aggregation of fibers (also known as flocculation).
