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

Revision 2580 - (show annotations)
Tue Aug 4 08:24:12 2009 UTC (9 years, 10 months ago) by gross
File MIME type: application/x-tex
File size: 18547 byte(s)
another example for the usage of matplot lib added to users guide.

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

Properties

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