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

Revision 1316 - (show annotations)
Tue Sep 25 03:18:30 2007 UTC (13 years, 2 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 % 2 % $Id$ 3 % 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5 % 6 % 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 \section{3D Wave Propagation} 18 \label{WAVE CHAP} 19 20 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 \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 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 \begin{eqnarray} \label{WAVE initial} 43 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 \end{eqnarray} 45 for all $x$ in the domain. 46 47 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 51 We use an explicit time integration scheme to calculate the displacement field $u$ at 52 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 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 56 \end{eqnarray} 57 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 58 \begin{eqnarray} \label{WAVE dyn 2b} 59 u\hackscore{,tt}=G(u) 60 \end{eqnarray} 61 where one sets $a^{(n)}=G(u^{(n-1)})$. 62 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 \begin{eqnarray} \label{WAVE natural at n} 69 \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 \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 \end{eqnarray} 75 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 79 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will 80 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 81 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is 82 \begin{equation}\label{WAVE dyn 100} 83 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 84 \end{equation} 85 The natural boundary condition 86 \begin{equation}\label{WAVE dyn 101} 87 n\hackscore{j}X\hackscore{ij} =0 88 \end{equation} 89 is used. 90 91 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 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 \end{array} 97 \end{equation} 98 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 103 The following script defines a the function \function{wavePropagation} which 104 implements the Verlet scheme to solve our wave propagation problem. 105 The argument \var{domain} which is a \Domain class object 106 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 \var{rho} are material properties. 109 \begin{python} 110 from esys.linearPDEs import LinearPDE 111 from numarray import identity,zeros,ones 112 from esys.escript import * 113 from esys.escript.pdetools import Locator 114 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 115 x=domain.getX() 116 # ... open new PDE ... 117 mypde=LinearPDE(domain) 118 mypde.setSolverMethod(LinearPDE.LUMPING) 119 kronecker=identity(mypde.getDim()) 120 # 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 # ... set initial values .... 129 n=0 130 # 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 t=0 135 # 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 while t

## Properties

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