/[escript]/trunk/doc/user/slip.tex
ViewVC logotype

Contents of /trunk/doc/user/slip.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (3 years ago) by jfenwick
File MIME type: application/x-tex
File size: 8206 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2018 by The University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Apache License, version 2.0
8 % http://www.apache.org/licenses/LICENSE-2.0
9 %
10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development 2012-2013 by School of Earth Sciences
12 % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 %
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15
16 \section{Slip on a Fault}\label{Slip CHAP}
17 \begin{figure}[ht]
18 \centerline{\includegraphics{Slip1}}
19 \caption{Domain $\Omega=[0,1]^2$ with a vertical fault of length $0.5$.}
20 \label{fig:slip.1}
21 \end{figure}
22 %
23 In this example we illustrate how to calculate the stress distribution around
24 a fault\index{fault} in the Earth's crust caused by a slip\index{slip} through
25 an earthquake.
26
27 To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with
28 a vertical fault in its center as illustrated in \fig{fig:slip.1}.
29 We assume that the slip distribution $s_{i}$ on the fault is known.
30 We want to calculate the distribution of the displacements $u_{i}$\index{displacement}
31 and stress $\sigma_{ij}$\index{stress} in the domain.
32 Further, we assume an isotropic, linear elastic material model of the form
33 \begin{eqnarray} \label{Slip stress}
34 \sigma_{ij} & = & \lambda u_{k,k} \delta_{ij} + \mu ( u_{i,j} + u_{j,i})
35 \end{eqnarray}
36 where $\lambda$ and $\mu$ are the Lam\'e coefficients\index{Lam\'e coefficients}
37 and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
38 On the boundary the normal stress is given by
39 \begin{eqnarray} \label{Slip natural fault}
40 \sigma_{ij}n_{j}=0
41 \end{eqnarray}
42 and normal displacements are set to zero:
43 \begin{eqnarray} \label{Slip constraint}
44 u_{i}n_{i} =0
45 \end{eqnarray}
46 The stress needs to fulfill the momentum equation\index{momentum equation}
47 \begin{eqnarray}\label{Slip general problem}
48 - \sigma_{ij,j}=0
49 \end{eqnarray}
50 This problem is very similar to the elastic deformation problem presented in \Sec{ELASTIC CHAP}.
51 However, we need to address an additional challenge: the displacement
52 $u_{i}$ is in fact discontinuous across the fault, but we are in the
53 lucky situation that we know the jump of the displacements across the fault.
54 This is in fact the given slip $s_{i}$.
55 So we can split the total distribution $u_{i}$ into a component
56 $v_{i}$ which is continuous across the fault and the known slip $s_{i}$
57 \begin{eqnarray}\label{Slip Split}
58 u_{i} = v_{i} + \frac{1}{2} s^{\pm}_{i}
59 \end{eqnarray}
60 where $s^{\pm}=s$ when right of the fault and $s^{\pm}=-s$ when left of the fault.
61 We assume that $s^{\pm}=0$ when sufficiently away from the fault.
62
63 We insert this into the stress definition in \eqn{Slip stress}
64 \begin{eqnarray} \label{Slip stress split}
65 \sigma_{ij} & = &
66 \sigma^c_{ij} +
67 \frac{1}{2} \sigma^s_{ij}
68 \end{eqnarray}
69 with
70 \begin{eqnarray} \label{Slip stress split 1}
71 \sigma^c_{ij} = \lambda v_{k,k} \delta_{ij} + \mu ( v_{i,j} + v_{j,i})
72 \end{eqnarray}
73 and
74 \begin{eqnarray} \label{Slip stress split 2}
75 \sigma^s_{ij} = \lambda s^{\pm}_{k,k} \delta_{ij} + \mu ( s^{\pm}_{i,j} + s^{\pm}_{j,i}).
76 \end{eqnarray}
77 In fact, $\sigma^s_{ij}$ defines a stress jump across the fault.
78 An easy way to construct this function is to use a function $\chi$ which is
79 $1$ on the right and $-1$ on the left side from the fault.
80 One can then set
81 \begin{eqnarray} \label{Slip stress split 23 }
82 \sigma^s_{ij} = \chi \cdot ( \lambda s_{k,k} \delta_{ij} + \mu ( s_{i,j} + s_{j,i}) )
83 \end{eqnarray}
84 assuming that $s$ is extended by zero away from the fault.
85 After inserting \eqn{Slip stress split} into (\ref{Slip general problem}) we
86 get the differential equation\index{momentum equation}
87 \begin{eqnarray}\label{Slip general problem 2 }
88 - \sigma^c_{ij,j}=\frac{1}{2} \sigma^s_{ij,j}
89 \end{eqnarray}
90 Together with the definition (\ref{Slip stress split 1}) we have a
91 differential equation for the continuous function $v_i$.
92 Notice that the boundary condition (\ref{Slip constraint}) and (\ref{Slip natural fault})
93 transfer to $v_i$ and $\sigma^c_{ij}$ as $s$ is zero away from the fault.
94 In \Sec{ELASTIC CHAP} we have discussed how this problem is solved using
95 the \LinearPDE class. We refer to this section for further details.
96
97 To define the fault we use the \class{FaultSystem} class introduced in \Sec{Fault System}.
98 The following statements define a fault system \var{fs} and add the fault \var{1} to the system:
99 \begin{python}
100 fs=FaultSystem(dim=2)
101 fs.addFault(fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
102 \end{python}
103 The fault added starts at point $(0.5,0.25)$ has length $0.5$ and points north.
104 The main purpose of the \class{FaultSystem} class is to define a
105 parameterization of the fault using a local coordinate system.
106 One can inquire the class to get the range used to parameterize a fault.
107 \begin{python}
108 p0,p1 = fs.getW0Range(tag=1)
109 \end{python}
110 Typically \var{p0} is equal to zero while \var{p1} is equal to the length of the fault.
111 The parameterization is given as a mapping from a set of local coordinates
112 onto a parameter range (in our case the range \var{p0} to \var{p1}).
113 For instance, to map the entire domain \var{mydomain} onto the fault one can
114 use
115 \begin{python}
116 x = mydomain.getX()
117 p,m = fs.getParametrization(x, tag=1)
118 \end{python}
119 Of course there is the problem that not all locations are on the fault.
120 For those locations which are on the fault \var{m} is set to 1, otherwise 0 is used.
121 So on return the values of \var{p} define the value of the fault parameterization
122 (typically the distance from the starting point of the fault along the fault)
123 where \var{m} is positive.
124 On all other locations the value of \var{p} is undefined.
125 Now \var{p} can be used to define a slip distribution on the fault via
126 \begin{python}
127 s = m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
128 \end{python}
129 Notice the factor \var{m} which ensures that \var{s} is zero away from the fault.
130 It is important that the slip is zero at the ends of the faults.
131
132 We can now put all components together to get the script:
133 \begin{python}
134 from esys.escript import *
135 from esys.escript.linearPDEs import LinearPDE
136 from esys.escript.models import FaultSystem
137 from esys.finley import Rectangle
138 from esys.weipa import saveVTK
139 from esys.escript.unitsSI import DEG
140
141 #... set some parameters ...
142 lam=1.
143 mu=1
144 slip_max=1.
145 mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16) # n1 needs to be a multiple of 4!
146 # .. create the fault system
147 fs=FaultSystem(dim=2)
148 fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
149 # ... create a slip distribution on the fault
150 p, m=fs.getParametrization(mydomain.getX(), tag=1)
151 p0,p1= fs.getW0Range(tag=1)
152 s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
153 # ... calculate stress according to slip:
154 D=symmetric(grad(s))
155 chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(), tag=1)
156 sigma_s=(mu*D+lam*trace(D)*kronecker(mydomain))*chi
157 #... open symmetric PDE ...
158 mypde=LinearPDE(mydomain)
159 mypde.setSymmetryOn()
160 #... set coefficients ...
161 C=Tensor4(0., Function(mydomain))
162 for i in range(mydomain.getDim()):
163 for j in range(mydomain.getDim()):
164 C[i,i,j,j]+=lam
165 C[j,i,j,i]+=mu
166 C[j,i,i,j]+=mu
167 # ... fix displacement in normal direction
168 x=mydomain.getX()
169 msk=whereZero(x[0])*[1.,0.] + whereZero(x[0]-1.)*[1.,0.] \
170 +whereZero(x[1])*[0.,1.] + whereZero(x[1]-1.)*[0.,1.]
171 mypde.setValue(A=C, X=-0.5*sigma_s, q=msk)
172 #... solve pde ...
173 mypde.getSolverOptions().setVerbosityOn()
174 v=mypde.getSolution()
175 # .. write the displacement to file:
176 D=symmetric(grad(v))
177 sigma=(mu*D+lam*trace(D)*kronecker(mydomain))+0.5*sigma_s
178 saveVTK("slip.vtu", disp=v+0.5*chi*s, stress=sigma)
179 \end{python}
180 The script creates the file \file{slip.vtu} which contains the total
181 displacements and stress.
182 These values are stored as cell-centered data.
183 %
184 \begin{figure} [ht]
185 \centerline{\includegraphics[width=\figwidth]{Slip2}}
186 \caption{Total Displacement after the slip event}
187 \label{fig:slip.2}
188 \end{figure}
189 %
190 See \fig{fig:slip.2} for a visualization of the result.
191

  ViewVC Help
Powered by ViewVC 1.1.26