This thesis is dedicated to the study of the turbulent fluid flow in rough-walled rock fractures. Fracture models were generated from 3D scans of fractured rock samples, while fluid flow was simulated numerically by means of FVM-based open-source CFD toolbox OpenFOAM, employing the high-performance computing cluster for the more demanding 3D models. First part of the thesis addresses the issue of fracture geometry. Realistic 2D and 3D fracture models were constructed from 3D scans of upper and lower halves of a fractured rock sample, taking both shear displacement and contact spots into account. Furthermore, we discuss the shortcomings of the available fracture aperture metrics and propose a new aperture metric based on the Hausdorff distance; imaging performance of the new metric is shown to be superior to the conventional vertical aperture, especially for rough fracture surfaces with abundant ridges and troughs. In the second part of the thesis we focus on the fluid flow through the rock fracture for both 2D and 3D cases. While previous studies were largely limited to the fully viscous Darcy or inertial Forchheimer laminar flow regimes, we chose to investigate across the widest possible range of Reynolds numbers from 0.1 to 10^6, covering both laminar and turbulent regimes, which called for a thorough investigation of suitable turbulence modeling techniques. Due to narrow mean aperture and high aspect ratio of the typical fracture geometry, meshing posed a particularly challenging problem. Taking into account limited computational resources and a sheer number of model geometries, we developed a highly-optimised workflow, employing the steady-state RANS simulation approach to obtain time-averaged flow fields. Our findings show that while flow fields remain mostly stationary and undisturbed for simpler contactless geometries, emergence of contact spots immediately triggers a transition to non-stationary flow starting from Re ∼ 10^2, which is reflected by the streamline tortuosity data. This transition disrupts the flow pattern across the fracture plane, causing strong channeling and large separation bubbles, with area of the latter being much larger than the generating contact spots. Adverse influence of the contact spots on the overall permeability is strong enough to override any benefits of aperture increase during shear and dilation. Contactless 3D models can to a certain degree be approximated by their 2D counterparts. Lastly, we investigate the influence of both shearing and contact spots on the overall permeability and friction factor of the fracture, drawing a parallel to the well-studied area of turbulent flow in rough-walled pipes and ducts. Unlike the latter, 3D curvilinear fracture geometries exhibit a gapless laminar–turbulent transition, behaving as a hydraulically rough channel in the turbulent range as the shear displacement increases.

The shear behaviour of rough rock joints was investigated by both laboratory testing and numerical simulation. The most powerful servo-controlled direct shear box apparatus in the world with normal forces up to 1000 kN, shear loading up to 800 kN and frequencies up to 40 Hz under full load was used to investigate the shear strength of schistose rock blocks with dimensions of up to 350 x 200 x 160 mm in length, width and height, respectively. The experiments were performed to study the behaviour of rough rock joints under constant normal load, constant normal stiffness and dynamic boundary conditions. The joint surface of rock specimen was scanned 3-dimensional at the initial stage before shearing by new 3D optical-scanning equipment. The 3D-scanner data were used to estimate the joint roughness coefficient (JRC) and to reconstruct rough surface of rock discontinuities in numerical models. Three dimensional numerical models were developed using FLAC3D to study the macro and micromechanical shear behaviour of the joints. Numerical simulation results were compared to experimental results. Three dimensional characteristic of the joint surface including micro-slope angle, aperture, contact area and normal stress distribution were determined and analyzed.

Increasing of the oil recovery is an essential task of the reservoir engineers. Many highly fractured carbonate reservoirs with low matrix porosity and low matrix permeability in the world have been depleted with low ultimate recovery due to the improper management. To obtain higher oil recovery and an optimum oil production, a three dimensional geological model for the reservoir is needed. The reservoir model should be calibrated via the history matching process which makes the model reliable for forecasting and development planning. An investigation has been carried out to find the most sensitive parameter which affects the matching between the real production data and the simulated production data. After the calibration of the model, some prediction scenarios have been run to realize the future performance of the reservoir. Three wells have been suggested and included in the simulation as producers in some prediction cases and as gas injectors in another other case. The results show that the three producers with low production rate will assist in increasing the recovery and by converting those to gas injectors will assist more especially when the oil production rate is higher. An oil field in Kurdistan region has been selected for this research.

