Flow and transport in porous media are essential in many processes in mechanical, chemical, and petrochemical industries. Despite the wide variety of applications and intensive research efforts, the complex hydrodynamics of these systems is still not fully understood, which renders their design and scale-up difficult. Most porous media have a particulate origin but some are composed of long particles/fibres and, therefore, are considered as fibrous media. They are encountered in a variety of modern technological applications, predominantly in the manufacturing of fibre-reinforced composites, with extensive use in the aerospace and automobile industries. The aim of this thesis is to further develop our understanding of the drag closures, i.e. the connection between microstructure (particle shape, orientation and arrangement) and macroscopic permeability/drag. To address this problem, we employ fully resolved finite element (FE) simulations of flow in static, regular and random arrays of cylinders (and other shapes) at low and moderate Reynolds numbers. Asymptotic analytical solutions at both dense and dilute limits are used to construct drag relations that are universal, i.e. valid for all porosities. Those relations are needed for coupling of the fluid and solid phases (particles) in multi-phase flow codes. The numerical experiments suggest a unique, scaling power law relationship between the permeability and the mean value of the shortest Delaunay triangulation edges, constructed using the centers of the fibres (which is identical to the averaged second nearest neighbor fibre distances). It is complemented by a closure relation that relates the effective microscopic channel lengths to the macroscopic porosity. This percolating network of narrow channels controls the macro flow properties. From our fully resolved FE results, for both ordered and random fibre arrays, we find that (i) the weak inertia correction to the linear Darcy relation is third power in superficial velocity, U, up to small Reynolds number, Re~1-5. When attempting to fit our data with a particularly simple relation, (ii) a non-integer power law in U performs astonishingly well up to the moderate Re~30. However, for randomly distributed arrays, (iii) a quadratic correction performs quite well as used in the Forchheimer (or Ergun) equation, from small to moderate Re. Finally, the universal fluid-particle drag relations have been incorporated into a coarse FE two-phase framework, based on coupling an unstructured FE mesh and a soft-sphere discrete element method (DEM) for moving particles. The mesh is a dynamic Delaunay triangulation based on the particle positions. This provides a framework for FE method discretization of the equations of fluid dynamics as well as a simple tool for detecting contacts between moving particles.