A numerical analysis of the seismic wave equation in different layers by the finite element method using fenics

A numerical analysis of the seismic wave equation in different layers by the finite

element method using fenics

by

DANIEL JAMES TARPLETT

THESIS for the degree of

MASTER OF SCIENCE

(Master i anvendt matematikk og mekanikk)

Faculty of Mathematics and Natural Sciences University of Oslo

December 2014

Det matematisk- naturvitenskapelige fakultet Universitetet i Oslo

Abstract

In this thesis we will investigate the seismic wave equation in different layers by using the finite element method in space and the finite difference method in time. The performance of the programming will be done by comparisons with analytical solutions by using test-solution methods, and convergence tests will be used for error control.

Acknowledgements

Firstly, I would like to thank Geir K. Pedersen and Mikael Mortensen for all their advice and support in the process of writing, and for the help in both the mathematics and programming. I would like to thank Miroslav Kuchta for all the help he has given in the FEniCS Q&A forum and in person. I would also like to thank Finn Løvholt and Valerie Maupin for their role as supervisors. I would like to thank my fellow students for their friendship, and the faculty staff for all the help they have given. Lastly, I would like to thank my friends and family for their continuous support during the process of writing.

Contents

1 Introduction 7

2 Theory 8

2.1 Governing equations . . . 8

2.2 The finite difference method . . . 9

2.3 The finite element method . . . 9

2.4 Discretizing the wave equation . . . 10

2.5 Discretizing the momentum equation . . . 11

2.6 Boundary conditions . . . 12

2.7 Sponge layers . . . 13

2.8 Error control, stability and convergence . . . 13

3 Waves on a sponge layer 14 3.1 An analytic solution . . . 15

3.2 Simulations and results . . . 16

3.3 Conclusion . . . 17

4 The Seismic Wave Equation with Test Solutions 20 4.1 P and S wave analytic solutions . . . 22

4.2 Simulations and results . . . 22

4.3 Conclusion . . . 23

5 Seismic test solutions with a given stress 26 5.1 P and S-wave analytic solutions . . . 27

5.2 Simulations and results . . . 28

5.3 Conclusion . . . 28

6 A Two layer model with vertical incidence 31 6.1 P-wave analytic solutions . . . 32

6.2 S-wave analytic solutions . . . 35

6.3 Simulations and results . . . 36

6.4 Conclusion . . . 37

7 A two layer model with an oblique angle 42 7.1 An Analytic solution with an incoming P-wave . . . 43

7.2 An analytic solution from an incoming S-wave . . . 45

8 Discussion 46 9 Appendix 48 9.1 Code for the sponge layer project . . . 48

9.2 Code for the seismic test solution with dirichlet conditions . . . . 51

9.3 Code for the seismic test solutions with given surface stress . . . 53

9.4 Code for the seismic waves on multiple layers . . . 57

List of Figures

1 The problem where waves travel with horizontal incidence into a sponge layer 15 2 Figure of the errors in the fluid domain for the run withL= 2,x_s= 1 and

a linear damping in the sponge layer. (a) shows the errors for the coarse mesh, (b) shows the errors for the finer mesh, and (c) shows the errors for the finest mesh . . . 18 3 Figure of the errors in the fluid domain for L = 3, x_s = 1 and a linear

damping in the sponge layer. (a) shows the errors for the coarse mesh, (b) shows the errors for the finer mesh, and (c) shows the errors for the finest mesh 19 4 Figure of the errors in the fluid domain withL= 3,x_s= 1 and a quadratic

damping function in the sponge layer. (a) shows the errors for the coarse mesh, (b) shows the errors for the finer mesh, and (c) shows the errors for the finest mesh . . . 20 5 The rectanguar domain used in the problem . . . 21 6 Errors for the x and z-components of displacement for a P-wave with an

angle of 71.57⁰with the x-axis. (a) and (b) show the x and z-displacements for a 24x24 mesh respectively, and a time step of 0.0075. figures (c) and (d) show the x and z-displacements for a 96x96 mesh respectivley, and a time step of 0.001875.. . . 24 7 Errors for the x and z-components of displacement for an S-wave with an

angle of 71.57⁰with the x-axis. (a) and (b) show the x and z-displacements for a 24x24 mesh respectively, and a time step of 0.0075. figures (c) and (d) show the x and z-displacements for a 96x96 mesh respectivley, and a time step of 0.001875.. . . 25 8 The problem with test solutions for dirichlet boundary conditions and a given

surface stress . . . 26 9 Figures of the displacement errors for a P-wave propagating with an angle

of θ = 71.57⁰ with respect to the x-axis. (a) and (b) show the x and z- displacements for a 24x24 mesh with a time step of 0.0075. (c) and (d) show the x and z-displacement errors for a 96x96 mesh with time step 0.0001875 . 29 10 Figures of the displacement errors for an S-wave propagating with an angle

of θ = 71.57⁰ with respect to the x-axis. (a) and (b) show the x and z- displacements for a 24x24 mesh with a time step of 0.0075. (c) and (d) show the x and z-displacement errors for a 96x96 mesh with time step 0.0001875 . 30 11 A two layer model for waves traveling at vertical incidence with the boundaries 32 12 A two layer model for P-waves traveling at vertical incidence with an internal

boundary and a free surface . . . 32 13 A two layer model for S-waves traveling at vertical incidence with an internal

boundary and a free surface . . . 35 14 Errors in the x and z components for P-waves hitting a solid-solid boundary.

Figure (a) and (b) shows the x and z-component errors for a 12x24 mesh respectively. Figures (c) and (d) shows the x and z-component errors for a 48x96 mesh respectively . . . 38

15 Errors in the x and z components for P-waves hitting a solid-liquid boundary.

