3D Wave Propagation

We want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation :

$$\rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0$$ (30)

in a three dimensional block of length $L$ in $x\hackscore{0}$ and $x\hackscore{1}$ direction and height $H$ in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location. $\sigma\hackscore{ij}$ is the stress field which in case of an isotropic, linear elastic material is given by

$$\sigma\hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})$$ (31)

where $\lambda$ and $\mu$ are the Lame coefficients and $\delta\hackscore{ij}$ denotes the Kronecker symbol. On the boundary the normal stress is given by

$$\sigma\hackscore{ij}n\hackscore{j}=0$$ (32)

for all time $t>0$.

At initial time $t=0$ the displacement $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given:

$$u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \m... ... at point charge, $x\hackscore{C}$} 0 & \mbox{elsewhere}\end{array} \right. u\hackscore{i,t}(0,x)=0 \;$$ and (33)

for all $x$ in the domain.

Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point $x\hackscore C$.

We use an explicit time-integration scheme to calculate the displacement field $u$ at 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 with constant time step size $h$ which is defined by

$$u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)}$$ (34)

for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form

$$u\hackscore{,tt}=G(u)$$ (35)

where one sets $a^{(n)}=G(u^{(n-1)})$.

In our case $a^{(n)}$ is given by

$$\rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}$$ (36)

and boundary conditions

$$\sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0$$ (37)

derived from Equation (2.33) where

$$\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}).$$ (38)

Now we have converted our problem for displacement, $u^{(n)}$, into a problem for acceleration, $a^(n)$, which now depends on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$.

In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will use the LinearPDE class of the esys.escript.linearPDEs to do so. The general form of the PDE defined through the LinearPDE class is discussed in Section 4.1. The form which is relevant here are

$$D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .$$ (39)

Implicitly, the natural boundary condition

$$n\hackscore{j}X\hackscore{ij} =0$$ (40)

is assumed.

With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$:

$$\begin{array}{ll} D\hackscore{ij}=\rho \delta\hackscore{ij}& X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \end{array}$$ (41)

The following script defines a the function wavePropagation which implements the Verlet scheme to solve our wave propagation problem. The argument domain which is a Domain class object defines the domain of the problem. h and tend are the step size and the time mark after which the simulation will be terminated respectively. lam, mu and rho are material properties.

Notice that all coefficients of the PDE which are independent from time $t$ are set outside the while loop. This allows the LinearPDE class to reuse information as much as possible when iterating over time.

We have seen most of the functions in previous examples but there are some new functions here: grad(u) returns the gradient $u\hackscore{i,j}$ of $u$ (in fact grad(g)[i,j]= $=u\hackscore{i,j}$). It is pointed out that in general there are restrictions on the argument of the grad function, for instance for esys.finley the statement grad(grad(u)) will raise an exception. trace(g) returns the sum of the main diagonal elements g[k,k] of g and transpose(g) returns the transposed of the matrix g (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$).

We initialize the values of u and u_last to be $U\hackscore 0$ for a small sphere of radius, src_radius around the point source located at $x\hackscore C$. These satisfy the initial conditions defined in Equation (2.34).

The output U_pc.out and vtu files are saved in a directory called data. The U_pc.out stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of the displacement vector $u$ at the point source. These can be plotted easily using any plotting package. In gnuplot the command: plot 'U_pc.out' u 1:2 title 'U_x' w l lw 2, 'U_pc.out' u 1:3 title 'U_y' w l lw 2, 'U_pc.out' u 1:4 title 'U_z' w l lw 2 will reproduce Figure 2.6.

The statement myPDE.setLumpingOn() switches on the usage of an aggressive approximation of the PDE operator, in this case formed by the coefficient D (actually the discrete version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of an additional error, very fast solving of the PDE when the solution is requested. In the general case, the presence of A, B or C setLumpingOn will produce wrong results.

One of the big advantages of the Verlet scheme is the fact that the problem to be solved in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping). The problem becomes so simple because we use the stress from the last time step rather then the stress which is actually present at the current time step. Schemes using this approach are called an explicit time integration scheme . The backward Euler scheme we have used in Chapter 2.2 is an example of an implicit schemes . In this case one uses the actual status of each variable at a particular time rather then values from previous time steps. This will lead to a problem which is more expensive to solve, in particular for non-linear problems. Although explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time steps then implicit schemes and can suffer from stability problems. Therefore they need a very careful selection of the time step size $h$.

An easy, heuristic way of choosing an appropriate time step size is the Courant condition which says that within a time step a information should not travel further than a cell used in the discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from

$$h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x$$ (42)

where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of the formula.

The following script uses the wavePropagation function to solve the wave equation for a point source located at the bottom face of a block. The width of the block in each direction on the bottom face is $10$km ( $x\hackscore 0$ and $x\hackscore 1$ directions ie. l0 and l1). ne gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions. The depth of the block is aligned with the $x\hackscore{2}$-direction. The depth (l2) of the block in the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the of the depth is chosen such that the elements become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the computation would be faster. Brick( $n\hackscore 0, n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$) is an esys.finley function which creates a rectangular mesh with $n\hackscore 0 \times n\hackscore 1 \times n\hackscore2$ elements over the brick $[0,l\hackscore 0] \times [0,l\hackscore 1] \times [0,l\hackscore 2]$.

The script is available as wave.py in the example directory . Before running the code make sure you have created a directory called data in the current working directory. To visualize the results from the data directory:

 mayavi -d usoln.1.vtu -m SurfaceMap &

You can rotate this figure by clicking on it with the mouse and moving it around. Again use Configure Data to move backwards and forward in time, and also to choose the results (acceleration, displacement or $u\hackscore x$) by using Select Scalar. Figure 2.7 shows the results for the displacement at various time steps.

Figure 2.6: Amplitude at Point source

t!

Figure: Selected time steps ( $n = 0,1,30,100,300$ and 2880) of a wave propagation over a $10$km$\times 10$km$\times 3.125$km block from a point source initially at $(5$km$,5$km$,0)$ with time step size $h=0.02083$. Color represents the displacement. Here the view is oriented onto the bottom face.

t