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

Revision 660 - (hide annotations)
Fri Mar 24 08:43:21 2006 UTC (17 years ago) by gross
File MIME type: application/x-tex
File size: 15028 byte(s)
that does not look to bad now although more stuff could be added.

 1 jgs 107 % $Id$ 2 gross 625 % 3 % Copyright © 2006 by ACcESS MNRF 4 % \url{http://www.access.edu.au 5 % Primary Business: Queensland, Australia. 6 % Licensed under the Open Software License version 3.0 7 8 % 9 jgs 121 \section{3D Wave Propagation} 10 jgs 107 \label{WAVE CHAP} 11 12 We want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation 13 \index{wave equation}: 14 \begin{eqnarray}\label{WAVE general problem} 15 \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0 16 \end{eqnarray} 17 in a three dimensional block of length $L$ in $x\hackscore{0}$ 18 and $x\hackscore{1}$ direction and height $H$ 19 in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location. 20 $\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by 21 \begin{eqnarray} \label{WAVE stress} 22 \sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 23 \end{eqnarray} 24 where $\lambda$ and $\mu$ are the Lame coefficients 25 \index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. 26 On the boundary the normal stress is given by 27 \begin{eqnarray} \label{WAVE natural} 28 \sigma\hackscore{ij}n\hackscore{j}=0 29 \end{eqnarray} 30 lkettle 581 for all time $t>0$. 31 32 At initial time $t=0$ the displacement 33 $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: 34 jgs 107 \begin{eqnarray} \label{WAVE initial} 35 lkettle 581 u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array} \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; 36 jgs 107 \end{eqnarray} 37 for all $x$ in the domain. 38 39 lkettle 581 Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we 40 set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point 41 $x\hackscore C$. 42 jgs 107 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 lkettle 581 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 50 jgs 107 \begin{eqnarray} \label{WAVE dyn 2b} 51 u\hackscore{,tt}=G(u) 52 \end{eqnarray} 53 gross 660 where one sets $a^{(n)}=G(u^{(n-1)})$. 54 jgs 107 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 gross 593 \begin{eqnarray} \label{WAVE natural at n} 61 jgs 107 \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 lkettle 581 \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 jgs 107 \end{eqnarray} 67 lkettle 581 Now we have converted our problem for displacement, $u^{(n)}$, into a problem for 68 acceleration, $a^(n)$, which now depends 69 on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. 70 jgs 107 71 In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will 72 gross 565 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 73 gross 625 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here are 74 jgs 107 \begin{equation}\label{WAVE dyn 100} 75 lkettle 581 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 76 jgs 107 \end{equation} 77 jgs 110 Implicitly, the natural boundary condition 78 jgs 107 \begin{equation}\label{WAVE dyn 101} 79 n\hackscore{j}X\hackscore{ij} =0 80 \end{equation} 81 gross 625 is assumed. 82 jgs 107 83 lkettle 581 With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: 84 jgs 107 \begin{equation}\label{WAVE dyn 6} 85 \begin{array}{ll} 86 D\hackscore{ij}=\rho \delta\hackscore{ij}& 87 X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ 88 \end{array} 89 \end{equation} 90 lkettle 581 91 92 %Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero 93 %(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. 94 jgs 107 95 lkettle 581 The following script defines a the function \function{wavePropagation} which 96 jgs 107 implements the Verlet scheme to solve our wave propagation problem. 97 The argument \var{domain} which is a \Domain class object 98 lkettle 581 defines the domain of the problem. \var{h} and \var{tend} are the step size 99 and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and 100 \var{rho} are material properties. 101 jgs 107 \begin{python} 102 from esys.linearPDEs import LinearPDE 103 lkettle 581 from numarray import identity,zeros,ones 104 jgs 107 from esys.escript import * 105 lkettle 581 from esys.escript.pdetools import Locator 106 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 107 jgs 107 x=domain.getX() 108 # ... open new PDE ... 109 mypde=LinearPDE(domain) 110 lkettle 581 mypde.setSolverMethod(LinearPDE.LUMPING) 111 jgs 110 kronecker=identity(mypde.getDim()) 112 lkettle 581 # spherical source at middle of bottom face 113 xc=[width/2.,width/2.,0.] 114 # define small radius around point xc 115 # Lsup(x) returns the maximum value of the argument x 116 src_radius = 0.1*Lsup(domain.getSize()) 117 print "src_radius = ",src_radius 118 dunit=numarray.array([1.,0.,0.]) # defines direction of point source 119 mypde.setValue(D=kronecker*rho) 120 jgs 107 # ... set initial values .... 121 jgs 110 n=0 122 lkettle 581 # initial value of displacement at point source is constant (U0=0.01) 123 # for first two time steps 124 u=U0*whereNegative(length(x-xc)-src_radius)*dunit 125 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit 126 jgs 107 t=0 127 lkettle 581 # define the location of the point source 128 L=Locator(domain,xc) 129 # find potential at point source 130 u_pc=L.getValue(u) 131 print "u at point charge=",u_pc 132 u_pc_x = u_pc[0] 133 u_pc_y = u_pc[1] 134 u_pc_z = u_pc[2] 135 # open file to save displacement at point source 136 u_pc_data=open('./data/U_pc.out','w') 137 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 138 jgs 107 while t

## Properties

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