Figure (a) and (b) shows the x and y-component errors for a 12x24 mesh respectively. Figures (c) and (d) shows the x and z-component errors for a 48x96 mesh respectively . . . 39 16 Errors in the x and z components for S-waves hitting a solid-solid boundary.

Figure (a) and (b) shows the x and z-component errors for a 12x24 mesh respectively. Figures (c) and (d) shows the x and z-component errors for a 48x96 mesh respectively . . . 40 17 Errors in the x and z components for S-waves hitting a solid-liquid boundary.

Figure (a) and (b) shows the x and z-component errors for a 12x24 mesh respectively. Figures (c) and (d) shows the x and z-component errors for a 48x96 mesh respectively . . . 41 18 The two layer domain for waves sent with an oblique angle . . . 43 19 The problem for a P-wave hitting the boundary between solid and fluid . . 44 20 The problem for an S-wave hitting the boundary between solid and fluid . . 45 21 The earthquake model for future study. The model includes the ocean, crust

and continent, where the earthquake has it´s source between the crust and continent. . . . 47

List of Tables

1 Table of numerical results for 3 different simulations. Lis the total length of the domain, x_s is the x coordinate of the boundary between fluid and sponge. b_landb_q denotes the linear and quadratic damping functions used.

∆x and ∆z are the element spacings in the x and z-directions, and ∆tis the time step. E_max and E_l2n are the maximum and L2 norm errors in the simulations, and C_max and C_l2n are the error reduction rates for the maximum and L2 norm errors with the respect to the previous simulation . 17 2 Table of the calculated amplitudes of the reflected waves from the sponge

layer in all 3 simulations. L denotes the length of the domain,x_sthe coor- dinate of the boundary between fluid and sponge, andb_landb_q the linear and quadratic damping functions respectivley.. . . 21 3 Table containing the numerical results of the simulations of the seismic wave

equation with a P wave test solution. The angleθgives the angle of prop- agation with the x-axis, ∆xand ∆zgive the element spacings in the x and y direction. ∆tis the time step. E_maxandE_L2denotes the maximum and L2 norm errors respectvely. C_maxandC_L2are the error reduction rates for the maximum and L2 norm errors with respect to the previous simulation.

A_rare the estimated amplitudes from the reflected waves . . . 23

4 Table containing the numerical results of the simulations of the seismic wave equation with an S wave test solution. The angleθgives the angle of prop- agation with the x-axis, ∆xand ∆zgive the element spacings in the x and z-direction. ∆tis the time step. E_maxandE_L2denotes the maximum and L2 norm errors respectivley. C_maxandC_L2are the error reduction rates for the maximum and L2 norm errors with respect to the previous simulation.

A_rare the estimated amplitudes of the reflected waves . . . 26 5 Table containing the numerical results of the simulations of the seismic wave

equation with P-wave test solutions. The angle θgives the angle of prop- agation with respect to the x-axis, ∆xand ∆z give the element spacings in the x and z direction. ∆tis the time step. E_max and E_L2 denotes the maximum and L2 norm errors. C_maxandC_L2are the error reduction rates with respect to the previous simulation. A_rare the estimated amplitudes of the reflected waves. . . 28 6 Table containing the numerical results of the simulations of the seismic wave

equation with S-wave test solutions. The angle θ gives the angle of prop- agation with respect to the x-axis, ∆xand ∆z give the element spacings in the x and z direction. ∆tis the time step. E_max and E_L2 denotes the maximum and L2 norm errors. C_maxandC_L2are the error reduction rates with respect to the previous simulation. A_rare the estimated amplitudes of the reflected waves. . . 31 7 Results for P-waves vertically incident on a solid-solid boundary and a free

surface. . . . 37 8 Results for S-waves vertically incident on a solid-solid boundary and a free

surface. . . . 37 9 Results for P-waves vertically incident on a solid-liqid boundary and a free

surface. . . . 37 10 Results for S-waves vertically incident on a solid-liqid boundary and a free

surface. . . . 42

1 Introduction

Elastic waves in the earth are commonly described as seismic waves, and are produced by earthquakes, explosions and similar events. The study of these waves are important in their own right for warning and detection purposes, but the mathematical theory can also be used in other applications of science. It is common to use potential theory when studying seismic waves and seismology, but in this here we will concentrate on more direct solutions of the seismic wave equation. Numerical experiments will be done by using the finite difference method in time, and the finite element method in space. The finite element method is chosen because of it’s ability to handle natural boundary conditions, but also because of it’s ability to handle more complex geometries. The im- plementation is done in python using the FEniCS software, as it contains a scripting enviornment and syntax close to the mathematical formalism in the finite element method. In the numerical testing, we will also introduce a con-

cept called test-solutions for simplifying analytic solutions. The overall goal of the thesis is to examine how FEniCS handles an implementation of the seismic wave equation with one and two layers of material. The work is divided into four seperate projects examining the different aspects of the method, and each with their own separate conclusions. We have also included a fifth section, where the mathematics for a further problem is discussed.

2 Theory

In this thesis, we will work with 2D functions in the x-z plane with the y axis pointing inward. We will use dyadic notation where boldface characters indicate vector quantities.

2.1 Governing equations

The scalar wave equation with a variable wave velocity and a damping term can be expressed by:

∂²u

∂t² +b∂u

∂t =∇(c∇u) (1)

whereu=u(x, z, t) is the displacement,b(x, z) is the damping term, andc(x, z) is the variable wave velocity. Under the continuum assumption as explained by Kundu and Cohen [2008, see pp. 4-5] the momentum equation for small particle displacements can be found from the momentum equation, as done by Stein and Wysession [2009], and is given by:

ρ∂²u

