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

Revision 6686 - (show 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 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2018 by The University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Apache License, version 2.0 8 9 % 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 \section{Wave Propagation} 17 \label{WAVE CHAP} 18 19 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 \begin{eqnarray}\label{WAVE general problem} 21 \rho u_{i,tt} - \sigma_{ij,j}=0 22 \end{eqnarray} 23 in a three dimensional block of length $L$ in $x_{0}$ 24 and $x_{1}$ direction and height $H$ in $x_{2}$ direction. 25 $\rho$ is the known density which may be a function of its location. 26 $\sigma_{ij}$ is the stress field\index{stress} which in case of an 27 isotropic, linear elastic material is given by 28 \begin{eqnarray}\label{WAVE stress} 29 \sigma_{ij} = \lambda u_{k,k} \delta_{ij} + \mu (u_{i,j} + u_{j,i}) 30 \end{eqnarray} 31 where $\lambda$ and $\mu$ are the Lam\'e coefficients\index{Lam\'e coefficients} 32 and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. 33 On the boundary the normal stress is given by 34 \begin{eqnarray}\label{WAVE natural} 35 \sigma_{ij}n_{j}=0 36 \end{eqnarray} 37 for all time $t>0$. 38 39 \begin{figure}[t!] 40 \centerline{\includegraphics[width=14cm]{waveu}} 41 \caption{Input Displacement at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).} 42 \label{WAVE FIG 1.2} 43 \end{figure} 44 45 \begin{figure} 46 \centerline{\includegraphics[width=14cm]{wavea}} 47 \caption{Input Acceleration at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).} 48 \label{WAVE FIG 1.1} 49 \end{figure} 50 51 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 In mathematical terms we use 55 \begin{eqnarray} \label{WAVE source} 56 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 \end{eqnarray} 58 for all $t\ge0$ where $\alpha$ is the width of the pulse and $t_{0}$ 59 is the time when the pulse changes from negative to positive. 60 In the simulations we will choose $\alpha=0.3$ and $t_{0}=2$ (see 61 \fig{WAVE FIG 1.2}) and apply the source as a constraint in a sphere of small 62 radius around the point $x_C$. 63 64 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 which is defined by 70 \begin{eqnarray} \label{WAVE dyn 2} 71 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 72 \end{eqnarray} 73 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 74 \begin{eqnarray} \label{WAVE dyn 2b} 75 u_{,tt}=G(u) 76 \end{eqnarray} 77 where one sets $a^{(n)}=G(u^{(n-1)})$. 78 79 In our case $a^{(n)}$ is given by 80 \begin{eqnarray}\label{WAVE dyn 3} 81 \rho a^{(n)}_{i}=\sigma^{(n-1)}_{ij,j} 82 \end{eqnarray} 83 and boundary conditions 84 \begin{eqnarray} \label{WAVE natural at n} 85 \sigma^{(n-1)}_{ij}n_{j}=0 86 \end{eqnarray} 87 derived from \eqn{WAVE natural} where 88 \begin{eqnarray} \label{WAVE dyn 3a} 89 \sigma_{ij}^{(n-1)} & = & \lambda u^{(n-1)}_{k,k} \delta_{ij} + \mu ( u^{(n-1)}_{i,j} + u^{(n-1)}_{j,i}). 90 \end{eqnarray} 91 We also need to apply the constraint 92 \begin{eqnarray} \label{WAVE dyn 3aa} 93 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 \end{eqnarray} 96 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 102 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 \begin{equation}\label{WAVE dyn 100} 107 D_{ij} a^{(n)}_{j} = - X_{ij,j}\; . 108 \end{equation} 109 The natural boundary condition 110 \begin{equation}\label{WAVE dyn 101} 111 n_{j}X_{ij} =0 112 \end{equation} 113 is used. 114 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 115 \begin{equation}\label{WAVE dyn 6} 116 \begin{array}{ll} 117 D_{ij}=\rho \delta_{ij}& 118 X_{ij}=-\sigma^{(n-1)}_{ij}\\ 119 \end{array} 120 \end{equation} 121 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 $x_C$ (we will use 3\% of the width of the domain): 124 \begin{equation}\label{WAVE dyn 6 1} 125 q_{i}(x) = 126 \left\{ 127 \begin{array}{l l} 128 1 & \quad \text{where $\|x-x_c\|\le R$}\\ 129 0 & \quad \text{otherwise.} 130 \end{array} 131 \right. 132 \end{equation} 133 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 \begin{python} 140 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 mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING) 149 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) 165 u_pc1.append(u_pc) 166 u_pc2.append(u_pc) 167 168 while t

## Properties

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