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

Revision 575 - (show annotations)
Fri Mar 3 03:33:07 2006 UTC (15 years, 6 months ago) by lkettle
File MIME type: application/x-tex
File size: 12827 byte(s)
I have changed some of the documentation and added more explanations for
the online reference guide for esys13. I have modified two of the
example source codes to write out the results for Helmholtz problem and
changed one variable name in the diffusion.py code to avoid confusion.


 1 % $Id$ 2 \section{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-2)} + 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 75 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 \linearPDEs 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 D\hackscore{ij} u\hackscore{j} = - X\hackscore{ij,j}\; . 80 \end{equation} 81 Implicitly, 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 constraints 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 function 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 # ... open new PDE ... 115 mypde=LinearPDE(domain) 116 mypde.setLumpingOn() 117 kronecker=identity(mypde.getDim()) 118 mypde.setValue(D=kronecker*rho, \ 119 q=x.whereZero()*kronecker[1,:]) 120 # ... set initial values .... 121 n=0 122 u=Vector(0,ContinuousFunction(domain)) 123 u_last=Vector(0,ContinuousFunction(domain)) 124 t=0 125 while t

## Properties

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