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

Revision 2104 - (show annotations)
Thu Nov 27 05:36:47 2008 UTC (13 years, 7 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 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2008 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 % 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 14 15 \section{3-D Wave Propagation} 16 \label{WAVE CHAP} 17 18 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 \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 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 \begin{eqnarray} \label{WAVE initial} 41 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 \end{eqnarray} 43 for all $x$ in the domain. 44 45 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 49 We use an explicit time integration scheme to calculate the displacement field $u$ at 50 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 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 54 \end{eqnarray} 55 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 56 \begin{eqnarray} \label{WAVE dyn 2b} 57 u\hackscore{,tt}=G(u) 58 \end{eqnarray} 59 where one sets $a^{(n)}=G(u^{(n-1)})$. 60 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 \begin{eqnarray} \label{WAVE natural at n} 67 \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 \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 \end{eqnarray} 73 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 77 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will 78 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 79 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is 80 \begin{equation}\label{WAVE dyn 100} 81 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 82 \end{equation} 83 The natural boundary condition 84 \begin{equation}\label{WAVE dyn 101} 85 n\hackscore{j}X\hackscore{ij} =0 86 \end{equation} 87 is used. 88 89 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 90 \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 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 101 The following script defines a the function \function{wavePropagation} which 102 implements the Verlet scheme to solve our wave propagation problem. 103 The argument \var{domain} which is a \Domain class object 104 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 \var{rho} are material properties. 107 \begin{python} 108 from esys.linearPDEs import LinearPDE 109 from numarray import identity,zeros,ones 110 from esys.escript import * 111 from esys.escript.pdetools import Locator 112 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 113 x=domain.getX() 114 # ... open new PDE ... 115 mypde=LinearPDE(domain) 116 mypde.setSolverMethod(LinearPDE.LUMPING) 117 kronecker=identity(mypde.getDim()) 118 # 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 # ... set initial values .... 127 n=0 128 # 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 t=0 133 # 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 139 u_pc_y = u_pc 140 u_pc_z = u_pc 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 while t

## Properties

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