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

Revision 3348 - (show annotations)
Fri Nov 12 05:55:50 2010 UTC (9 years ago) by caltinay
File MIME type: application/x-tex
File size: 7619 byte(s)
Added weipa import statements to the relevant parts in user guide and cookbook.
Documented the deprecation of esys.escript.saveVTK in the changes chapter.


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