ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2651 by jfenwick, Mon Sep 7 03:39:45 2009 UTC revision 2654 by gross, Tue Sep 8 07:11:12 2009 UTC
# Line 1  Line 1 
1    \begin{figure} [ht]
2    \centerline{\includegraphics[width=\figwidth]{figures/Slip1}}
3    \caption{Domain $\Omega=[0,1]^2$ with a vertical fault of length $0.5$.}
4    \label{fig:slip.1}
5    \end{figure}
7  \section{Slip on a Fault}  \section{Slip on a Fault}
8  \label{Slip CHAP}  \label{Slip CHAP}
10  In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation:  In this example we illustrate how to calculate the stress distribution around a fault \index{fault} in the
11  \index{wave equation}  Earth's crust caused by a slip \index{slip} through an earthquake.
12  \begin{eqnarray}\label{WAVE general problem fault}  
13  \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0  To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with
14  \end{eqnarray}  a vertical fault in its center. We assume that the slip distribution $s\hackscore{i}$ on the fault is known. We want to calculate the distribution of the displacements $u\hackscore{i}$\index{displacement} and stress $\sigma_{ij}$\index{stress} in the domain. We assume an isotropic, linear elastic material model of the form
15  in a three dimensional block of length $L$ in $x\hackscore{0}$  \begin{eqnarray} \label{Slip  stress}
 and $x\hackscore{1}$ direction and height $H$  
 in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location.  
 $\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by  
 \begin{eqnarray} \label{WAVE stress fault}  
16  \sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})  \sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
17  \end{eqnarray}  \end{eqnarray}
18  where $\lambda$ and $\mu$ are the Lame coefficients  where $\lambda$ and $\mu$ are the Lame coefficients
19  \index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.  \index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
20  On the boundary the normal stress is given by  On the boundary the normal stress is given by
21  \begin{eqnarray} \label{WAVE natural fault}  \begin{eqnarray} \label{Slip natural fault}
22  \sigma\hackscore{ij}n\hackscore{j}=0  \sigma\hackscore{ij}n\hackscore{j}=0
23  \end{eqnarray}  \end{eqnarray}
 for all time $t>0$.  
