VGeST - Virtual Geoscience Simulation Tools

Fracture modelling

This research was undertaken in 2012 by Qinghua Lei, Liwei Guo and Dr.Jiansheng Xiang under the supervision of Dr. J-P Latham

2D Fracture models

(1) The Fracture Model 

A combined single and smeared crack model is applied in the context of the combined finite-discrete element method (FEMDEM) to simulate fracturing for quasi-brittle materials both in tensile and compressive stress fields. It is based on the accurate representation of the stress-strain curves obtained in physical experiments. A Mohr-Coulomb with tensile strength cut-off failure criterion has been implemented and a material heterogeneity option is available. 

(2) Geomechanical modelling of analogue fracture networks 

Advanced numerical technologies for the behaviour of fractured rock masses have application to reservoir flow simulation, geothermal power and geological disposal of nuclear wastes. In this study, ~10 m scale analogue fracture patterns, i.e. Kilve 6×6 m and Bristol 18×8 m networks, are extracted from the 2D outcrop maps of natural fracture systems that involve joint sets with complicated intersections, bends and segmentations. In situ stresses of the strike-slip faulting systems are calculated according to the data of reservoir overburden, pore fluid pressure and Poison’s ratio, etc. The response of the fractured rock models to the application of biaxial effective stresses is modelled using the in-house code Y2D based on the combined finite-discrete element method (FEMDEM). The two fractured reservoir models exhibit significant fracture-dependent stress heterogeneity, as the contours of principal stresses shown below. 

Kilve 6x6m Oblique Sets

Bristol Kilve Bed 3 18x8m ~ Orthogonal Sets


Maximum Principal Stress (Pa)

Minimum Principal Stress (Pa)

Maximum Principal Stress Kilve Bed 3 (Pa)

An important advantage of this modelling methodology is its capability of realistically describing the distributions of normal aperture and shear displacement between fracture walls. Further post-processing using the rose diagram technique provides more details about the magnitudes of aperture and shear displacement in different directions. 

Normal apertures under boundary stresses - Kilve

Shear displacement under boundary stresses - Kilve

Rose plot normal aperture

Rose plot shear displacement

Fracture apertures respond to in-situ stress states. Flow simulation based on mapped analogues or geomechanically modelled fracture networks that accommodate aperture variability should more realistically reflect the permeability properties of natural reservoirs and permit more reliable prediction of potential oil recovery than conventional approaches that assume a constant aperture value for all fractures. 

(3) Modelling fracture network growth 

Natural fractured reservoirs around the world are a very significant source of hydrocarbons. The difficulty in characterising such reservoirs is mainly attributed to the lack of sufficient sub-surface data to create realistic fracture models. Such fracture models can be generated stochastically, but these have limited abilities to capture realistic cross-cutting, branching and truncation relations. To overcome these limitations, the fracture model in the FEMDEM method has been used to create suitable fracture networks from virgin rock. 

Multilayer of limestones and shales at Kilve

Fractures in multilayer systems of limestones and shales have intrigued structural geologists for decades.  There appears to be a spacing to thickness relationship that may be governed by an underlying mechanical principle. Field measurements may be able to tell us something about the behaviour of natural rock systems. Many studies investigating mechanism of failure using analytical models, model experiments or numerical simulation  consider the response to layer-parallel extension. 

One application of VGeST  was to use Y2D in plane strain. First, a shale-limestone vertical section model was investigated by Latham, Guo, Wang and Xiang  (2011). The limestone layer is sandwiched by thicker shale layers. This three-layer model is confined by pressure in the horizontal direction and compressed in the vertical direction. Material heterogeneity is introduced into the model according to a Weibull distribution.

Three layer sandwich - homogeneous

Three layer sandwich - heterogeneous

 With regard to the homogeneous material, it can be seen that the shear fractures in shale layers are formed at 45 degrees, and in the limestone layer the vertical tensile fractures are evenly distributed. In this case, the fracture spacing to layer thickness ratio in the limestone layer is approximately 0.8-1.2. In contrast, for the heterogeneous material although the shear fractures in the shale layers are similar to the ones in the homogeneous case, the vertical tensile fractures in the limestone layer are distributed more randomly. Note that pure shear fracturing is captured in the central area of the limestone layer. Some other factors such as friction coefficients on bedding planes and confining pressures, have also been investigated. The relation between fracture saturation and confining pressure for the homogeneous material is shown in the plot below.  

saturation of fractures and the effect of confinement

A seven-layer all limestone multilayer with some variability of properties has been modelled. This model is stretched in the horizontal direction. Different fracture spacing to layer thickness ratios have been obtained in alternative layers. Furthermore, the normal apertures have been extracted in the post-processing as shown  by different colours in the figure. Note that in some areas the fractures are linked across bedding planes which is a common occurence in some natural systems and has significance for potential pathways for fluid flow.

Seven layer limestone sequence - normal aperture (m)

3D Fracture Modelling

The original two-dimensional fracture model in the FEMDEM method has been developed into three dimensions. A few validation tests are presented here. The first example is a three-point bending test. The size of this beam is 580 mm × 70 mm × 30 mm. The two supporting wedges move up at a constant velocity and the top wedge moves down at the same velocity. The material used in this test is typical concrete. It can be seen that the fracture develops in the middle of the beam, and the concentration of horizontal tensile stresses around the fracture tip is precisely captured.

3 Point bend test in concrete beam

horizontal tensile stress distribution in the three-point bending test

 The second example is a Brazilian test. A cylinder disc of typical rock properties is compressed by two steel platens. The diameter of the disc is 40 mm and the thickness is 17 mm. The failure mechanism includes tensile splitting from the centre of the disc along the loading axis and crushing near the loading platens.

Brazilian test of a rock specimen

One example of collision phenomena has also been modelled using the 3D fracture model in the FEMDEM method. In the examples, two identical concrete spheres move towards each other at an initial horizontal velocity of 10 m/s. The diameter of the spheres is 24 mm. During the collision process, most of the energy is released by the fracturing of spheres into small pieces.

Test of two spheres in direct collision resulting in fracture and fragmentation