∂t² =∇ ·σ+f (2)

where u =u(x, z, t) is the velocity, ρ = ρ(x, z) is the density, σ is the stress tensor andf=f(x, z, t) denotes the body forces. Equation (2) can also be called the navieres primitive equation of motion. By studying the strain of a material in 3 dimensions as done by Stein and Wysession [2009, pp. 49-51], we can find the stress tensor

σ=λ(∇ ·u)I+µ(∇u+∇u^T) (3) where we assume the material to be linear elastic, isotropic and that the stresses are symmetric. σ is the stress tensor, u is the displacement vector, I is the identity matrix,λis lame´es first constant, andµis the shear modulus. Inserting equation (3) into equation (2), we get

u_tt=(λ+µ)

ρ ∇(∇ ·u) +µ

ρ∇²u+f (4)

which is the seismic wave equation.

2.2 The finite difference method

The classic definitions for discretizing derivatives can be found in multiple text- books and multiple websites. Tveito and Winther [2005, pp. 46] gives a good derivation by using taylor series. We invoke the notation uⁿ = u(x, y, z, t), uⁿ⁻¹ =u(x, y, z, t−∆t) and uⁿ⁺¹ = u(x, y, z, t+ ∆t). We approximate first derivatives by using the midpoint rule:

ut≈uⁿ⁺¹−uⁿ⁻¹

2∆t +O(∆t²) (5)

and second derivatives by the central difference formula:

utt≈uⁿ⁺¹−2uⁿ+uⁿ⁻¹

∆t² +O(∆t²) (6)

where we notice that both approximations have an second order error in time.

2.3 The finite element method

The finite element method is a vast collection of mathematical principals and ideas put together in a comprehensive framework for solving differential equa- tions and boundary value problems. The full detail of the method is beyond the scope of this thesis, but we review the basic idea as given by Anders Logg [2012, pp. 77-94]. We divide the domain into triangles for two dimensional domains, and tetrahedrons for three dimensional domains and call these subdomains for elements. We then seek polynomial approximations to the unknown in each element and then assemble all the parts together to find the global system. We assume that our function can be approximated by the sum:

u(x) =

N

X

j=0

cjψj(x) (7)

wherecj are unknown constants, xdenotes the spatial coordinates andψ_j are given functions of an arbitrary degree. The functionsψ_j are commonly refered to as basis functions or weight functions. Suppose our problem is to approximate our solution u with a function f. This gives the simple solution:

u(x, y)≈f(x, y) (8)

And the difference between these two give a residual:

R(x, y) =f(x, y)−u(x, y) (9) The point is now to minimize this residual as much as possible, and this can be done by methods including the interpolation, least squares or weighted residuals method as explained by Langtangen [1999, see pp. 142-144]. We will focus on

the latter method, as this is used by the FEniCS software. We define a function space that is spanned by the basis functions:

Vˆ = span{ψj} And seek weight functions:

v∈Vˆ

such that the inner product of the residual and the test function is zero:

Z

Ω

R(x, y)vdΩ = 0 ∀v∈Vˆ (10)

Inserting the expression for R from equation (9) into the inner product in equa- tion (10) we get the equation:

Z

Ω

uvdΩ = Z

Ω

f vdΩ (11)

Equation (11) is the variational form of the problem, and constitutes a linear system of equations. The point of the finite element method is to solve this system using one of many integration methods, including LU solvers and krylov solvers. We end the review of the finite element method here, and interested readers can read the fenics book Anders Logg [2012] or many other good publi- cations on the topic. The rest of this thesis will focus on the variational forms while FEniCS handles the rest.

2.4 Discretizing the wave equation

We first apply the finite difference scheme for time using equations (5) and (6) for the time derivatives in equation (1) and get the explicit formula in time:

uⁿ⁺¹−2uⁿ+uⁿ⁻¹

∆t² +buⁿ⁺¹−uⁿ⁻¹

2∆t =∇(c∇uⁿ) +fⁿ (12) By further introducing the help functions:

A= 1

1 + ^b∆t₂ B=b∆t

2 −1 We get the explicit formula for the time stepping:

uⁿ⁺¹= 2Auⁿ+ABuⁿ⁻¹+A∇ ·(c∇uⁿ) +A∆t²fⁿ (13) The space variables are then solved by using the finite element method. Using the chain rule for the laplace term:

∇ ·(c∇uⁿv) =∇ ·(c∇uⁿ)v+c∇uⁿ∇v

and applying green’s theorem, as done by Tveito and Winther [2005, see]:

Z

Ω

∇ ·(c∇uⁿv)dΩ = Z

Γ

n·c∇uⁿvdΩ

The variational form of equations (13) is:

Z

Ω

uⁿ⁺¹vdΩ = 2 Z

Ω

AuⁿvdΩ + Z

Ω

ABuⁿ⁻¹vdΩ

− Z

Ω

cA∇uⁿ∇vdΩ + Z

Γ

An· ∇uⁿvdΓ + ∆t²

Z

Ω

AfⁿvdΩ

(14)

2.5 Discretizing the momentum equation

The momentum equation is vector valued, and has components in the x,y, and z directions. The weight functions must therefore also have components in the x,y,z direction. In our two dimentional description, we get the velocity vector

u=ui+wk (15)

In all the projects, we will work with the same nodes foruandv. we use local form functions NI where I is the global node number, and we use the local weight functionsw_I =NI. the vector weight function has the form:

w=axNIi+azNIk (16) where ax = 1 andaz = 0 gives the x-component of the variational form, and ax= 0 and az = 1 gives the z-component. Using the chain rule on the stress tensor as we did for the wave equation, we get

∇ ·(σ·w) = (∇ ·σ)·w+σ:∇w And applying green’s theorem