24    and normal displacements are set to zero:
25    \begin{eqnarray} \label{Slip constraint}
26    u\hackscore{i}n\hackscore{i} =0
27    \end{eqnarray}
28    The stress needs to fulfill the momentum equation
29    \index{momentum equation}
30    \begin{eqnarray}\label{Slip general problem}
31    - \sigma\hackscore{ij,j}=0
32    \end{eqnarray}
33    This problem is very similar to the elastic deformation problem presented in Section~\ref{ELASTIC CHAP}.
35    However, we need to address an additional challenge: The displacement $u\hackscore{i}$ is in fact discontinuous across the fault
36    but we are in the lucky situation that we know the jump of the displacements across the fault. This is
37    in fact the the given slip $s\hackscore{i}$. So we can split the total distribution $u\hackscore{i}$ into a component $v\hackscore{i}$
38    which is continuous across the fault and the known slip $s\hackscore{i}$
39    \begin{eqnarray}\label{Slip Split}
40    u\hackscore{i} = v\hackscore{i} + \frac{1}{2} s^{\pm}\hackscore{i}
41    \end{eqnarray}
42    where $s^{\pm}=s$ when right of the fault and $s^{\pm}=-s$
43    when left of the fault. We assume that $s^{\pm}=0$ sufficiently away from the fault.
45    We insert this into the stress definition~\ref{Slip  stress}
46    \begin{eqnarray} \label{Slip  stress split}
47    \sigma\hackscore{ij} & = &
48    \sigma^c\hackscore{ij} +
49    \frac{1}{2} \sigma^s\hackscore{ij}
50    \end{eqnarray}
51     with
52    \begin{eqnarray} \label{Slip  stress split 1 }
53    \sigma^c\hackscore{ij} = \lambda v\hackscore{k,k} \delta\hackscore{ij} + \mu ( v\hackscore{i,j} + v\hackscore{j,i})
54    \end{eqnarray}
55    and
56    \begin{eqnarray} \label{Slip  stress split 2 }
57    \sigma^s\hackscore{ij} = \lambda s^{\pm}\hackscore{k,k} \delta\hackscore{ij} + \mu ( s^{\pm}\hackscore{i,j} + s^{\pm}\hackscore{j,i})
58    \end{eqnarray}
59    In fact $\sigma^s\hackscore{ij}$ defines a stress jump across the fault. In easy way to construct this
60    function is to us a function $\chi$ which is $1$ on the right and $-1$ on the left side from the fault. One can then
61    set
62    \begin{eqnarray} \label{Slip  stress split 23 }
63    \sigma^s\hackscore{ij} = \chi \cdot  ( \lambda s\hackscore{k,k} \delta\hackscore{ij} + \mu ( s\hackscore{i,j} + s\hackscore{j,i}) )
64    \end{eqnarray}
65    assuming that $s$ is extended by zero away from the fault. After inserting~\ref{Slip  stress split} into~\ref{Slip general problem} we get the differential equation
66    \index{momentum equation}
67    \begin{eqnarray}\label{Slip general problem 2 }
68    - \sigma^c\hackscore{ij,j}=\frac{1}{2} \sigma^s\hackscore{ij,j}
69    \end{eqnarray}
70    Together with the definition~\ref{Slip  stress split 1} we have a differential equation for the
71    continuous function $v_i$. Notice that the boundary condition~\ref{Slip constraint}
72    and~\ref{Slip natural fault} transfer to $v_i$ and $\sigma^c\hackscore{ij}$ as
73    $s$ is zero away from the fault. In Section~\ref{ELASTIC CHAP} we have discussed how this problem
74    is solved using the \LinearPDE class. We refer to this section for further details.
76    To define the fault we use the \class{FaultSystem} class introduced in Section~\ref{Fault System}.
77    The following statements define a fault system \var{fs} and add the fault \var{1}
78    to the system.
79    \begin{python}
80    fault_start=[0.5,0.25]
81    fault_end=[0.5,0.75]
82    fs=FaultSystem(dim=2)
83    fs.addFault(top=[fault_start, fault_end], tag=1)
84    \end{python}
85    The main purpose of the \class{FaultSystem} class is to define a parameterization
86    of the fault using a local coordinate system. One can enquire the class
87    to get the range used to parameterize a fault.
88    \begin{python}
89    p0,p1= fs.getW0Range(tag=1)
90    \end{python}
91    Typically \var{p0} is equal to zero while \var{p1} is equal to the length of the fault. The parameterization
92    is given as a mapping from a set of local coordinates onto parameter range (in our case
93    the range \var{p0} to \var{p1}). For instance to map the entire domain \var{mydomain} anto the fault
94    one can use
95    \begin{python}
96    x=mydomain.getX()
97    p, m=fs.getParametrization(x,tag=1)
98    \end{python}
99    Of course there is the problem that not all location are on the fault. For those locations which are on the
100    fault \var{m} is set to $1$ otherwise $0$ is used. So on return the values of \var{p} are defining
101    the value of the fault parameterization (typically the distance from the starting point of the fault along the fault)
102    where \var{m} is positive. On all other locations the value if \var{p} is undefined. Now \var{p} can
103    be used to define a slip distribution on the fault via
104    \begin{python}
105    s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
106    \end{python}
107    Notice that the factor \var{m} which makes sure that \var{s} is zero away from the fault. It is important
108    that the slip is zero at the ends of the faults.
110    \begin{figure} [ht]
111    \centerline{\includegraphics[width=\figwidth]{figures/Slip2}}
112    \caption{Total Displacement after slip event.}
113    \label{fig:slip.2}
114    \end{figure}
116    We can now put all components together to get the script:
117    \begin{python}
118    from esys.escript import *
119    from esys.escript.linearPDEs import LinearPDE
120    from esys.escript.models import FaultSystem
121    from esys.finley import Rectangle
122    #... set some parameters ...
123    lam=1.
124    mu=1
125    slip_max=1.
126    fault_start=[0.5,0.25]
127    fault_end=[0.5,0.75]
129    mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16)  # n1 need to be multiple of 4!!!
130    # .. create the fault system
131    fs=FaultSystem(dim=2)
132    fs.addFault(top=[fault_start, fault_end], tag=1)
133    # ... create a slip distribution on the fault:
134    p, m=fs.getParametrization(mydomain.getX(),tag=1)
135    p0,p1= fs.getW0Range(tag=1)
136    s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
137    # ... calculate stress according to slip:
138    D=symmetric(grad(s))
139    chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(),tag=1)
140    sigma_s=(mu*D+lam*trace(D)*kronecker(mydomain))*chi
141    #... open symmetric PDE ...
142    mypde=LinearPDE(mydomain)
143    mypde.setSymmetryOn()
144    #... set coefficients ...
145    C=Tensor4(0.,Function(mydomain))
146    for i in range(mydomain.getDim()):
147      for j in range(mydomain.getDim()):
148         C[i,i,j,j]+=lam
149         C[j,i,j,i]+=mu
150         C[j,i,i,j]+=mu
151    # ... fix displacement in normal direction
152    x=mydomain.getX()
153    msk=whereZero(x[0])*[1.,0.] + whereZero(x[0]-1.)*[1.,0.] \
154       +whereZero(x[1])*[0.,1.] + whereZero(x[1]-1.)*[0.,1.]
155    mypde.setValue(A=C,X=-0.5*sigma_s,q=msk)
156    #... solve pde ...
157    mypde.getSolverOptions().setVerbosityOn()
158    v=mypde.getSolution()
159    # .. write the displacement to file:
160    D=symmetric(grad(v))
161    sigma=(mu*D+lam*trace(D)*kronecker(mydomain))+0.5*sigma_s
162    saveVTK("slip.vtu",disp=v+0.5*chi*s, stress= sigma)
163    \end{python}
164    The script creates the \file{slip.vtu} giving the total displacements and stress. These values are stored as cell
165    centered data. See Figure~\ref{fig:slip.2} for the result.

Removed from v.2651  
changed lines
  Added in v.2654

  ViewVC Help
Powered by ViewVC 1.1.26