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

Revision 2663 - (hide annotations)
Tue Sep 15 06:04:44 2009 UTC (11 years, 4 months ago) by gross
File MIME type: application/x-tex
File size: 7882 byte(s)
finally there is a 3D version for the fault system class

 1 gross 2654 \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} 6 7 gross 2647 \section{Slip on a Fault} 8 \label{Slip CHAP} 9 10 gross 2654 In this example we illustrate how to calculate the stress distribution around a fault \index{fault} in the 11 Earth's crust caused by a slip \index{slip} through an earthquake. 12 13 To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with 14 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 \begin{eqnarray} \label{Slip stress} 16 gross 2647 \sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 17 \end{eqnarray} 18 where $\lambda$ and $\mu$ are the Lame coefficients 19 \index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. 20 On the boundary the normal stress is given by 21 gross 2654 \begin{eqnarray} \label{Slip natural fault} 22 gross 2647 \sigma\hackscore{ij}n\hackscore{j}=0 23 \end{eqnarray} 24 gross 2654 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}. 34 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. 44 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. 75 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 fs=FaultSystem(dim=2) 81 gross 2663 fs.addFault(fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1) 82 gross 2654 \end{python} 83 gross 2663 The fault added starts at point $(0.5,0.25)$ has length $0.5$ and is pointing north. 84 gross 2654 The main purpose of the \class{FaultSystem} class is to define a parameterization 85 of the fault using a local coordinate system. One can enquire the class 86 to get the range used to parameterize a fault. 87 \begin{python} 88 p0,p1= fs.getW0Range(tag=1) 89 \end{python} 90 Typically \var{p0} is equal to zero while \var{p1} is equal to the length of the fault. The parameterization 91 is given as a mapping from a set of local coordinates onto parameter range (in our case 92 the range \var{p0} to \var{p1}). For instance to map the entire domain \var{mydomain} anto the fault 93 one can use 94 \begin{python} 95 x=mydomain.getX() 96 p, m=fs.getParametrization(x,tag=1) 97 \end{python} 98 Of course there is the problem that not all location are on the fault. For those locations which are on the 99 fault \var{m} is set to $1$ otherwise $0$ is used. So on return the values of \var{p} are defining 100 the value of the fault parameterization (typically the distance from the starting point of the fault along the fault) 101 where \var{m} is positive. On all other locations the value if \var{p} is undefined. Now \var{p} can 102 be used to define a slip distribution on the fault via 103 \begin{python} 104 s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.] 105 \end{python} 106 Notice that the factor \var{m} which makes sure that \var{s} is zero away from the fault. It is important 107 that the slip is zero at the ends of the faults. 108 109 \begin{figure} [ht] 110 \centerline{\includegraphics[width=\figwidth]{figures/Slip2}} 111 \caption{Total Displacement after slip event.} 112 \label{fig:slip.2} 113 \end{figure} 114 115 We can now put all components together to get the script: 116 \begin{python} 117 from esys.escript import * 118 from esys.escript.linearPDEs import LinearPDE 119 from esys.escript.models import FaultSystem 120 from esys.finley import Rectangle 121 gross 2663 from esys.escript.unitsSI import DEG 122 123 gross 2654 #... set some parameters ... 124 lam=1. 125 mu=1 126 slip_max=1. 127 mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16) # n1 need to be multiple of 4!!! 128 # .. create the fault system 129 fs=FaultSystem(dim=2) 130 gross 2663 fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1) 131 gross 2654 # ... create a slip distribution on the fault: 132 p, m=fs.getParametrization(mydomain.getX(),tag=1) 133 p0,p1= fs.getW0Range(tag=1) 134 s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.] 135 # ... calculate stress according to slip: 136 D=symmetric(grad(s)) 137 chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(),tag=1) 138 sigma_s=(mu*D+lam*trace(D)*kronecker(mydomain))*chi 139 #... open symmetric PDE ... 140 mypde=LinearPDE(mydomain) 141 mypde.setSymmetryOn() 142 #... set coefficients ... 143 C=Tensor4(0.,Function(mydomain)) 144 for i in range(mydomain.getDim()): 145 for j in range(mydomain.getDim()): 146 C[i,i,j,j]+=lam 147 C[j,i,j,i]+=mu 148 C[j,i,i,j]+=mu 149 # ... fix displacement in normal direction 150 x=mydomain.getX() 151 msk=whereZero(x)*[1.,0.] + whereZero(x-1.)*[1.,0.] \ 152 +whereZero(x)*[0.,1.] + whereZero(x-1.)*[0.,1.] 153 mypde.setValue(A=C,X=-0.5*sigma_s,q=msk) 154 #... solve pde ... 155 mypde.getSolverOptions().setVerbosityOn() 156 v=mypde.getSolution() 157 # .. write the displacement to file: 158 D=symmetric(grad(v)) 159 sigma=(mu*D+lam*trace(D)*kronecker(mydomain))+0.5*sigma_s 160 saveVTK("slip.vtu",disp=v+0.5*chi*s, stress= sigma) 161 \end{python} 162 The script creates the \file{slip.vtu} giving the total displacements and stress. These values are stored as cell 163 centered data. See Figure~\ref{fig:slip.2} for the result.