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

Revision 593 - (show annotations)
Tue Mar 14 02:50:27 2006 UTC (14 years, 3 months ago) by gross
File MIME type: application/x-tex
File size: 15053 byte(s)
updates on escript documentation (unfinished)

 1 % $Id$ 2 \section{3D Wave Propagation} 3 \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 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 \begin{eqnarray} \label{WAVE initial} 28 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 \end{eqnarray} 30 for all $x$ in the domain. 31 32 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 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 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 41 \end{eqnarray} 42 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 43 \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 at n} 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 \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 \end{eqnarray} 60 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 64 In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will 65 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 66 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The forms which are relevant here are 67 \begin{equation}\label{WAVE dyn 100} 68 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 69 \end{equation} 70 Implicitly, the natural boundary condition 71 \begin{equation}\label{WAVE dyn 101} 72 n\hackscore{j}X\hackscore{ij} =0 73 \end{equation} 74 is assumed. The \LinearPDE object allows defining constraints of the form 75 \begin{equation}\label{WAVE dyn 102} 76 u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0 77 \end{equation} 78 where $r$ and $q$ are given functions. However in this problem we don't need to define 79 any constraints for our problem. 80 81 With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: 82 \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 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 93 The following script defines a the function \function{wavePropagation} which 94 implements the Verlet scheme to solve our wave propagation problem. 95 The argument \var{domain} which is a \Domain class object 96 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 \begin{python} 100 from esys.linearPDEs import LinearPDE 101 from numarray import identity,zeros,ones 102 from esys.escript import * 103 from esys.escript.pdetools import Locator 104 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 105 x=domain.getX() 106 # ... open new PDE ... 107 mypde=LinearPDE(domain) 108 mypde.setSolverMethod(LinearPDE.LUMPING) 109 kronecker=identity(mypde.getDim()) 110 # 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 # ... set initial values .... 119 n=0 120 # 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 t=0 125 # 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 while t

## Properties

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