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

Revision 3379 - (show annotations)
Wed Nov 24 04:48:49 2010 UTC (8 years, 2 months ago) by gross
File MIME type: application/x-tex
File size: 17170 byte(s)
some clarification on lumping

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

## Properties

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