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

Revision 6686 - (hide annotations)
Thu Jun 14 04:01:43 2018 UTC (2 years, 4 months ago) by aellery
File MIME type: application/x-tex
File size: 17768 byte(s)
Fixed some minor typos

 1 caltinay 5293 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 jfenwick 6651 % Copyright (c) 2003-2018 by The University of Queensland 4 caltinay 5293 5 % 6 % Primary Business: Queensland, Australia 7 jfenwick 6112 % Licensed under the Apache License, version 2.0 8 9 caltinay 5293 % 10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC) 11 % Development 2012-2013 by School of Earth Sciences 12 % Development from 2014 by Centre for Geoscience Computing (GeoComp) 13 % 14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 15 16 gross 4537 \section{Wave Propagation} 17 \label{WAVE CHAP} 18 ksteube 1811 19 jfenwick 3295 In this next example we want to calculate the displacement field $u_{i}$ for any time $t>0$ by solving the wave equation:\index{wave equation} 20 jgs 107 \begin{eqnarray}\label{WAVE general problem} 21 jfenwick 3295 \rho u_{i,tt} - \sigma_{ij,j}=0 22 jgs 107 \end{eqnarray} 23 jfenwick 3295 in a three dimensional block of length $L$ in $x_{0}$ 24 and $x_{1}$ direction and height $H$ in $x_{2}$ direction. 25 caltinay 3280 $\rho$ is the known density which may be a function of its location. 26 jfenwick 3295 $\sigma_{ij}$ is the stress field\index{stress} which in case of an 27 caltinay 3280 isotropic, linear elastic material is given by 28 \begin{eqnarray}\label{WAVE stress} 29 jfenwick 3295 \sigma_{ij} = \lambda u_{k,k} \delta_{ij} + \mu (u_{i,j} + u_{j,i}) 30 jgs 107 \end{eqnarray} 31 sshaw 5284 where $\lambda$ and $\mu$ are the Lam\'e coefficients\index{Lam\'e coefficients} 32 jfenwick 3295 and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. 33 jgs 107 On the boundary the normal stress is given by 34 caltinay 3280 \begin{eqnarray}\label{WAVE natural} 35 jfenwick 3295 \sigma_{ij}n_{j}=0 36 jgs 107 \end{eqnarray} 37 gross 2502 for all time $t>0$. 38 lkettle 581 39 gross 2502 \begin{figure}[t!] 40 caltinay 3280 \centerline{\includegraphics[width=14cm]{waveu}} 41 jfenwick 3295 \caption{Input Displacement at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).} 42 gross 2502 \label{WAVE FIG 1.2} 43 \end{figure} 44 45 caltinay 3280 \begin{figure} 46 \centerline{\includegraphics[width=14cm]{wavea}} 47 jfenwick 3295 \caption{Input Acceleration at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).} 48 gross 2502 \label{WAVE FIG 1.1} 49 \end{figure} 50 51 jfenwick 3295 Here we are modelling a point source at the point $x_C$ in the 52 $x_{0}$-direction which is a negative pulse of amplitude 53 $U_{0}$ followed by the same positive pulse. 54 caltinay 3280 In mathematical terms we use 55 gross 2502 \begin{eqnarray} \label{WAVE source} 56 jfenwick 3295 u_{0}(x_C,t)= U_{0} \; \sqrt{2} \frac{(t-t_{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t_{0})^2}{\alpha^2}} \ 57 jgs 107 \end{eqnarray} 58 jfenwick 3295 for all $t\ge0$ where $\alpha$ is the width of the pulse and $t_{0}$ 59 caltinay 3280 is the time when the pulse changes from negative to positive. 60 jfenwick 3295 In the simulations we will choose $\alpha=0.3$ and $t_{0}=2$ (see 61 caltinay 3280 \fig{WAVE FIG 1.2}) and apply the source as a constraint in a sphere of small 62 jfenwick 3295 radius around the point $x_C$. 63 jgs 107 64 caltinay 3280 We use an explicit time integration scheme to calculate the displacement field 65 $u$ at certain time marks $t^{(n)}$, where $t^{(n)}=t^{(n-1)}+h$ with time 66 step size $h>0$. 67 In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. 68 We use the Verlet scheme\index{Verlet scheme} with constant time step size $h$ 69 jgs 107 which is defined by 70 \begin{eqnarray} \label{WAVE dyn 2} 71 lkettle 575 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 72 jgs 107 \end{eqnarray} 73 lkettle 581 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 74 jgs 107 \begin{eqnarray} \label{WAVE dyn 2b} 75 jfenwick 3295 u_{,tt}=G(u) 76 jgs 107 \end{eqnarray} 77 gross 660 where one sets $a^{(n)}=G(u^{(n-1)})$. 78 jgs 107 79 In our case $a^{(n)}$ is given by 80 \begin{eqnarray}\label{WAVE dyn 3} 81 jfenwick 3295 \rho a^{(n)}_{i}=\sigma^{(n-1)}_{ij,j} 82 jgs 107 \end{eqnarray} 83 and boundary conditions 84 gross 593 \begin{eqnarray} \label{WAVE natural at n} 85 jfenwick 3295 \sigma^{(n-1)}_{ij}n_{j}=0 86 jgs 107 \end{eqnarray} 87 derived from \eqn{WAVE natural} where 88 \begin{eqnarray} \label{WAVE dyn 3a} 89 jfenwick 3295 \sigma_{ij}^{(n-1)} & = & \lambda u^{(n-1)}_{k,k} \delta_{ij} + \mu ( u^{(n-1)}_{i,j} + u^{(n-1)}_{j,i}). 90 jgs 107 \end{eqnarray} 91 gross 2502 We also need to apply the constraint 92 \begin{eqnarray} \label{WAVE dyn 3aa} 93 jfenwick 3295 a^{(n)}_{0}(x_C,t)= U_{0} 94 \frac{\sqrt(2.)}{\alpha^2} (4\frac{(t-t_{0})^3}{\alpha^3}-6\frac{t-t_{0}}{\alpha})e^{\frac{1}{2}-\frac{(t-t_{0})^2}{\alpha^2}} 95 gross 2502 \end{eqnarray} 96 caltinay 3280 which is derived from equation~\ref{WAVE source} by calculating the second 97 order time derivative (see \fig{WAVE FIG 1.1}). 98 Now we have converted our problem for displacement, $u^{(n)}$, into a problem 99 for acceleration, $a^{(n)}$, which depends on the solution at the previous two 100 time steps $u^{(n-1)}$ and $u^{(n-2)}$. 101 jgs 107 102 caltinay 3280 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, 103 and we will use the \LinearPDE class of the \linearPDEs package to do so. 104 The general form of the PDE defined through the \LinearPDE class is discussed 105 in \Sec{SEC LinearPDE}. The form which is relevant here is 106 jgs 107 \begin{equation}\label{WAVE dyn 100} 107 jfenwick 3295 D_{ij} a^{(n)}_{j} = - X_{ij,j}\; . 108 jgs 107 \end{equation} 109 ksteube 1316 The natural boundary condition 110 jgs 107 \begin{equation}\label{WAVE dyn 101} 111 jfenwick 3295 n_{j}X_{ij} =0 112 jgs 107 \end{equation} 113 ksteube 1316 is used. 114 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 115 caltinay 3280 \begin{equation}\label{WAVE dyn 6} 116 jgs 107 \begin{array}{ll} 117 jfenwick 3295 D_{ij}=\rho \delta_{ij}& 118 X_{ij}=-\sigma^{(n-1)}_{ij}\\ 119 jgs 107 \end{array} 120 \end{equation} 121 caltinay 3280 Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. 122 We will apply the constraint on a small sphere of radius $R$ around 123 jfenwick 3295 $x_C$ (we will use 3\% of the width of the domain): 124 caltinay 3280 \begin{equation}\label{WAVE dyn 6 1} 125 jfenwick 3295 q_{i}(x) = 126 gross 2502 \left\{ 127 caltinay 3280 \begin{array}{l l} 128 1 & \quad \text{where $\|x-x_c\|\le R$}\\ 129 caltinay 5296 0 & \quad \text{otherwise.} 130 gross 2502 \end{array} 131 \right. 132 \end{equation} 133 caltinay 3280 The following script defines the function \function{wavePropagation} which 134 implements the Verlet scheme to solve our wave propagation problem. 135 The argument \var{domain} which is a \Domain class object defines the domain 136 of the problem. 137 \var{h} and \var{tend} are the time step size and the end time of the simulation. 138 \var{lam}, \var{mu} and \var{rho} are material properties. 139 jgs 107 \begin{python} 140 caltinay 3280 def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0): 141 # lists to collect displacement at point source which is returned 142 # to the caller 143 ts, u_pc0, u_pc1, u_pc2 = [], [], [], [] 144 145 x=domain.getX() 146 # ... open new PDE ... 147 mypde=LinearPDE(domain) 148 sshaw 4821 mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING) 149 caltinay 3280 kronecker=identity(mypde.getDim()) 150 dunit=numpy.array([1., 0., 0.]) # defines direction of point source 151 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit) 152 # ... set initial values .... 153 n=0 154 # for first two time steps 155 u=Vector(0., Solution(domain)) 156 u_last=Vector(0., Solution(domain)) 157 t=0 158 # define the location of the point source 159 L=Locator(domain, xc) 160 # find potential at point source 161 u_pc=L.getValue(u) 162 print("u at point charge=",u_pc) 163 ts.append(t) 164 u_pc0.append(u_pc[0]) 165 u_pc1.append(u_pc[1]) 166 u_pc2.append(u_pc[2]) 167 gross 2580 168 caltinay 3280 while t

## Properties

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