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

Revision 625 - (show annotations)
Thu Mar 23 00:41:25 2006 UTC (13 years, 5 months ago) by gross
File MIME type: application/x-tex
File size: 15047 byte(s)
some updates and linearPDE.tex switched off

 1 % $Id$ 2 % 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 \section{3D Wave Propagation} 10 \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 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 \begin{eqnarray} \label{WAVE initial} 35 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 \end{eqnarray} 37 for all $x$ in the domain. 38 39 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 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 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 48 \end{eqnarray} 49 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 50 \begin{eqnarray} \label{WAVE dyn 2b} 51 u\hackscore{,tt}=G(u) 52 \end{eqnarray} 53 where one sets $a^{(n)}=G(u^{(n-1)})$, see \Ref{Mora94}. 54 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 \begin{eqnarray} \label{WAVE natural at n} 61 \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 \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 \end{eqnarray} 67 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 71 In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will 72 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 73 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here are 74 \begin{equation}\label{WAVE dyn 100} 75 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 76 \end{equation} 77 Implicitly, the natural boundary condition 78 \begin{equation}\label{WAVE dyn 101} 79 n\hackscore{j}X\hackscore{ij} =0 80 \end{equation} 81 is assumed. 82 83 With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: 84 \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 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 95 The following script defines a the function \function{wavePropagation} which 96 implements the Verlet scheme to solve our wave propagation problem. 97 The argument \var{domain} which is a \Domain class object 98 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 \begin{python} 102 from esys.linearPDEs import LinearPDE 103 from numarray import identity,zeros,ones 104 from esys.escript import * 105 from esys.escript.pdetools import Locator 106 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 107 x=domain.getX() 108 # ... open new PDE ... 109 mypde=LinearPDE(domain) 110 mypde.setSolverMethod(LinearPDE.LUMPING) 111 kronecker=identity(mypde.getDim()) 112 # 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 # ... set initial values .... 121 n=0 122 # 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 t=0 127 # 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 while t

## Properties

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