Z

Ω

∇ ·(σ·w)dΩ = Z

Γ

n·σ·wdΓ

we get the variational form of equation (2) Z

Ω

ρuⁿ⁺¹·wdΩ = 2 Z

Ω

ρuⁿ·wdΩ− Z

Ω

ρuⁿ⁻¹·wdΩ

+ ∆t² Z

Γ

n·σⁿ·wdΓ−∆t² Z

Ω

σⁿ:∇wdΩ + ∆t²

Z

Ω

fⁿ·wdΩ

(17)

2.6 Boundary conditions

In this thesis, we will give 4 different boundary conditions that are valid for seismic waves and their interactions between solids, liquids and air.

Fixed boundary

At the fixed boundary, the velocity or displacement is known at the boundary node I.w_I is not used and the variational form in equation (17) is not solved.

Instead, A value is directly inserted into the node points at the boundary:

u=U(x, z, t) (18)

whereUis a given boundary function.

Free boundary

The free boundary condition gives a known stress at the boundary, making the boundary integral term in (17) solvable.

n·σ=σ_n (19)

Whereσis the stress tensor,nis the normal vector andσ_n is a given function for the stress at the boundary. σ_n is often set to zero to model free surface boundary conditions..

Internal solid-solid boundary

The solid solid boundary condition describes a type of interaction between two solid media, like the Moho discontinuity discussed by Stein and Wysession [2009, see pp. 122] at the crust-mantle boundary. In the solid-solid interface, all velocity component and tractions must be continuous.

σ⁽¹⁾ =σ⁽²⁾

u⁽¹⁾ =u⁽²⁾ (20)

where u⁽¹⁾ and u⁽²⁾ are the velocity vectors in layers 1 and 2, and σ⁽¹⁾ and σ⁽²⁾ are the shear stresses in layers 1 and 2. In the finite element method, the solid-solid boundary gives duplicate nodes at the boundary, and are assembled into the global system.

Internal solid-liquid boundary

The solid-liquid boundary condition describes the interactions between solid and liquid media, like the sea floor and ocean. Due to the vanishing shear stress, the normal tractions and displacements need to be continuous. The shear stress

in the solid vanishes at the boundary, and there is no restriction on the shear displacements.

σ⁽¹⁾_n =σ⁽²⁾_n σ⁽¹⁾_s = 0 u⁽¹⁾_n =u⁽²⁾_n u⁽¹⁾_s 6=u⁽²⁾_s

(21)

whereσ_ndenotes the normal stress,σ_sis the shear stress,u_nis the normal dis- placements, andu_sdenotes the shear displacements. The solid-liquid boundary produces duplicate nodes at the boundary as for the solid-solid boundary, and are assembled into the global system.

2.7 Sponge layers

In the finite element method, boundaries are forced on the domain. If no bound- ary is specified as a essential boundary condition, the natural boundary condi- tions are applied. This gives difficulties if one wants the solution to flow out of the domain. One solution to this is by using sponge layers. The sponge layer is a type of damping layer often used to curb solutions to rest. We present two types of sponge layers: The damping function and the input method. The damping function can be implemented by inserting:

d=b∂u

∂t (22)

into the differential equation. This causes natural damping where

b=b(x;α1, ..., alphaN) is the damping function. the valuesα1, ..., αN are con- stants that depend on the problem and domain. Large values of b cause a larger damp effect. The damping function is easily applied to simple geometries, but finding a functionb(x) in more complex boundaries can be difficult. In the input method we force the solution to be reduced by setting

u=µu (23)

for every time step in the domain considered. µ ∈ (0,1) gives the damping, where 0 is absolute damping and 1 is no damping effect. The input method is easily applied to more complex geometries, but the method itself can produce large discontinuities in the domain, giving total reflections instead of dampings.

2.8 Error control, stability and convergence

The combination of the finite difference and finite element method gives a ex- plicit set of equations to be solved at each time step, and by this method we also impose stability conditions on the numerical scheme. Although important, the mathematics is involved, and left for further analysis, yet we will keep in mind the existence of stability in our programming. Another important property of

the numerical scheme is the existense of numerical dispersion. For waves with an angular frequencyω, the numerical scheme produces a numerical frequency ˆ

ω where ω 6= ˆω. Such an analysis is also quite involved in the finite element method, and is also left for further study, yet Langtangen [1999, see pp. 656]

gives a nice review of the method for a finite difference scheme. In the numerical testing, we will have analytic solutions to compare our simulations with, and we put an emphasis on investigation of errors. The L2 norm error can be defined as

EL2= s

PN

i=0(ue−u)²

N (24)

where EL2 is the L2 norm error, ue is the exact solution, u is the numerical solution and N is the number of nodes. For P1 elements we get a second order error in the spatial coordinates. Combined with the second order errors in the finite difference schemes for the time discretization, we get the error in the scheme

E1=Ax(∆x)²+Az(∆z)²+At(∆t)² where we notice that halving this error gives

E2=Ax(∆x

2 )²+Az(∆z

2 )²+At(∆t 2 )² and that the ratio between the errors are

E2

E1

= 0.25

This shows that the error is reduced by a factor 4 when halving spatial and time steps. We will call the number 0.25 the error reduction rate. The spatial and time steps can be collected into a common parameterh, such that the error is given by

E=Ch² (25)

where E is the error, C is some constant andh=h(∆x,∆z,∆t) is a common parameter for the spatial and time steps. The exponent is commonly referred to as the convergence rate.

3 Waves on a sponge layer

