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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 107 - (hide annotations)
Thu Jan 27 06:21:48 2005 UTC (14 years, 5 months ago) by jgs
Original Path: trunk/esys2/doc/user/wave.tex
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 jgs 107 % $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