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

Revision 581 - (hide annotations)
Wed Mar 8 05:48:51 2006 UTC (15 years, 7 months ago) by lkettle
File MIME type: application/x-tex
File size: 15046 byte(s)
changes to section 2.3 on 3D wave propagation in online reference guide and added some figures


 1 jgs 107 % $Id$ 2 jgs 121 \section{3D Wave Propagation} 3 jgs 107 \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 lkettle 581 for all time $t>0$. 24 25 At initial time $t=0$ the displacement 26 $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: 27 jgs 107 \begin{eqnarray} \label{WAVE initial} 28 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 \; 29 jgs 107 \end{eqnarray} 30 for all $x$ in the domain. 31 32 lkettle 581 Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we 33 set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point 34 $x\hackscore C$. 35 jgs 107 36 We use an explicit time-integration scheme to calculate the displacement field $u$ at 37 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$ 38 which is defined by 39 \begin{eqnarray} \label{WAVE dyn 2} 40 lkettle 575 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 41 jgs 107 \end{eqnarray} 42 lkettle 581 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 43 jgs 107 \begin{eqnarray} \label{WAVE dyn 2b} 44 u\hackscore{,tt}=G(u) 45 \end{eqnarray} 46 where one sets $a^{(n)}=G(u^{(n-1)})$, see \Ref{Mora94}. 47 48 In our case $a^{(n)}$ is given by 49 \begin{eqnarray}\label{WAVE dyn 3} 50 \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j} 51 \end{eqnarray} 52 and boundary conditions 53 \begin{eqnarray} \label{WAVE natural} 54 \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0 55 \end{eqnarray} 56 derived from \eqn{WAVE natural} where 57 \begin{eqnarray} \label{WAVE dyn 3a} 58 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}). 59 jgs 107 \end{eqnarray} 60 lkettle 581 Now we have converted our problem for displacement, $u^{(n)}$, into a problem for 61 acceleration, $a^(n)$, which now depends 62 on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. 63 jgs 107 64 In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will 65 gross 565 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 66 lkettle 581 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The forms which are relevant here are 67 jgs 107 \begin{equation}\label{WAVE dyn 100} 68 lkettle 581 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 69 jgs 107 \end{equation} 70 jgs 110 Implicitly, the natural boundary condition 71 jgs 107 \begin{equation}\label{WAVE dyn 101} 72 n\hackscore{j}X\hackscore{ij} =0 73 \end{equation} 74 jgs 110 is assumed. The \LinearPDE object allows defining constraints of the form 75 jgs 107 \begin{equation}\label{WAVE dyn 101} 76 u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0 77 \end{equation} 78 lkettle 581 where $r$ and $q$ are given functions. However in this problem we don't need to define 79 any constraints for our problem. 80 jgs 107 81 lkettle 581 With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: 82 jgs 107 \begin{equation}\label{WAVE dyn 6} 83 \begin{array}{ll} 84 D\hackscore{ij}=\rho \delta\hackscore{ij}& 85 X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ 86 \end{array} 87 \end{equation} 88 lkettle 581 89 90 %Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero 91 %(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. 92 jgs 107 93 lkettle 581 The following script defines a the function \function{wavePropagation} which 94 jgs 107 implements the Verlet scheme to solve our wave propagation problem. 95 The argument \var{domain} which is a \Domain class object 96 lkettle 581 defines the domain of the problem. \var{h} and \var{tend} are the step size 97 and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and 98 \var{rho} are material properties. 99 jgs 107 \begin{python} 100 from esys.linearPDEs import LinearPDE 101 lkettle 581 from numarray import identity,zeros,ones 102 jgs 107 from esys.escript import * 103 lkettle 581 from esys.escript.pdetools import Locator 104 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 105 jgs 107 x=domain.getX() 106 # ... open new PDE ... 107 mypde=LinearPDE(domain) 108 lkettle 581 mypde.setSolverMethod(LinearPDE.LUMPING) 109 jgs 110 kronecker=identity(mypde.getDim()) 110 lkettle 581 # spherical source at middle of bottom face 111 xc=[width/2.,width/2.,0.] 112 # define small radius around point xc 113 # Lsup(x) returns the maximum value of the argument x 114 src_radius = 0.1*Lsup(domain.getSize()) 115 print "src_radius = ",src_radius 116 dunit=numarray.array([1.,0.,0.]) # defines direction of point source 117 mypde.setValue(D=kronecker*rho) 118 jgs 107 # ... set initial values .... 119 jgs 110 n=0 120 lkettle 581 # initial value of displacement at point source is constant (U0=0.01) 121 # for first two time steps 122 u=U0*whereNegative(length(x-xc)-src_radius)*dunit 123 u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit 124 jgs 107 t=0 125 lkettle 581 # define the location of the point source 126 L=Locator(domain,xc) 127 # find potential at point source 128 u_pc=L.getValue(u) 129 print "u at point charge=",u_pc 130 u_pc_x = u_pc 131 u_pc_y = u_pc 132 u_pc_z = u_pc 133 # open file to save displacement at point source 134 u_pc_data=open('./data/U_pc.out','w') 135 u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 136 jgs 107 while t

## Properties

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