| 1 |
jgs |
107 |
% $Id$ |
| 2 |
jgs |
121 |
\section{3D Wave Propagation} |
| 3 |
jgs |
107 |
\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 |
lkettle |
575 |
u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ |
| 48 |
jgs |
107 |
\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 |
gross |
565 |
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through |
| 77 |
jgs |
107 |
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} |
Parent Directory
| 
Revision Log