In this first project, the performance of a sponge layer will be tested for a simple wave problem on a rectangular domain. Waves are sent into the sponge layer, and it’s ability to damp out the motion will be analyzed. We assume a rectangular domain Ω with lengthLand heightH. The domain is divided into two sub domains Ω1 and Ω2 divided by a vertical line at the pointx=xS, We give the first and second domain the lengths Lp and Ls respectively, and the height of both domains areH. the subscripts p and s are short for p-wave and sponge layer. The problem is shown in figure 1. Each domain is divided into np×mandns×melements respectively.

u=U(z, t)

u= 0

∂u

∂z = 0

∂u

∂x = 0

b= 0 b=b(x)

Fluid Sponge

Figure 1: The problem where waves travel with horizontal incidence into a sponge layer

3.1 An analytic solution

In the fluid layer we have no damping and a constant wave velocityc1. In the sponge layer we apply a damping coefficient only dependent on x and a constant wave velocityc2. Equation (1) then reduces to:

∂²u1

∂t² =c²₁∇²u x∈(0, Lp) (26)

∂²u2

∂t² +b(x)∂u2

∂t =c²₂∇u x∈(Lp, Ls) (27) For the fluid and sponge respectivly. u1 is the displacement in the fluid layer, andu2is the displacement in the sponge layer. The boundary value problem is subject to 4 boundary conditions in the domain. At the topy =H we assume no displacements. At the bottom y = 0 and at the right x = Ls we assume Neumann boundary conditions. At the left hand boundaryx= 0 we have an inflow condition. All four boundary conditions are stated as

u1(x, H, t) = 0 (28)

∂u1(x,0, t)

∂z = 0 (29)

u2(L, z, t)

∂x = 0 (30)

u1(0, z, t) =U(z, t) (31)

This boundary value problem has an analytical solution by solving equation (26) by separation of variables. The calculations are not done in this thesis, but the solution can be on the form

u1(x, z, t) =Asin(ωt−kx) cos(lz) (32) provided the dispersion relation is satisfied.

c²= ω²

k²+l² (33)

equation (31) needs to satisfy equations (26), (28) and (29), and a reasonable ansatz is a solution on the same form as equation (32). We assume

U(0, z, t) =Asin(ωt) cos(l(z+B)) (34) where A is the amplitude of the incoming waves, andl and B are determined by the boundary conditions. By inserting equation (34) into equation (29), it is shown thatB = 0 for non trivial solutions. By applying equation (34) into (28) the constants from equation (33) get the values:

lk = π

2h(1 +k)

where k takes the integer values 0,1,2,.. The resulting inflow condition is:

U(z, t) =Asin(ωt) cos(πz

2h(1 +k)) (35)

3.2 Simulations and results

For the convergence tests, we run three simulations with a total simulation time of T=10s, and with equally spaced time and spatial resolutions. We use p1 elements, and the implementation is given in section 9.1. the time and spatial values specified as

• ∆t= 0.01, ∆x= 1/24, ∆z= 1/24

• ∆t= 0.005, ∆x= 1/48, ∆z= 1/48

• ∆t= 0.0025, ∆x= 1/96, ∆z= 1/96

Run L x_s b(x) ∆x ∆z ∆t E_max E_l2n C_max C_l2n

1 2 1 b_l 1/24 1/24 0.01 0.06645 0.02588 - -

2 2 1 b_l 1/48 1/48 0.005 0.02225 0.00970 0.335 0.375

3 2 1 b_l 1/96 1/96 0.0025 0.01647 0.00662 0.740 0.682

1 3 1 b_l 1/24 1/24 0.01 0.06350 0.02390 - -

2 3 1 b_l 1/48 1/48 0.005 0.01614 0.00594 0.254 0.250

3 3 1 b_l 1/96 1/96 0.0025 0.0093 0.00301 0.575 0.520

1 3 1 b_q 1/24 1/24 0.01 0.07274 0.02659 - -

2 3 1 b_q 1/48 1/48 0.005 0.02215 0.00793 0.304 0.298

3 3 1 b_q 1/96 1/96 0.0025 0.0093 0.00354 0.419 0.447

Table 1: Table of numerical results for 3 different simulations. Lis the total length of the domain,x_sis the x coordinate of the boundary between fluid and sponge. b_landb_q denotes the linear and quadratic damping functions used. ∆xand ∆zare the element spacings in the x and z-directions, and ∆tis the time step. E_max andE_l2nare the maximum and L2 norm errors in the simulations, andC_maxandC_l2nare the error reduction rates for the maximum and L2 norm errors with the respect to the previous simulation

We test the sponge by using a linear and a quadratic function each given by

bl(x;Lp) = 10(x−Lp) (36)

bq(x;Lp) = 10(x²−2Lpx+L²_p) (37) The linear function is continuous in the point Lp, and the quadratic function has the function value and the first derivative continous atLp. The valuesk= 0 andω= 10 are chosen, so that the constantslk andkk get the forms:

l0= π 2h k0±

rω² c² − π²

4h²

The three simulations are run with the following domains and damping func- tions.

• L= 2 andLp= 1 with the damping coefficient in equation (36).

• L= 3 andLp= 1 with the damping coefficient in equation (36).

• L= 3 andLp= 1 with the damping coefficient in equation (37).

The results from the simulations are given in figure 2, 3, 4 and table 1.

3.3 Conclusion

An analysis of the scheme shows that when halving the time steps and spatial steps, the maximum error and the L2 norm error from equation (25) should have an error reduction factor around 0.25. Table 1 shows a reduction of the L2 norm and maximum errors, but not with the correct factor. The second simula- tion with a larger sponge layer gives a slightly better result. The L2 norm and

(a)

(b)

(c)

Figure 2: Figure of the errors in the fluid domain for the run withL= 2,x_s = 1 and a linear damping in the sponge layer. (a) shows the errors for the coarse mesh, (b) shows the errors for the finer mesh, and (c) shows the errors for the finest mesh

