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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 11 months 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 caltinay 5293
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland
4 caltinay 5293 % http://www.uq.edu.au
5     %
6     % Primary Business: Queensland, Australia
7 jfenwick 6112 % Licensed under the Apache License, version 2.0
8     % http://www.apache.org/licenses/LICENSE-2.0
9 caltinay 5293 %
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 caltinay 3291 \section{Slip on a Fault}\label{Slip CHAP}
17     \begin{figure}[ht]
18     \centerline{\includegraphics{Slip1}}
19 gross 2654 \caption{Domain $\Omega=[0,1]^2$ with a vertical fault of length $0.5$.}
20     \label{fig:slip.1}
21     \end{figure}
22 caltinay 3291 %
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 gross 2654
27     To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with
28 caltinay 3291 a vertical fault in its center as illustrated in \fig{fig:slip.1}.
29 jfenwick 3295 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 caltinay 3291 and stress $\sigma_{ij}$\index{stress} in the domain.
32     Further, we assume an isotropic, linear elastic material model of the form
33 gross 2654 \begin{eqnarray} \label{Slip stress}
34 jfenwick 3295 \sigma_{ij} & = & \lambda u_{k,k} \delta_{ij} + \mu ( u_{i,j} + u_{j,i})
35 gross 2647 \end{eqnarray}
36 sshaw 5284 where $\lambda$ and $\mu$ are the Lam\'e coefficients\index{Lam\'e coefficients}
37 jfenwick 3295 and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
38 gross 2647 On the boundary the normal stress is given by
39 gross 2654 \begin{eqnarray} \label{Slip natural fault}
40 jfenwick 3295 \sigma_{ij}n_{j}=0
41 gross 2647 \end{eqnarray}
42 gross 2654 and normal displacements are set to zero:
43     \begin{eqnarray} \label{Slip constraint}
44 jfenwick 3295 u_{i}n_{i} =0
45 gross 2654 \end{eqnarray}
46 caltinay 3291 The stress needs to fulfill the momentum equation\index{momentum equation}
47 gross 2654 \begin{eqnarray}\label{Slip general problem}
48 jfenwick 3295 - \sigma_{ij,j}=0
49 gross 2654 \end{eqnarray}
50 caltinay 3291 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 jfenwick 3295 $u_{i}$ is in fact discontinuous across the fault, but we are in the
53 caltinay 3291 lucky situation that we know the jump of the displacements across the fault.
54 jfenwick 3295 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 gross 2654 \begin{eqnarray}\label{Slip Split}
58 jfenwick 3295 u_{i} = v_{i} + \frac{1}{2} s^{\pm}_{i}
59 gross 2654 \end{eqnarray}
60 caltinay 3291 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 gross 2654
63 caltinay 3291 We insert this into the stress definition in \eqn{Slip stress}
64     \begin{eqnarray} \label{Slip stress split}
65 jfenwick 3295 \sigma_{ij} & = &
66     \sigma^c_{ij} +
67     \frac{1}{2} \sigma^s_{ij}
68 gross 2654 \end{eqnarray}
69 caltinay 3291 with
70     \begin{eqnarray} \label{Slip stress split 1}
71 jfenwick 3295 \sigma^c_{ij} = \lambda v_{k,k} \delta_{ij} + \mu ( v_{i,j} + v_{j,i})
72 gross 2654 \end{eqnarray}
73 caltinay 3291 and
74     \begin{eqnarray} \label{Slip stress split 2}
75 jfenwick 3295 \sigma^s_{ij} = \lambda s^{\pm}_{k,k} \delta_{ij} + \mu ( s^{\pm}_{i,j} + s^{\pm}_{j,i}).
76 gross 2654 \end{eqnarray}
77 jfenwick 3295 In fact, $\sigma^s_{ij}$ defines a stress jump across the fault.
78 caltinay 3291 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 gross 2654 \begin{eqnarray} \label{Slip stress split 23 }
82 jfenwick 3295 \sigma^s_{ij} = \chi \cdot ( \lambda s_{k,k} \delta_{ij} + \mu ( s_{i,j} + s_{j,i}) )
83 gross 2654 \end{eqnarray}
84 caltinay 3291 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 gross 2654 \begin{eqnarray}\label{Slip general problem 2 }
88 jfenwick 3295 - \sigma^c_{ij,j}=\frac{1}{2} \sigma^s_{ij,j}
89 gross 2654 \end{eqnarray}
90 caltinay 3291 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 jfenwick 3295 transfer to $v_i$ and $\sigma^c_{ij}$ as $s$ is zero away from the fault.
94 caltinay 3291 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 gross 2654
97 caltinay 3291 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 gross 2654 \begin{python}
100 caltinay 3291 fs=FaultSystem(dim=2)
101     fs.addFault(fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
102 gross 2654 \end{python}
103 caltinay 3331 The fault added starts at point $(0.5,0.25)$ has length $0.5$ and points north.
104 caltinay 3291 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 gross 2654 \begin{python}
108 caltinay 3291 p0,p1 = fs.getW0Range(tag=1)
109 gross 2654 \end{python}
110 caltinay 3291 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 gross 2654 \begin{python}
116 caltinay 3291 x = mydomain.getX()
117     p,m = fs.getParametrization(x, tag=1)
118 gross 2654 \end{python}
119 caltinay 3291 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 gross 2654 \begin{python}
127 caltinay 3291 s = m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
128 gross 2654 \end{python}
129 caltinay 3291 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 gross 2654
132 caltinay 3291 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 caltinay 3348 from esys.weipa import saveVTK
139 caltinay 3291 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 gross 2654 \begin{figure} [ht]
185 caltinay 3279 \centerline{\includegraphics[width=\figwidth]{Slip2}}
186 caltinay 3291 \caption{Total Displacement after the slip event}
187 gross 2654 \label{fig:slip.2}
188     \end{figure}
189 caltinay 3291 %
190     See \fig{fig:slip.2} for a visualization of the result.
191 gross 2654

  ViewVC Help
Powered by ViewVC 1.1.26