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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 110 - (hide annotations)
Mon Feb 14 04:14:42 2005 UTC (14 years, 6 months ago) by jgs
File MIME type: application/x-tex
File size: 12831 byte(s)
*** empty log message ***

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 jgs 110
75 jgs 107 In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will
76     use the \LinearPDE class of the \linearPDEsPack to do so. The general form of the PDE defined through
77     the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are
78     \begin{equation}\label{WAVE dyn 100}
79 jgs 110 D\hackscore{ij} u\hackscore{j} = - X\hackscore{ij,j}\; .
80 jgs 107 \end{equation}
81 jgs 110 Implicitly, the natural boundary condition
82 jgs 107 \begin{equation}\label{WAVE dyn 101}
83     n\hackscore{j}X\hackscore{ij} =0
84     \end{equation}
85 jgs 110 is assumed. The \LinearPDE object allows defining constraints of the form
86 jgs 107 \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 jgs 110 \var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time:
108 jgs 107 \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     # ... open new PDE ...
115     mypde=LinearPDE(domain)
116     mypde.setLumpingOn()
117 jgs 110 kronecker=identity(mypde.getDim())
118 jgs 107 mypde.setValue(D=kronecker*rho, \
119 jgs 110 q=x[0].whereZero()*kronecker[1,:])
120 jgs 107 # ... set initial values ....
121 jgs 110 n=0
122 jgs 107 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 jgs 110 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
129     # ... get new acceleration ....
130     mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])
131 jgs 107 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     t+=h
138 jgs 110 n+=1
139     print n,"-th time step t ",t
140     print "a=",inf(a),sup(a)
141     print "u=",inf(u),sup(u)
142     # ... save current acceleration in units of gravity
143     if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10))
144 jgs 107 \end{python}
145 jgs 110 Notice that
146     all coefficients of the PDE which are independent from time $t$ are set outside the \code{while}
147     loop. This allows the \LinearPDE class to ruse information as much as possible
148     when iterating over time.
149    
150     We have seen most of the functions in previous examples but there are some new functions here:
151 jgs 107 \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
152     It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance
153     for \finley the statement \code{grad(grad(u))} will raise an exception.
154     \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
155     and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie.
156     $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
157    
158     The initialization of \var{u} and \var{u_last} needs some explanation. The statement
159     \begin{python}
160 jgs 110 u=Vector(0.,ContinuousFunction(domain))
161 jgs 107 \end{python}
162 jgs 110 declares \var{u} as a vector-valued function of its coordinates
163     in the domain (i.e. with two or three components in the case of
164     a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal
165     to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}
166     specifies the `smoothness` of the function. In our case we want the displacement field to be
167     a continuous function over the \Domain \var{domain}. On other option would be
168     \var{Function(domain)} in which case continuity is not garanteed. In fact the
169     argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
170     returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.
171 jgs 107
172 jgs 110 The statement \code{myPDE.setLumpingOn()}
173 jgs 107 switches on the usage of an aggressive approximation of the PDE operator, in this case
174     formed by the coefficient \var{D} (actually the discrete
175     version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of
176     an additional error, very fast
177     solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var
178     {C} \code{setLumpingOn} will produce wrong results.
179 jgs 110
180 jgs 107 It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent
181     with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time
182     direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the
183     constraint defined by \eqn{WAVE impact} fulfills
184     the initial conditions in \eqn{WAVE initial}.
185    
186    
187 jgs 110 One of the big advantages of the Verlet scheme is the fact that the problem to be solved
188     in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping).
189     The problem becomes so simple because we use the stress from the last time step rather then the stress which
190     actually is present at the current time step. Schemes using this approach are called an explicit time integration
191     scheme \index{explicit scheme} \index{time integration!explicit}. The
192     backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
193     an example of an implicit schemes
194     \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
195     each variable at a particular time rather then values from previous time steps. This will lead to problem which is
196     more expensive to solve, in particular for non-linear problems.
197     Although
198     explicit time integration schemes are cheep to finalize a single time step, they need a significant smaller time
199     step then implicit schemes and can suffer from stability problems. Therefor they need a
200     very careful selection of the time step size $h$.
201 jgs 107
202 jgs 110 An easy, heuristic way of choosing an appropriate
203 jgs 107 time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
204     which says that within a time step a information should not travel further than a cell used in the
205     discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
206     $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
207     \begin{eqnarray}\label{WAVE dyn 66}
208     h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
209     \end{eqnarray}
210     where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristic of
211     the formula.
212    
213     The following script uses the \code{wavePropagtion} function to solve the
214 jgs 110 wave equation over a block in crust of the earth. The depth of the block is
215     $10km$ and the width in each direction is $100km$. The depth of the block is alligned
216     with the $x\hackscore{0}$-direction. \var{ne} gives the number of elements in $x\hackscore{0}$-direction. The number
217     of elements in $x\hackscore{1}$- and $x\hackscore{2}$-direction is chosen such that the elements
218     become cubic. The impact is $2m$ acting within a time of about $10sec$:
219 jgs 107 \begin{python}
220 jgs 110 ne=10 # number of cells in x_0-direction
221     depth=10000. # length in x_0-direction
222     width=100000. # length in x_1 and x_2 direction
223     lam=3.462e9
224     mu=3.462e9
225 jgs 107 rho=1154.
226 jgs 110 tau=10.
227     umax=2.
228     tend=60
229     h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)
230     def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)
231    
232     print "time step size = ",h
233     mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
234     wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
235 jgs 107 \end{python}
236     The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}.
237    
238    
239     % \begin{figure}
240     % \centerline{\includegraphics[width=\figwidth]{DiffusionDomain}}
241     % \caption{Temperature Diffusion Problem with Circular Heat Source}
242     % \label{WAVE FIG 1}
243     % \end{figure}
244    
245 jgs 110 \begin{figure}
246 jgs 107 % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
247     % \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
248     % \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
249     %\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
250 jgs 110 \caption{Selected time steps of a wave propagation over a $10km \times 100km \times 100km$
251     with time step size $h=0.0666$. Color represents the acceleration.}
252     \label{WAVE FIG 2}
253     \end{figure}
254 jgs 107
255    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26