The shear behaviour of rough rock joints was investigated by both laboratory testing and numerical simulation. The most powerful servo-controlled direct shear box apparatus in the world with normal forces up to 1000 kN, shear loading up to 800 kN and frequencies up to 40 Hz under full load was used to investigate the shear strength of schistose rock blocks with dimensions of up to 350 x 200 x 160 mm in length, width and height, respectively. The experiments were performed to study the behaviour of rough rock joints under constant normal load, constant normal stiffness and dynamic boundary conditions. The joint surface of rock specimen was scanned 3-dimensional at the initial stage before shearing by new 3D optical-scanning equipment. The 3D-scanner data were used to estimate the joint roughness coefficient (JRC) and to reconstruct rough surface of rock discontinuities in numerical models. Three dimensional numerical models were developed using FLAC3D to study the macro and micromechanical shear behaviour of the joints. Numerical simulation results were compared to experimental results. Three dimensional characteristic of the joint surface including micro-slope angle, aperture, contact area and normal stress distribution were determined and analyzed.

The anisotropic behavior of a rock mass with persistent and planar joint sets is mainly governed by the geometrical and mechanical characteristics of the joints. The aim of the study is to develop a continuum-based approach for simulation of multi jointed geomaterials. There are two available numerical techniques for the strain-stress analysis of rock masses: continuum-based methods and discontinuum based methods. Joints are simulated explicitly in discontinuous methodology. This technique provides a more accurate description for the behavior of a rock mass. However, in some projects, the explicit definition becomes impractical, especially with increasing number of joints. Besides, the calculation efficiency will be significant reduced as the number of joints increases within the model. Considering the above mentioned shortcomings of the discontinuous method, the continuum-based approach is widely used in rock mechanics. Within the continuum methods, the discontinuities are regarded as smeared cracks in an implicit manner and all the joint parameters are incorporated into the equivalent constitutive equations. A new equivalent continuum model, called multi-joint model, is developed for jointed rock masses which may contain up to three arbitrary persistent joint sets. The Mohr-Coulomb yield criterion is used to check failure of the intact rock and the joints. The proposed model has solved the issue of multiple plasticity surfaces involved in this approach combined with multiple failure mechanisms. The multi-joint model is implemented into FLAC and is verified against the distinct element method (UDEC), analytical solutions, and experimental data. Uniaxial compression tests with artificial rock-like material (gypsum) are carried out in the laboratory in order to verify the developed constitutive model and to investigate the behavior of jointed specimen. Samples with two crossing joints covering more than 20 angle configurations and two different property sets were prepared and tested. Simulation results are in good agreement with experimental observations. The developed model is applied to two potential practical applications: the stability analysis of a slope and a tunnel under different stress conditions. Finally, the main achievements of the whole PhD study are summarized and future research work is proposed.

Mining operations produce a multiscale network of fractures in the coal seams. Permeability evolution in rocks is important for coal bed methane (CBM) and shale gas exploitation as well as for greenhouse gas storage. Therefore, this work presents laboratory tests and a coupled model using PFC3D and FLAC3D to simulate the stress induced permeability evolution in coal samples. Basic mechanical properties are determined via lab testing. The spatial distributions of different components inside the reconstructed samples produce a significant heterogeneity based on CT technique. A newly developed experimental system is employed to perform 3-dimensional loading and to measure the flow rate simultaneously. The evolution process is described by 5 distinct phases in terms of permeability and deformation. Triaxial tests are simulated with PFC3D using a novel flexible wall boundary method. Gas seepage simulations are performed with FLAC3D. Relations between hydraulic properties and fracture data are established. Permeability and volumetric strain show good nonlinear exponential relation after a newly introduced expansion point. Piecewise relations fit the whole process, the expansion point can be treated as critical point. The structural characteristics of the samples influence this relation before and after the expansion point significantly.