References to cite
2018: Slip on wavy frictional faults: Is the 3rd dimension a sticking point?
2017: Stress concentrations around voids in three dimensions: The roots of failure
- Introduction
- 1. How to run
- 2. Analytical tests
- 3. Files and data format
- 4. Notation in Code
- 5. Elements used in 2D and 3D
- 6. Code highlights and applications
- 7. Acknowledgements
Gif of a fault array shearing due to a remote stress. Computed using the 3D code
This software is for academic use. Do not use this in a commercial environment. There are plenty of commercial Boundary Element codes in circulation. It does not require compiling, simply run the script 'AddFilePaths.m' then load the file 'Mainframe.m' in MATLAB (either from the 3D or 2D directory) then start the code by pressing the ▶ (run) symbol in the editor tab of MATLAB, alternatively call "run MainFrame.m" in the cmd window. Once the code starts it should run spitting out diagrams, showing progressbars and supplying results.
This code uses the Boundary Element Method (BEM), specifically the Displacement Discontinuity Method (DDM). Only fault surfaces or closed contours of bodies need to be digitised with boundary conditions placed on these elements. This assumes: The material is isotropic, a linear elastic and that infinitesimal deformation applies.
When calculating an unknown such as slip due a prescribed stress the method is dependent on sampling of the surfaces. Users must increase sampling until they are satisfied the result is static or has a changes in precision that are acceptable. This can be seen clearly changing the sampling of surfaces or lines in the regression tests.
This code was created in windows MATLAB versions:
2015a:2016b:2017a
It’s also been tested in Octave GUI windows versions:
4.0.0:4.0.3
although there are sometimes 'bugs' in the Octave Windows version it tends to run fine. It tends to fail in certain places, such as dealing with strings and drawing transparent objects so slight modifications may be required.
It has also been run on both on Linux and Mac versions of MATLAB (2016b) and works fine. For these versions errors may exist when loading the file paths in newer scripts not tested on these platforms i.e \ /.
I have tried to comment lines as much as possible. References to 'Pollard' are :'Pollard, D.D. and Fletcher, R.C., 2005. Fundamentals of structural geology. Cambridge University Press.'
I have attempted to keep copy-write tags on files from other places. If you spot somewhere I have missed this please email ASAP so I can fix this. Places where I was unsure if I could use such files a info.txt has been left with a link to the file download location. Download the files and place these in the file containing 'info.txt'.
Codes for 2D and 3D are laid out in the same manner. To run:
- Download the file from GitHub
- Unzip the folder and move it to the directory you would like to work from
- If you rename the folder make sure the new name still contains: 'CutAndDisplace'
- Open MATLAB
- Run the file inside the folder called ‘AddFilePaths.m’
- Then run using and editing 'MainFrame.m' in the respective 2D or 3D directory
The code has analytical regression tests built in to test changes have not introduced bugs before a commit. Note that all images below have been created directly using the software and are drawn like this when the regression tests run. These are:
Test by running the file runallregressions.m in 2D/AnalyticalTests folder or the individual scripts further inside this directory.
-
Half space formulation from:
Martel SM. Boundary element scripts http://www.soest.hawaii.edu/martel/Martel.BEM_dir/, University Hawaii at Manoa, 2003 tested against the 3D analytical solution of:
Nikkhoo, M. and Walter, T.R., 2015. Triangular dislocation: an analytical, artefact-free solution. Geophysical Journal International, 201(2), pp.1117-1139.
-
Comparing a basic dislocation from the (Crouch and Starfield, 1983) code to a dislocation described in:
Barber, J., 2010. Elasticity (ed., Vol. 172).(G. Gladwell, Ed.) Michigan.
and
Pollard, D.D. and Fletcher, R.C., 2005. Fundamentals of structural geology. Cambridge University Press.
Numerical vs analytical comparison of Cartesian stress components around a right lateral shear dislocation lying along the x-axis -
Comparing the DDM match of slip and stress surrounding a planar crack to the analytical solutions of
Pollard, D.D. and Segall, P., 1987. Theoretical displacements and stresses near fractures in rock: with applications to faults, joints, veins, dikes, and solution surfaces. Fracture mechanics of rock, 277(349), pp.277-349.
Numerical (dots) vs analytical (green line) comparison of crack wall displacements due to a line crack loaded with a tensile stress
Numerical vs analytical comparison of Cartesian stress components around a 2D dislocation lying along the x-axis opening due to a tensile stress. -
Comparing the match of the DDM to the (Kirsch, 1898) solution for a pressurised hole. This tests the elemental locking used to stabilise the DDM result. The (Kirsch, 1898) formulas are taken from:
Pollard, D.D. and Fletcher, R.C., 2005. Fundamentals of structural geology. Cambridge University Press.
Numerical vs analytical comparison of Cartesian stress components around a pressurised hole -
Comparing the frictional formulation for a planar crack in 2D and its match to the analytical solution for the slip profile of a crack with increasing cohesion towards the crack tips. Analytical slip profile given in:
Bürgmann, R., Pollard, D.D. and Martel, S.J., 1994. Slip distributions on faults: effects of stress gradients, inelastic deformation, heterogeneous host-rock stiffness, and fault interaction. Journal of Structural Geology, 16(12), pp.1675-1690.
-
Comparing the DDM solution for gravitational stress under valley described in (Martel and Muller, 2000) to the analytical solution of:
Savage, W.Z., Powers, P.S. and Swolfs, H.S., 1984. In situ geomechanics of crystalline and sedimentary rocks; Part V, RVT, a Fortran program for an exact elastic solution for tectonics and gravity stresses in isolated symmetric ridges and valleys (No. 84-827). US Geological Survey,.
Numerical vs analytical comparison of Cartesian stress components due to valley incision -
Comparing the DDM result of stresses across inhomogeneous elastics to the analytical solution for the stresses across a dual elastic annulus given in:
Crouch, S.L. and Starfield, A.M., 1982. Boundary element methods in solid mechanics: with applications in rock mechanics and geological engineering. Allen & Unwin.
Thanks to Dr Robert Simpson of Glasgow University for sitting down and taking the time walk through the (Crouch and Starfield, 1983) book with me and explain how this is done in other BEM formulations. It was very enlightening and helped me understand how to formulate this in the code.
Comparison of the radial stress tensor components from numerical (dots) and analytical (lines) solutions to the dual material annulus problem.
Test by running the file runallregressions.m in 3D/AnalyticalTests folder or the individual scripts further inside this directory.
-
Comparing the superposition of two of the triangular dislocations of Nikkhoo, M. and Walter, T.R., 2015. to:
the (Okada, 1985) solution coded by François Beauducel.
This is done in the Nikkhoo paper but it’s good to check his code is not changed.
-
Comparing the 3D DDM to the slip profile of a penny shaped crack of (Sneddon, 1967) given concisely in:
Segall, P., 2010. Earthquake and volcano deformation. Princeton University Press.
Comparison of analytical (lines) and numerical (dots) crack wall displacement profiles for a 3D penny shaped crack subject to shear and tensile loads -
Comparing the 3D DDM to the ground surface displacement from the (Mogi, 1958) point source approximation given in:
Segall, P., 2010. Earthquake and volcano deformation. Princeton University Press. This uses the elemental locking much like the (Kirsch, 1898) test in 2D.
Analytical (lines) vs numerical (dots) comparison of ground displacement due to a 'Mogi' point source. -
Comparing the 3D DDM to the 2D analytical solutions for:
- Friction
- Gravitational stresses
- Inhomogeneous materials
Described above for 2D.
Comparison of the numerical and analytical results of crack wall displacement. The elliptical frictionless crack wall displacement profile is shown in blue and the displacements when the crack has a linear frictional profile is shown in green. Black dots show the numerical result. -
Comparing the 3D DDM solution to the (Eshelby, 1959) solution of a void. Using the code from:
Meng, C., Heltsley, W. and Pollard, D.D., 2012. Evaluation of the Eshelby solution for the ellipsoidal inclusion and heterogeneity. Computers & Geosciences, 40, pp.40-48.
Formats that the code uses:
The 2D code uses line data, each straight segment of a line is used as an element. These can be imported as: Shape files(.shp) with multiple lines or alternatively you can define lines your self.
The 3D code uses triangulated surfaces, each triangle then acts as an element (see gif above). These can be imported using these file types: Stl (.stl) : only the non binary format, option available as export from the OS Meshlab software. http://meshlab.sourceforge.net/ gocad ascii (.ts) : a common format in geological software. Easy to read the text files. You could attempt to make 3D surfaces yourself if you have lots of time.
Once these are defined in Step 1 then the user can continue.
The code is run using the files 'Mainframe'. This is a list of steps that call functions as they go, these give an indication of how to use the software. If you put a breakpoint at the beginning of each step you should be able to slowly work your way through the code and work out how it works. The steps are summarised below:
-
Step 1. Import the line/surface, defining the boundary and what it represents (i.e. locked elements, interfaces) for the boundary element model. Importing from file or defining directly.
-
Step 2. Define fullspace/halfspace and elastic constants of the material. This goes on to create the necessary data structure for you.
-
Step 3. Define boundary conditions on elements, displacement, stress or traction at element centres, friction etc. This then calculates the slip due to these conditions and draws slip distributions or allows you to animate this as shown in the gif above.
-
Step 4. Define observation points in the elastic medium. Defining points at which you would like to evaluate stress and displacement (or secondary fault surfaces to look at Coulomb stress change on faces).
-
Step 5. Calculate stress and displacement at the observation points.
-
Step 6. Compute error on observation points due to the infinitesimal strain assumption.
-
Step 7. Visualise and analyse stress at the observation points, draw coulomb stress change, principal stresses, distortion maps, iso-contours etc.
Its recommended that you look in the analytical test directories for more information on setting up complex problems.
This part of the document goes through the notation used in the code so you can look up variables you see while reading the code.
S = σ
E = ε
U = υ
T = t
P = Ρ
x & X = x direction
y & Y = y direction
z & Z = z direction
r & R = radial distance/radial direction when a tensor component subscript
t * T = angle θ from y-axis, circumferential direction when a tensor component subscript
p * P = angle φ from z-axis, circumferential direction when a tensor component subscript
Shear modulus (G in some texts):
mu = μLamés constant:
lambda = λPoisson's ratio:
nu = νYoungs modulus:
E = Ε
[ Sxx Sxy ]
[ Sxy Syy ]
[ Sxx Sxy Sxz ]
[ Sxy Syy Syz ]
[ Sxz Syz Szz ]
[ Exx Exy ]
[ Exy Eyy ]
[ Exx Exy Exz ]
[ Exy Eyy Eyz ]
[ Exz Eyz Ezz ]
(for 2D scripts ignore the z components and shear components will be a single s)
Diagram of angles ax ay az and lengths Nx Ny Nz
[ CosAx ]
[ CosAy ]
[ CosAz ]
Diagram of tractions Tx Ty Tz. These components are converted to normal and shear traction components in relation to the element orientation.
3D 2D
[ Tnn ] [ Tn ]
[ Tss ] [ Ts ]
[ Tds ]
Diagram of displacement Dn (a and b) and Ds (d and c) of an element.
3D 2D
[ Dn ] [ Dn ]
[ Dss ] [ Ds ]
[ Dds ]
Diagram of principal stresses around a magma chamber in 2D.
3D 2D
[ S1 ] [ S1 ]
[ S2 ] [ S2 ]
[ S3 ]
3D 2D
[ E1 ] [ E1 ]
[ E2 ] [ E2 ]
[ E3 ]
3D 2D
[ S1Dir ] [ S1Dir ]
[ S2Dir ] [ S2Dir ]
[ S3Dir ]
Matrix A: Table containing stress values of how much a discontinuity of magnitude 1 from the jth element (table’s column number) effects the stress at the midpoint of every other element (i) (table’s row number).
[ A ]
In 2d will be composed as:
[ DnTn DsTn ]
[ DnTs DsTs ]
Where DnTs represents element j opening by a magnitude of 1 and its effect on the shear traction at element i.
Matrix B: Table containing displacement values of how much a displacement of 1 from the jth element causes the midpoint of every other element (i) to displace.
[ B ]
Matrix C: Inverse of matrix A, the summed influence of all elements displacing the amount described in each column of matrix [C] will cause a traction of one stress unit, t on element (i) in the direction defined by the subscript.
[ C ]
[ I ]
infmatrix - Stress
Dinfmatrix - Displacement
StressInf - Stress
DispInf - Displacement
NUM = Number of elements
sz = Dimension of a vector
dimx = Number of cols of matrix;
dimy = Number of rows of matrix;
dimz = Size of matrix in 3rd dimension;
Triangles = Triangles is a list where each row contains 3 row index locations in array 'Points' which contains the XYZ location of each corner of the triangle
Points = Cols 2 3 and 4 are the xyz locations of one the corner points of a triangle.
MidPoint = Each row is the xyz location of the midpoint of each triangle
P1 P2 P3 = Each row is the corner point of one triangle (xyz). Each correponding row in P1 P2 P3 are the 3 corner points of a triangle.
FaceNormalVector = Direction cosines of each tris normal. CosAx is simply "CosAx=FaceNormalVector(:,1);"
Lines = Lines is a list where each row contains 2 row index locations in array 'Points' which contains the xy location of each lines end points.
Points = Cols 2 3 and 4 are the xyz locations of one the corner points of a triangle.
MidPoint = Each row is the xy location of the midpoint of each triangle
P1 P2 = Each row is the end point of one line in (xy). Each correponding row in P1 P2 is the end points of a line.
LineNormalVector = Direction cosines of each lines normal. CosAx is simply "CosAx=LineNormalVector(:,1);"
Halflength = The half length of each 2D element
Cohesion = Cohesive strength of material
CSC = Coulomb stress change
Mu = Coefficient of friction (not to be confused with shear modulus mu)
E1 = Elastic in body 1 (Ω1)
FB = Free boundary, no continity conditions apply at these elements
IF = Interface, traction and displacement must be equal at the two coincident elements at the interface between two bodies.
halfspace = Compute coefficients and stress/disps in a half-space.
Fdisp = Elements midpoints that have a displacement set to 0. I.e this point in the body is locked.
bad = Typically used to remove elements or nan elements. I.e. "X(bad)=[];"
Uniform = Flag to say if the points are on a uniformly spaced grid.
SecondSurface = Flag that says there is a second surface we compute CSC on.
StressTTotal = 3(2d) or 6(3D) vectors for tensors due to remote stresses and movement of the dislocations.
StressTChg = 3(2d) or 6(3D) tensor at each point due to the movement of dislocations.
StressTReg = 3(2d) or 6(3D) tensor at each point that is the remote stress in the body (stress at infinity).
The basic boundary elements are from:
2D: Crouch, S.L. and Starfield, A.M., 1982. Boundary element methods in solid mechanics: with applications in rock mechanics and geological engineering. Allen & Unwin.
-
Frictional effects across fault surfaces, both coefficient of friction and cohesion. Using the formulations described in:
Ritz, E., Mutlu, O. and Pollard, D.D., 2012. Integrating complementarity into the 2D displacement discontinuity boundary element method to model faults and fractures with frictional contact properties. Computers & Geosciences, 45, pp.304-312. Kaven, J.O., Hickman, S.H., Davatzes, N.C. and Mutlu, O., 2012. Linear complementarity formulation for 3D frictional sliding problems. Computational Geosciences, 16(3), pp.613-624.
-
Gravitational effects of topography as described in:
-
Elemental locking as described in 2D by (useful for magma chamber modelling):
Crouch, S.L. and Starfield, A.M., 1982. Boundary element methods in solid mechanics: with applications in rock mechanics and geological engineering. Allen & Unwin.
-
Inhomogeneous elastics (interface between two types): Latter sections of:
Crouch, S.L. and Starfield, A.M., 1982. Boundary element methods in solid mechanics: with applications in rock mechanics and geological engineering. Allen & Unwin.
-
Calculating finite strain (and error): Equations are given in:
Pollard, D.D. and Fletcher, R.C., 2005. Fundamentals of structural geology. Cambridge University Press.
This uses two calculations. For evenly gridded data finite strain is calculated using a modified version of: Dirk-Jan Kroon's MATLAB script 'finite strain' and the finite strain calculation of: Cardozo, N. and Allmendinger, R.W., 2009. SSPX: A program to compute strain from displacement/velocity data. Computers & Geosciences, 35(6), pp.1343-1357. -
Calculating stress intensity factors (2D and 3D): Equations are given in:
Thanks to:
Prof Eleonora Rivalta, University of Bologna/GFZ (PhD supervisor).
Dr David Healy, University of Aberdeen (MSc supervisor).
Dr Mehdi Nikkhoo, GFZ.
Dr Juliet Crider, University of Washington.
Dr Robert Simpson, Glasgow University.
And others who have replied to emails regarding code, solutions etc.
Questions? Contact "timothy.davis[at]bristol.ac.uk"