maximum error is reduced by almost a factor of 0.25 between simulation 1 and 2, but is only reduced by a factor 0.5 between simulations 2 and 3. The errors in with the quadratic damping function are worse than for the linear damping function for the same length of the sponge, The convergence is also worse be- tween the first and second run, but is slightly better between the second and third run. In all cases, it seems that the errors from the sponge become more dominant for better resolutions. figures 2, 3 and 4 show a periodic behaviour

(a)

(b)

(c)

Figure 3: Figure of the errors in the fluid domain forL= 3,x_s= 1 and a linear damping in the sponge layer. (a) shows the errors for the coarse mesh, (b) shows the errors for the finer mesh, and (c) shows the errors for the finest mesh

of the error, indicating that the sponge layer is producing reflected waves with a certain amplitude. In table 2 we have approximated values of the amplitudes from the reflected waves by subtracing the largest and smallest errors in fig- ures and taking the square root2, 3 and 4. The amplitdes are large for poor resolutions, but are reduced with finer resolutions.

(a)

(b)

(c)

Figure 4: Figure of the errors in the fluid domain withL= 3, x_s = 1 and a quadratic damping function in the sponge layer. (a) shows the errors for the coarse mesh, (b) shows the errors for the finer mesh, and (c) shows the errors for the finest mesh

4 The Seismic Wave Equation with Test Solu- tions

In this project, an implementation of the momentum equation will be tested by simple analytic solutions, and the boundary value problem will be simplified by a technique we call test solutions. Assume a rectangular domain Ω of lengthL

Run L x_s b(x) A_r

1 2 1 b_l 0.245

2 2 1 b_l 0.144

3 2 1 b_l 0.116

1 3 1 b_l 0.239

2 3 1 b_l 0.120

3 3 1 b_l 0.074

1 3 1 b_q 0.244

2 3 1 b_q 0.128

3 3 1 b_q 0.075

Table 2: Table of the calculated amplitudes of the reflected waves from the sponge layer in all 3 simulations. L denotes the length of the domain,x_sthe coordinate of the boundary between fluid and sponge, andb_landb_qthe linear and quadratic damping functions respectivley.

H

L

Figure 5: The rectanguar domain used in the problem

and heightH, as given in figure 5. The domain is divided inton×melements in the x and z directions respectively. We assume no body forces in this problem, so equations (2) and (3) reduce to

ρ∂²u

∂t² =∇ ·σ in Ω (38)

σ=λ(∇ ·u)I+µ(∇u+∇u^T) in Ω (39) in the domain. We consider the problem at the times t = t0, t1, ..., tn, and assume that we have an analytic soluion u_e on the whole domain for all t. In the test solution method,ue is applied as initial and boundary conditions. we then have

u(x, z, t) =u_e(x, z, t) at t=t0 (40) u(x, z, t) =u_e(x, z, t) at t=t1 (41) u(x, z, t) =u_e(x, z, t) on Γ (42) By using this method, the need to find more complex solutions by separation of variables or other teqniques are eliminated, and the programs ability to maintain an analytic solution for a given time is tested.

4.1 P and S wave analytic solutions

Known simple solutions of the seismic wave equation are compression and shear waves, denoted as P an S waves. P and S-waves can be divided into further categories as done in Stein and Wysession [2009], but we will concentrate on the coupled P-SV waves in our 2d analysis. A P-wave in the x-z plane can be defined as:

up=Ane^{i(kn·r−ωt)} (43)

where A is the amplitude of the wave, k is the wave number,ω is the angular frequency, t is the time, andris the spatial coordinate vector,

r=xi+zk

andnis the unit normal vector of the wave, given by:

n=nxi+nzk satisfying

|n|= 1 An S-wave in the x-z plane can be defined by:

us=B(n×j)e^{i(kn·r−ωt)} (44)

wherejis the direction along the positive y-axis. The real part of equation (43) is on the form:

u=A(nxi+nzk) cos(knxx+knzz−ωt) (45) And this is a valid solution of equation 4 provided

ω²= (λ+ 2µ)

ρ k² (46)

is satisfied. The real part of the S wave from equation (44) is

u=A(nzi−nxk) cos(knxx+knzz−ωt) (47) and is a solution of equation (4) provided

ω²= µ

ρk² (48)

is satisfied.

4.2 Simulations and results

The program is run with the P and S wave test solutions from equations (45) and (47). The variational form of the problem is given in (17) and we use p1 elements. The implementation is given in section 9.2. For both test solutions,

P θ ∆x ∆z ∆t E_{M ax} E_L2 C_max C_L2 A_r

1 0 1/24 1/24 0.0075 1.71e-7 6.56e-8 - - 0.0003

2 0 1/48 1/48 0.00375 4.07e-8 1.64e-8 0.238 0.250 0.0001 3 0 1/96 1/96 0.001875 9.94e-9 4.09e-9 0.244 0.249 6e-5

1 26.57 1/24 1/24 0.0075 5.41e-7 2.14e-7 - - 0.0005

2 26.57 1/48 1/48 0.00375 1.30e-7 5.41e-8 0.240 0.252 0.0002 3 26.57 1/96 1/96 0.001875 3.19e-8 1.35e-8 0.246 0.250 0.0001

1 71.57 1/24 1/24 0.0075 7.65e-7 2.96e-7 - - 0.0005

2 71.57 1/48 1/48 0.00375 1.83e-7 7.44e-8 0.239 0.251 0.0003 3 71.57 1/96 1/96 0.001875 4.48e-8 1.86e-8 0.245 0.250 0.0001

1 90 1/24 1/24 0.0075 1.71e-7 6.56e-8 - - 0.0003

