Home

Awesome

CutAndDisplace

References to cite

DOI

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

Index

FaultsDisplacingUnderShearStress
Gif of a fault array shearing due to a remote stress. Computed using the 3D code

Citation

Davis, T., 2017. A new open source boundary element code and its application to geological deformation: Exploring stress concentrations around voids and the effects of 3D frictional distributions on fault surfaces (M.Sc thesis. Aberdeen University)

Introduction

Boundary Element MATLAB code. Modelling faults and deformation

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'.

1. How to run

Index ^

Codes for 2D and 3D are laid out in the same manner. To run:

2. Analytical tests

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:

2.1. 2D Analytical tests

Test by running the file runallregressions.m in 2D/AnalyticalTests folder or the individual scripts further inside this directory.

2.2. 3D Analytical tests

Test by running the file runallregressions.m in 3D/AnalyticalTests folder or the individual scripts further inside this directory.

3. Files and data format

Index ^

Formats that the code uses:

3.1. 2D files and data

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.

3.2. 3D files and data

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.

3.3. Code structure

Index ^

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:

Its recommended that you look in the analytical test directories for more information on setting up complex problems.

4. Notation

Index ^

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.

4.1. Basics

Stress:

S = σ

Strain:

E = ε

Displacement:

U = υ

Traction:

T = t

Pressure:

P = Ρ

Cartesian coordinates:

x & X = x direction

y & Y = y direction

z & Z = z direction

Radial, cylindrical or spherical coordinates:

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

Elastic constants:

Shear modulus (G in some texts):
mu = μ

Lamés constant:
lambda = λ

Poisson's ratio:
nu = ν

Youngs modulus:
E = Ε

4.2. Tensors

2*2 2D stress tensor

Diagram of 2D stress tensor

    [ Sxx Sxy ]    
    [ Sxy Syy ]

3*3 3D stress tensor

Diagram of 3D stress tensor

    [ Sxx Sxy Sxz ]    
    [ Sxy Syy Syz ]    
    [ Sxz Syz Szz ]

2*2 2D strain tensor

    [ Exx Exy ]    
    [ Exy Eyy ]

3*3 3D strain tensor

    [ Exx Exy Exz ]    
    [ Exy Eyy Eyz ]    
    [ Exz Eyz Ezz ]

4.3. Vectors

(for 2D scripts ignore the z components and shear components will be a single s)

Direction cosines of surfaces normal, often known as: nx ny nz.

Diagram of angles ax ay az and lengths Nx Ny Nz

    [ CosAx ]   
    [ CosAy ]     
    [ CosAz ]

Traction magnitude in the surface normal, strike-slip and dip-slip directions. Tn and Ts for 2D

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 ]

Discontinity or burgers vector movement in the surface normal, strike-slip and dip-slip directions

Diagram of displacement Dn (a and b) and Ds (d and c) of an element.

      3D                       2D
    [ Dn ]                   [ Dn ]   
    [ Dss ]                  [ Ds ]       
    [ Dds ]

Principal stresses magnitudes (convention is that S1 is the most tensile stress)

Diagram of principal stresses around a magma chamber in 2D.

      3D                       2D
    [ S1 ]                   [ S1 ]      
    [ S2 ]                   [ S2 ]    
    [ S3 ]

Principal strain

      3D                       2D
    [ E1 ]                   [ E1 ]       
    [ E2 ]                   [ E2 ]      
    [ E3 ]

Principal directions (stress)

      3D                       2D
    [ S1Dir ]                [ S1Dir ]      
    [ S2Dir ]                [ S2Dir ]    
    [ S3Dir ]

4.4. Matrices

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 ]

Matrix I: Square matrix with ones on the main diagonal

    [ I ]

Influence matrix: Matrix of zeros to be filled with influence coefficients.

    infmatrix       - Stress  
    Dinfmatrix      - Displacement

Structure array containing influence matricies

    StressInf       - Stress  
    DispInf         - Displacement

4.5. Element geometry

Sizes:

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;

Triangulation (3D):

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);"

Line parts (2D):

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

Mechanics terms:

Cohesion = Cohesive strength of material

CSC = Coulomb stress change

Mu = Coefficient of friction (not to be confused with shear modulus mu)

Multiple elastic bodies (inhomogenous material interfaces):

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.

4.6. Flags

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.

4.7. Tables

Tensor component tables (strain tensors are the same but with the word straina the front).

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).

5. Elements used in 2D and 3D

Index ^

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.

3D: Nikkhoo, M. and Walter, T.R., 2015. Triangular dislocation: an analytical, artefact-free solution. Geophysical Journal International, 201(2), pp.1117-1139.

6. Code highlights and applications

Index ^

7. Acknowledgements

Index ^

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"