| 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 |
|
Parent Directory
| 
Revision Log