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

Revision 1316 - (hide annotations)
Tue Sep 25 03:18:30 2007 UTC (12 years, 10 months ago) by ksteube
File MIME type: application/x-tex
File size: 14990 byte(s)
Quickly edited chapters 1 and 2 of the User Guide, but it needs more work.
Ran entire document through spell checker.


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

## Properties

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