2 90 1/48 1/48 0.00375 4.07e-8 1.64e-8 0.238 0.250 0.0001 3 90 1/96 1/96 0.001875 9.94e-9 4.09e-9 0.244 0.249 6e-5

Table 3: Table containing the numerical results of the simulations of the seismic wave equation with a P wave test solution. The angleθgives the angle of propagation with the x-axis, ∆xand ∆zgive the element spacings in the x and y direction. ∆tis the time step.

E_max andE_L2 denotes the maximum and L2 norm errors respectvely. C_max andC_L2 are the error reduction rates for the maximum and L2 norm errors with respect to the previous simulation.A_rare the estimated amplitudes from the reflected waves

the lengthL= 1, heightH = 1 and a total simulation time ofT = 5 are chosen.

For the material, the constants λ= 1, µ= 1 and ρ = 1 are used. The wave parameters are A = 1 and ω = 0.5. A convergence test is made by running 3 different simulations for both test solutions with the time and spatial steps evenly distributed

• ∆t= 0.0075, ∆x= 1/24, ∆z= 1/24

• ∆t= 0.00375, ∆x= 1/48, ∆z= 1/48

• ∆t= 0.001875, ∆x= 1/96, ∆z= 1/96

Some results of the simulations are given in tables 3 and 4. the component errors for the p-wave simulation with a propagation angle ofθ= 71.57⁰with the x-axis is given in figure 6. The component errors for an S-wave with a propagation angle ofθ= 71.57⁰ with the x-axis is given in figure 7.

4.3 Conclusion

Tables 3 and 4 show the different simulations for different propagation angles for the P and S-wave test solutions. In all cases the error reduction rates are slightly better than 0.25 which we found in equation (25). From figure 6 we see that the errors in x-displacements are larger in the center of the mesh and close to the corner points, and kept to machine precision at the boundaries. The errors in z-displacements are largest at the center of the mesh, and decreases towards the boundaries, where the error is kept to machine precision. In figure 7, all displacements have their maximum error in the center of the mesh, and decrease towards the boundaries where the errors are kept to machine precision.

In all cases, the errors are kept small, even for the coarsest time and element

(a)

(b)

(c)

(d)

Figure 6: Errors for the x and z-components of displacement for a P-wave with an angle of 71.57⁰ with the x-axis. (a) and (b) show the x and z-displacements for a 24x24 mesh respectively, and a time step of 0.0075. figures (c) and (d) show the x and z-displacements for a 96x96 mesh respectivley, and a time step of 0.001875.

(a)

(b)

(c)

(d)

Figure 7: Errors for the x and z-components of displacement for an S-wave with an angle of 71.57⁰ with the x-axis. (a) and (b) show the x and z-displacements for a 24x24 mesh respectively, and a time step of 0.0075. figures (c) and (d) show the x and z-displacements for a 96x96 mesh respectivley, and a time step of 0.001875.

S θ ∆x ∆z ∆t E_{M ax} E_L2 C_max C_L2 A_r

1 0 1/24 1/24 0.0075 4.48e-7 1.71e-7 - - 0.0004

2 0 1/48 1/48 0.00375 1.07e-7 4.28e-8 0.238 0.250 0.0002 3 0 1/96 1/96 0.001875 2.65e-8 1.07e-8 0.248 0.250 0.0001

1 26.57 1/24 1/24 0.0075 2.81e-6 1.43e-6 - - 0.0012

2 26.57 1/48 1/48 0.00375 6.47e-7 3.49e-7 0.230 0.244 0.0006 3 26.57 1/96 1/96 0.001875 1.60e-7 8.67e-8 0.247 0.249 0.0003

1 71.57 1/24 1/24 0.0075 3.02e-6 1.44e-6 - - 0.0012

2 71.57 1/48 1/48 0.00375 6.98e-7 3.53e-7 0.231 0.245 0.0006 3 71.57 1/96 1/96 0.001875 1.73e-7 8.77e-8 0.248 0.249 0.0003

1 90 1/24 1/24 0.0075 4.48e-7 1.71e-7 - - 0.0004

2 90 1/48 1/48 0.00375 1.07e-7 4.28e-8 0.238 0.250 0.0002 3 90 1/96 1/96 0.001875 2.65e-8 1.07e-8 0.248 0.250 0.0001

Table 4: Table containing the numerical results of the simulations of the seismic wave equation with an S wave test solution. The angleθgives the angle of propagation with the x-axis, ∆xand ∆z give the element spacings in the x and z-direction. ∆tis the time step.

E_max andE_L2 denotes the maximum and L2 norm errors respectivley. C_maxand C_L2 are the error reduction rates for the maximum and L2 norm errors with respect to the previous simulation.A_rare the estimated amplitudes of the reflected waves

Γf

Γd

Γd

Γd

H

L

Figure 8: The problem with test solutions for dirichlet boundary conditions and a given surface stress

spacing. By looking at the tables equation 3, 4, the convergence formula (25) and our choices for ∆x, ∆z and ∆t, we see that the constant C in equation (25) must be smaller than one for the simulations. We also keep in mind that a numerical dispersion analysis has not been made, implying that C could be even smaller. In our simulations, we see that the error has a periodic behaviour, implying that the boundaries are producing reflected waves into the domain.

The amplitudes are estimated by taking the square of the L2 norm errors in tables 3 and 4, and we see that the amplitudes decrease for better resolutions of the mesh.

5 Seismic test solutions with a given stress

In this project, we aim at implementing the seismic wave equation with test solutions, as we did for the previous project, however in this project we apply a given stress to one of the boundaries instead of a given displacement. This gives

