A numerical tool has been developed for particle-level simulations of fiber suspension flows, particularly forming of the fiber network structure of paper sheets in the paper machine. The model considers inert fibers of various equilibrium shapes, and finite stiffness, interacting with each other through normal, frictional, and lubrication forces, and with the surrounding fluid medium through hydrodynamic forces. Fiber–fluid interactions in the non-creeping flow regime are taken into account, and the two-way coupling between the solids and the fluid phases is included by enforcing momentum conservation between phases. The incompressible three-dimensional Navier-Stokes equations are employed to model the motion of the fluid medium. The validity of the model has been tested by comparing simulation results with experimental data from the literature. It was demonstrated that the model predicts well the motion of isolated fibers in shear flow over a wide range of fiber flexibilities. It was also shown that the model predicts details of the orientation distribution of multiple, straight, rigid fibers in a sheared suspension. Furthermore, model predictions of the shear viscosity and first normal stress difference were in fair agreement with experimental data found in the literature. Since the model is based solely on first principles physics, quantitative predictions could be made without any parameter fitting. Based on these validations, a series of simulations have been performed to investigate the basic mechanisms responsible for the development of the stress tensor components for monodispersed, non-Brownian fibers suspended in a Newtonian fluid in shear flow. The effects of fiber aspect ratio, concentration, and inter-particle friction, as well as the tendency of fiber agglomeration, were examined in the non-concentrated regimes. For the case of well dispersed suspensions, semi-empirical relationships were found between the aforementioned fiber suspension properties, and the steady state apparent shear viscosity, and the first/second normal stress differences. Finally, simulations have been conducted for the development of paper structures in the forming section of the paper machine. The conditions used for the simulations were retrieved from pilot-scale forming trial data in the literature, and from real pulp fiber analyses. Dewatering was simulated by moving two forming fabrics toward each other through a fiber suspension. Effects of the jet-to-wire speed difference on the fiber orientation anisotropy, the mass density distribution, and three-dimensionality of the fiber network, were investigated. Simulation results showed that the model captures well the essential features of the forming effects on these paper structure parameters, and also posed new questions on the conventional wisdom of the forming mechanics.