Non-Brownian fibers commonly flocculate in flowing suspensions at relatively low concentrations ( < 1% by weight). We have developed a particle-level simulation technique modeling fibers as chains of rods connected by hinges to probe fiber flocculation. The model incorporates fiber flexibility, irregular fiber equilibrium shapes, and frictional fiber interactions. Model fibers reproduce known orbits of isolated rigid and flexible fibers in shear flow. Simulation predictions of first normal stress differences in homogeneously dispersed, dilute flexible fiber suspensions agree with experimental data. Fiber features such as flexibility and irregular equilibrium shapes strongly impact single fiber and suspension behavior. Fibers aggregate in simulations with interfiber friction, in the absence of attractive forces between fibers. Strong flocculation is observed in suspensions of stiff fibers with irregular equilibrium shapes. Flocs contain many fibers with three or more contact points, and derive cohesiveness from elastic energy held in fibers-consistent with the elastic interlocking mechanism of flocculation. At higher concentrations (nL(3) approximate to 100, where n is the fiber number density and L is the fiber length), coherent fiber networks form in simulations. Average numbers of contacts per fiber and contact force magnitudes in sheared and static networks are compared with existing fiber network theory predictions.