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

Revision 2503 - (show annotations)
Wed Jul 1 05:05:08 2009 UTC (9 years, 9 months ago) by jfenwick
File MIME type: application/x-tex
File size: 15859 byte(s)
Fixing doc build so examples tar and zip are built properly.
Fixed image filenames in user guide.

Please do not put extensions on image files.


 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 \begin{figure}[t!] 39 \centerline{\includegraphics[angle=-90,width=4.in]{figures/waveu}} 40 \caption{Displacement at Source Point} 41 \label{WAVE FIG 1.2} 42 \end{figure} 43 44 \begin{figure}[t!] 45 \centerline{\includegraphics[angle=-90,width=4.in]{figures/wavea}} 46 \caption{Acceleration at Source Point} 47 \label{WAVE FIG 1.1} 48 \end{figure} 49 50 Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$-direction 51 which raise from zero to a maximum displacement $U\hackscore 0$ with in $\alpha$ seconds and then falls back to zero over the same period. In mathematical terms we use 52 \begin{eqnarray} \label{WAVE source} 53 u\hackscore 0(x\hackscore C,t)= U\hackscore 0 \frac{t^2}{\alpha^2} e^{1-\frac{t^2}{\alpha^2}} \ 54 \end{eqnarray} 55 for all $t\ge0$. In the simulations we will choose $\alpha=0.3$, see Figure~\ref{WAVE FIG 1.2} 56 and will apply the source as a constraint in a in a sphere of small radius around the point 57 $x\hackscore C$. 58 59 We use an explicit time integration scheme to calculate the displacement field $u$ at 60 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$ 61 which is defined by 62 \begin{eqnarray} \label{WAVE dyn 2} 63 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 64 \end{eqnarray} 65 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 66 \begin{eqnarray} \label{WAVE dyn 2b} 67 u\hackscore{,tt}=G(u) 68 \end{eqnarray} 69 where one sets $a^{(n)}=G(u^{(n-1)})$. 70 71 In our case $a^{(n)}$ is given by 72 \begin{eqnarray}\label{WAVE dyn 3} 73 \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j} 74 \end{eqnarray} 75 and boundary conditions 76 \begin{eqnarray} \label{WAVE natural at n} 77 \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0 78 \end{eqnarray} 79 derived from \eqn{WAVE natural} where 80 \begin{eqnarray} \label{WAVE dyn 3a} 81 \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}). 82 \end{eqnarray} 83 We also need to apply the constraint 84 \begin{eqnarray} \label{WAVE dyn 3aa} 85 a^{(n)}\hackscore{0}(x\hackscore C,t)= U\hackscore{0} 86 \frac{2}{\alpha^2} \left( 1- 5 \frac{t^2}{\alpha^2} +2 \frac{t^4}{\alpha^4} 87 \right) e^{1-\frac{t^2}{\alpha^2}} 88 \end{eqnarray} 89 which is derived from equation~\ref{WAVE source} by calculating the second order time derivative, 90 see~\ref{WAVE FIG 1.2}. Now we have converted our problem for displacement, $u^{(n)}$, into a problem for 91 acceleration, $a^{(n)}$, which now depends 92 on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. 93 94 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will 95 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 96 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is 97 \begin{equation}\label{WAVE dyn 100} 98 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . 99 \end{equation} 100 The natural boundary condition 101 \begin{equation}\label{WAVE dyn 101} 102 n\hackscore{j}X\hackscore{ij} =0 103 \end{equation} 104 is used. 105 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 106 \begin{equation}\label{WAVE dyn 6} 107 \begin{array}{ll} 108 D\hackscore{ij}=\rho \delta\hackscore{ij}& 109 X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ 110 \end{array} 111 \end{equation} 112 Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. We will apply 113 the constraint on a small sphere of radius $R$ around $x\hackscore C$ (we will use 3p.c. of the width of the domain): 114 \begin{equation}\label{WAVE dyn 6 1} 115 r\hackscore{i}(x) = 116 \left\{ 117 \begin{array}{rc} 118 1 & \|x-x_c\|\le R \\ 119 0 & \mbox{otherwise} 120 \end{array} 121 \right. 122 \end{equation} 123 The following script defines a the function \function{wavePropagation} which 124 implements the Verlet scheme to solve our wave propagation problem. 125 The argument \var{domain} which is a \Domain class object 126 defines the domain of the problem. \var{h} and \var{tend} are the time step size 127 and the end time of the simulation. \var{lam}, \var{mu} and 128 \var{rho} are material properties. 129 \begin{python} 130 def wavePropagation(domain,h,tend,lam,mu,rho,U0): 131 x=domain.getX() 132 # ... open new PDE ... 133 mypde=LinearPDE(domain) 134 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 135 kronecker=identity(mypde.getDim()) 136 137 # spherical source at middle of bottom face 138 139 xc=[width/2.,width/2.,0.] 140 # define small radius around point xc 141 # Lsup(x) returns the maximum value of the argument x 142 src_radius = 0.03*width 143 print "src_radius = ",src_radius 144 145 dunit=numpy.array([1.,0.,0.]) # defines direction of point source 146 147 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit) 148 # ... set initial values .... 149 n=0 150 # initial value of displacement at point source is constant (U0=0.01) 151 # for first two time steps 152 u=Vector(0.,Solution(domain)) 153 u_last=Vector(0.,Solution(domain)) 154 t=0 155 156 # define the location of the point source 157 L=Locator(domain,xc) 158 # find potential at point source 159 u_pc=L.getValue(u) 160 print "u at point charge=",u_pc 161 # open file to save displacement at point source 162 u_pc_data=FileWriter('./data/U_pc.out') 163 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2])) 164 165 while t

Properties

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