/[escript]/trunk/esys2/doc/user/wave.tex
ViewVC logotype

Contents of /trunk/esys2/doc/user/wave.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 107 - (show annotations)
Thu Jan 27 06:21:48 2005 UTC (14 years, 6 months ago) by jgs
File MIME type: application/x-tex
File size: 11512 byte(s)
commit of branch dev-01 back to main trunk on 2005-01-27

1 % $Id$
2 \chapter{3D Wave Propagation}
3 \label{WAVE CHAP}
4
5 We want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation
6 \index{wave equation}:
7 \begin{eqnarray}\label{WAVE general problem}
8 \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0
9 \end{eqnarray}
10 in a three dimensional block of length $L$ in $x\hackscore{0}$
11 and $x\hackscore{1}$ direction and height $H$
12 in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location.
13 $\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by
14 \begin{eqnarray} \label{WAVE stress}
15 \sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
16 \end{eqnarray}
17 where $\lambda$ and $\mu$ are the Lame coefficients
18 \index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
19 On the boundary the normal stress is given by
20 \begin{eqnarray} \label{WAVE natural}
21 \sigma\hackscore{ij}n\hackscore{j}=0
22 \end{eqnarray}
23 for all time $t>0$. At initial time $t=0$ the displacement
24 $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are give:
25 \begin{eqnarray} \label{WAVE initial}
26 u\hackscore{i}(0,x)=0 & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \;
27 \end{eqnarray}
28 for all $x$ in the domain.
29
30 The points on the left face of the domain, ie. points $x=(x\hackscore{i})$ with $x\hackscore{0}=0$,
31 are moved into the $x\hackscore{2}$ direction by $u^{max}$. This is expressed by the constraint
32 \begin{eqnarray} \label{WAVE impact}
33 u\hackscore{2}(x,t)=s(t) \mbox{ for } x=(x\hackscore{i}) with x\hackscore{0}=0
34 \end{eqnarray}
35 where
36 \begin{eqnarray} \label{WAVE impact s}
37 s(t)=u^{max} (1-e^{-(\tau^{-1}t)^3})
38 \end{eqnarray}
39 $\tau$ measures the time needed to reach the maximum displacement $u^{max}$ (
40 actually $\tau (ln 2)^{\frac{1}{3}}) \approx 0.9 \tau$ is the time until $\frac{u^{max}}{2}$ is reached).
41 The smaller $\tau$ the faster $u^{max}$ is reached.
42
43 We use an explicit time-integration scheme to calculate the displacement field $u$ at
44 certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$
45 which is defined by
46 \begin{eqnarray} \label{WAVE dyn 2}
47 u^{(n)}=2u^{(n-1)}-u^{(n)} + h^2 a^{(n)} \\
48 \end{eqnarray}
49 for all $n=1,2,\ldots$. It is designed to solve a system of equations of the form
50 \begin{eqnarray} \label{WAVE dyn 2b}
51 u\hackscore{,tt}=G(u)
52 \end{eqnarray}
53 where one sets $a^{(n)}=G(u^{(n-1)})$, see \Ref{Mora94}.
54
55 In our case $a^{(n)}$ is given by
56 \begin{eqnarray}\label{WAVE dyn 3}
57 \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}
58 \end{eqnarray}
59 and boundary conditions
60 \begin{eqnarray} \label{WAVE natural}
61 \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0
62 \end{eqnarray}
63 derived from \eqn{WAVE natural} where
64 \begin{eqnarray} \label{WAVE dyn 3a}
65 \sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i})
66 \end{eqnarray}.
67 Additionally $a^{(n)}$ has to fullfill the constraint
68 \begin{eqnarray}\label{WAVE dyn 4}
69 a^{(n)}\hackscore{2}(x)= s\hackscore{,tt}(t^{(n)})\mbox{ where } s\hackscore{,tt}(t)=\frac{u^{max}}{\tau^2}
70 (6\tau^{-1}t- 9(\tau^{-1}t)^4) e^{-(\tau^{-1}t)^3}
71 \end{eqnarray}
72 derived from \eqn{WAVE impact}.
73
74 In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will
75 use the \LinearPDE class of the \linearPDEsPack to do so. The general form of the PDE defined through
76 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are
77 \begin{equation}\label{WAVE dyn 100}
78 D\hacksco
79 re{ij} u\hackscore{j} = - X\hackscore{ij,j}\; .
80 \end{equation}
81 Impliciltly, the natural boundary condition
82 \begin{equation}\label{WAVE dyn 101}
83 n\hackscore{j}X\hackscore{ij} =0
84 \end{equation}
85 is assumed. The \LinearPDE object allows defining constriants of the form
86 \begin{equation}\label{WAVE dyn 101}
87 u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0
88 \end{equation}
89 where $r$ and $q$ are given functions.
90
91 With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$, $r$ and $q$:
92 \begin{equation}\label{WAVE dyn 6}
93 \begin{array}{ll}
94 D\hackscore{ij}=\rho \delta\hackscore{ij}&
95 X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
96 r\hackscore{i}=\delta\hackscore{i2} \cdot s\hackscore{,tt}(t^{(n)}) & q\hackscore{i}=\delta\hackscore{i2} \cdot x\hackscore{0}.\texttt{whereZero()}
97 \end{array}
98 \end{equation}
99 Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero
100 (up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero.
101
102 In the following script defines a the function \function{wavePropagation} which
103 implements the Verlet scheme to solve our wave propagation problem.
104 The argument \var{domain} which is a \Domain class object
105 defines the domain of the problem. \var{h} and \var{tend} is the step size
106 and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and
107 \var{rho} are material properties. \var{s_tt} is a fucntion which returns $s\hackscore{,tt}$ at a given time:
108 \begin{python}
109 from esys.linearPDEs import LinearPDE
110 from numarray import identity
111 from esys.escript import *
112 def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
113 x=domain.getX()
114 ndim=domain.getDim()
115 # ... open new PDE ...
116 mypde=LinearPDE(domain)
117 mypde.setLumpingOn()
118 kronecker=identity(ndim)
119 mypde.setValue(D=kronecker*rho, \
120 q=x[0].whereZero()*kronecker[ndim-1,:])
121 # ... set initial values ....
122 u=Vector(0,ContinuousFunction(domain))
123 u_last=Vector(0,ContinuousFunction(domain))
124 t=0
125 while t<tend:
126 # ... get current stress ....
127 g=grad(u)
128 stress=trace(g)*lam*kronecker+mu*(g+transpose(g))
129 # ... get new acceleration ....
130 mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[ndim-1,:])
131 a=mypde.getSolution()
132 # ... get new displacement ...
133 u_new=2*u-u_last+h**2*a
134 # ... shift displacements ....
135 u_last=u
136 u=u_new
137 # ... next time step ...
138 t+=h
139 \end{python}
140 We have seen most of the functions in previous examples. However,
141 there are some new functions here:
142 \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
143 It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance
144 for \finley the statement \code{grad(grad(u))} will raise an exception.
145 \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
146 and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie.
147 $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
148
149 The initialization of \var{u} and \var{u_last} needs some explanation. The statement
150 \begin{python}
151 u=Vector(0,ContinuousFunction(domain))
152 \end{python}
153 creates a new object which is of class \Data, see \Sec{SEC ESCRIPT DATA}. The object represnts a vector-valued
154 function of the location coordinates. Vector-value means that it has two or three components on a two or three dimensional domain. respectively. The second argument The first argument (here $=0$) is the initial value
155 assigned to any spatial point.
156
157
158
159 All coefficients of the PDE which are independent from time are set outside the \code{while}
160 loop. This allows the \LinearPDE class to ruse information as much as possible
161 when iterating over time. The statement \code{myPDE.setLumpingOn()}
162 switches on the usage of an aggressive approximation of the PDE operator, in this case
163 formed by the coefficient \var{D} (actually the discrete
164 version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of
165 an additional error, very fast
166 solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var
167 {C} \code{setLumpingOn} will produce wrong results.
168 =========
169 It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent
170 with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time
171 direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the
172 constraint defined by \eqn{WAVE impact} fulfills
173 the initial conditions in \eqn{WAVE initial}.
174
175
176 ================
177 The advantage of the Verlat scheme is that the PDE we have to solve in each time step is very
178 simple as we are using the stress based on the displacements of last time step. One could, similar to the
179 backward Euler scheme we have used in \Chap{DIFFUSION CHAP}, use the stress calculated with current
180 displacement which would lead to a more complex PDE involving a coefficient \var{A}. These so-called implicit schemes
181 \index{implicit scheme} \index{time integration!implicit} is more expensive in each time then the Verlat scheme which
182 is a so called explicit scheme \index{explicit scheme} \index{time integration!explicit}.
183
184 Explicit time integration schemes need a very careful selection of the time step size $h$ otherwise
185 if $h$ is too big the scheme can be unstable. An easy, heuristic way of choosing an appropriate
186 time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
187 which says that within a time step a information should not travel further than a cell used in the
188 discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
189 $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
190 \begin{eqnarray}\label{WAVE dyn 66}
191 h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
192 \end{eqnarray}
193 where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristic of
194 the formula.
195
196 The following script uses the \code{wavePropagtion} function to solve the
197 wave equation over a block in crust of the earth. The
198 depth is $l\hackscore{0}=20km$ and the length is $l\hackscore{1}=l\hackscore{2}=200km$.
199 The time of the impact is about $2sec$ and maximum displacement is $u^{max}=15m$ which corresponds to an heavy
200 earthquake. We are using a $6 \times 60 \times 60$ mesh:
201 \begin{python}
202 ne=6
203 lame_lambda=3.462e9
204 lame_mu=3.462e9
205 rho=1154.
206 tau=2.
207 umax=15.
208 xc=[0,1000,1000]
209 r=1.
210 tend=10.
211 dt=1./5.*sqrt(rho/(lame_lambda+2*lame_mu))(20000./ne)
212 print "step size = ",dt
213 mydomain=finely.Brick(ne,10*ne,10*ne,l0=20000,l1=200000,l2=200000)
214 wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax)
215 \end{python}
216 The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}.
217
218
219 % \begin{figure}
220 % \centerline{\includegraphics[width=\figwidth]{DiffusionDomain}}
221 % \caption{Temperature Diffusion Problem with Circular Heat Source}
222 % \label{WAVE FIG 1}
223 % \end{figure}
224
225 % \begin{figure}
226 % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
227 % \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
228 % \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
229 %\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
230 % \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
231 % \label{WAVE FIG 2}
232 %\end{figure}
233
234

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26