A model is proposed for simulating the motion of flexible fibres in fluid flow. Care has been taken to include typical papermaking conditions into the validity range of the model. Fibres are modelled as chains of fibre segments, whose motion is governed by Newton's second law. The fluid motion is calculated from the three-dimensional incompressible Navier-Stokes equations. By enforcing momentum conservation, the two-way coupling between the solids and fluid phase is taken into account. Fibre-fibre interactions as well as self-interactions include normal, frictional and lubrication forces. Furthermore, the model considers nonlaminar fibre-fluid interactions and particle inertia. Simulation results were compared with experimental data found in the literature. The model predicts very well the orbit period of rigid fibre motion in shear flow. Quantitative predictions were made for the amount of bending of flexible fibres in shear flow. It was also possible to reproduce the different regimes of motion of flexible fibres in shear flow, ranging from rigid motion to coiled motion and self-entanglement.