insight as to how FEniCS handles boundary integrals and natural boundary conditions. We assume a rectangular domain, as given in figure 8, with the lengthLand heightH. The domain is divided intol×melements in the x and z-directions respectivley. As for the previous project, we neglect body forces for this implementation, giving the the equations of motion and stress:

ρ∂²u

∂t² =∇ ·σ in Ω (49)

σ=λ(∇ ·u)I+µ(∇u+∇u^T) in Ω (50) Again, we assume an analytic solutionue, and solve the problem for the times t=t0, t1, ..., tn. We apply our analytic solution as boundary and initial condi- tions so that

u(x, z, t0) =u_e(x, z, t0) on Ω (51) u(x, z, t1) =u_e(x, z, t0) on Ω (52) u(x, z, t) =u_e(x, z, t) on Γd (53)

σ(u) =σ(ue) on Γf (54)

5.1 P and S-wave analytic solutions

As for the previous project, the P and S-waves from equations (45) and (47) are solutions of the momentum equation provided the dispersion relations from equations (46) and (48) are satisfied respectivley. These solutions are applied as boundary conditions on Γd. On Γs, we apply the given surface stress.

σ_n=n·σ

=k·(σxxii+σxzik+σzxki+σzzkk)

=σzxi+σzzk (55)

The components of stress are found from equation (3) σzx=µ(∂w

∂x +∂u

∂z) σzz=λ(∂u

∂x +∂w

∂z) + 2µ∂w

∂z

(56)

For the P-wave, the components of stress at Γf are:

σzz=−λAk(n²_x+n²_z) sin(knxx+knzz−ωt)

−2µAkn²_ysin(knxx+knzz−ωt) σzx=−2µAknxnzsin(knxx+knzz−ωt)

(57)

And for the S-wave, the components of stress at Γf are:

σzx=µAk(n²_x−n²_z) sin(knxx+knzz−ωt)

σzz=−2µAknxnzsin(knxx+knzz−ωt) (58)

P θ ∆x ∆z ∆t E_{M ax} E_L2 C_max C_L2 A_r

1 0 1/24 1/24 0.0075 1.60e-6 2.77e-7 - - 0.0005

2 0 1/48 1/48 0.00375 3.95e-7 6.61e-8 0.248 0.239 0.0003 3 0 1/96 1/96 0.001875 9.89e-8 1.62e-8 0.249 0.245 0.0001

1 26.57 1/24 1/24 0.0075 4.87e-6 8.49e-7 - - 0.0009

2 26.57 1/48 1/48 0.00375 1.20e-6 2.01e-7 0.246 0.237 0.0004 3 26.43 1/96 1/96 0.001875 2.97e-7 4.91e-8 0.248 0.244 0.0002

1 71.57 1/24 1/24 0.0075 1.11e-6 2.25e-7 - - 0.0005

2 71.57 1/48 1/48 0.00375 2.79e-7 5.66e-8 0.253 0.251 0.0002 3 71.57 1/96 1/96 0.001875 6.93e-8 1.41e-8 0.248 0.250 0.0001

1 90 1/24 1/24 0.0075 5.50e-7 1.81e-7 - - 0.0004

2 90 1/48 1/48 0.00375 1.41e-7 4.68e-8 0.257 0.258 0.0002 3 90 1/96 1/96 0.001875 3.53e-8 1.18e-8 0.250 0.252 0.0001

Table 5: Table containing the numerical results of the simulations of the seismic wave equation with P-wave test solutions. The angleθgives the angle of propagation with respect to the x-axis, ∆xand ∆zgive the element spacings in the x and z direction. ∆tis the time step.E_maxandE_L2denotes the maximum and L2 norm errors. C_maxandC_L2are the error reduction rates with respect to the previous simulation. A_r are the estimated amplitudes of the reflected waves

5.2 Simulations and results

The variational form is given in equation (17) and we use p1 elements. The implementation is given in section 9.3. We run 3 simulations for both the P wave and the S wave test solutions with the length L = 1, height h= 1 and a total simulation time of T = 5. For the material, we choose the constants λ= 1,µ = 1 and ρ= 1. We also choose the parametersA = 1 and ω = 0.5.

The convergence tests are made by varying the evenly distributed element and time spacings

• ∆t= 0.0075, ∆x= 1/24, ∆z= 1/24

• ∆t= 0.00375, ∆x= 1/48, ∆z= 1/48

• ∆t= 0.001875, ∆x= 1/96, ∆z= 1/96

The results for the simulations are given in tables 5 and 6. The component errors for the simulations with an angle ofθ= 71.57⁰with the x-axis are given in figures 9 and 10.

5.3 Conclusion

Tables 5 and 6 show that the error reduction rates for both the P and S-wave test solutions are close to the values estimated from equation (25), yet they are slightly worse than for the previous project for some of the simulations.

Figure 9 shows the x and z-component errors for a wave propagating with an angle of θ = 71.57⁰ with the x-axis. from the figure, 4we see that the larger errors are found at Γf. Local error maximums are also found in parts of the inner domain, while the errors at Γd are kept to machine precision. Figure 10

(a)

(b)

(c)

(d)

Figure 9: Figures of the displacement errors for a P-wave propagating with an angle of θ= 71.57⁰with respect to the x-axis. (a) and (b) show the x and z-displacements for a 24x24 mesh with a time step of 0.0075. (c) and (d) show the x and z-displacement errors for a 96x96 mesh with time step 0.0001875

(a)

(b)

(c)

(d)

Figure 10: Figures of the displacement errors for an S-wave propagating with an angle of θ= 71.57⁰with respect to the x-axis. (a) and (b) show the x and z-displacements for a 24x24 mesh with a time step of 0.0075. (c) and (d) show the x and z-displacement errors for a 96x96 mesh with time step 0.0001875