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

Revision 2104 - (hide annotations)
Thu Nov 27 05:36:47 2008 UTC (13 years, 8 months ago) by jfenwick
File MIME type: application/x-tex
File size: 14972 byte(s)
Fixed non-compilage.
Removed mysterious t and t! from output.

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

## Properties

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