uI
uR uP
uT uF
θI θR
θP
θT
θT θF
Fluid
Solid Figure 20: The problem for an S-wave hitting the boundary between solid and fluid
which is a form of snell´s law. Inserted into the rest of the boundary conditions, the system of equations determining the amplitudes are found, and given in equation (97):
F =−T e2ik2Hcos(θT) cos(θI)(I−R) + sin(θs)S= cos(θT)(T−F)
k1sin(2θI)(I−R) =kscos(2θs)S
Sµkssin(2θs) =k1(λ+ 2µcos2(θI))(I+R)−κk2(T+F) (97)
Notice that when θI = 0, the system of equations reduces to the results in equations (64), (66) and (68) from the previous project. The system (97) is complicated, and the hand calculations are not done in this thesis. However, numerical methods can be used to solve for the amplitudes by using the complex linear system solver in the scipy module for python, explained by the documen-tation by Jones et al.. To verify the results for the closed system, conservation of energy can be examined by
|E1|=|E1| (98)
where 1 denotes layer 1 and 2 deotes layer 2. In the numerical solver, this equality must be correct to machine precision.
assumed thatcp> cs> cf wherecpis the P-wave velocity in the solid,csis the S-wave velocity in the solid, andcf is the P-wave velocity in the fluid.
uIs=Is(−cos(θIs)i+ sin(θIs)k)ei(k1xsin(θIs)+k1zcos(θIs)−ωt) uRs=Rs(cos(θRs)i+ sin(θRs)k)ei(k1xsin(θRs)−k1zcos(θRs)−ωt) uP s=Ps(sin(θP s)i−cos(θP s)k)ei(kpxsin(θP s)−kpzcos(θP s)−ωt) uT s=Ts(sin(θT s)i+ cos(θT s)k)ei(k2xsin(θT s)+k2zcos(θT s)−ωt) uF s=Fs(sin(θF s)i−cos(θF s)k)ei(k2xsin(θF s)−k2zcos(θF s)−ωt)
(99)
Again, by physical observations it is known that θI =θR and θT = θF. The free surface boundary condition is in this problem also equal to the case with an incoming P-wave, and given from equation (91). Continuity of vertical dis-placement at the internal boundary again forces the angles to follow the type of snell´s law:
k1sin(θIs) =kP ssin(θP s) =k2sin(θT s) (100) The set of equations determining the amplitude ratios are found from the bound-ary conditions (92), (93) and (94) and gives the system of equations determining the amplitude ratios provided I is known in equation (101).
Fs=−Tse2ik2Hcos(θT)
sin(θIs)(Is+Rs) = (Ts−Fs) cos(θT s) +Pscos(θP s) k1cos(2θIs)(Is+Rs) =−Pskpsin(2θP s)
k2κ(Ts+Fs) =PskP(λ+ 2µcos2(θP s)) + (Is−Rs)k1µsin(2θIs) (101) Notice that forθI = 0, the system is unsolvable because no waves are transmit-ted into the fluid, and instead we use the results from the previous project with Ts= 0,Fs= 0, Ps= 0 andRs=Is. We also notice that in this case we have a critical angle at
θ1= kp
k1
where no reflected P-wave is produced at the internal boundary. The equations in (101) are then solved withP = 0. The system is solved in the same manner as for the P-wave solution. Again, the closed system is verified by conservation of energy, giving
|E1|=|E2| (102)
for layer 1 and 2 respectively, and these need to be correct to machine precision when solved numerically.
8 Discussion
At the beginning of this thesis, the goal was to build a model to to solve an earthquake problem and the following P-SV wave propagations in the sea floor,
Ocean Continet
Crust Source
Figure 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.
continent and sea. The projects in this thesis where originally intended to be exercises to test the different parts of the software before a final implementation was attempted. However problems occured in the two-layer model. We have seen that FEniCS handles single domains in a sufficient way by the test solution process. The free surface imposes more errors, but convergence is still main-tained. More difficulties are seen with multiple layers. The sponge layer model has a nice convergence at more coarse resolutions, but this is lost as the reso-lutions improve as the errors from the reflected waves become more dominant.
In The two layer model with vertical incidence, we have nice convergence rates for the P-waves on the solid-solid and solid-liquid problems, but larger errors are found in the solid-liquid boundary. The S-waves have a nice convergence in the solid-solid problem, but we lose convergence for incoming S-waves in a solid-liquid boundary, as large chaotic displacement errors are found in the the fluid domain. In all cases, except the latter, we see a periodic behaviour of the errors, and this shows that the single and multiple layer test-solution process produces small reflected waves at the boundaries. In future researh, a numeri-cal dispersion analysis of the model should be performed to better understand the behaviour of the different simulations. The sponge layer we have used is easily implemented for simple geometries and boundaries found in this thesis, but finding a function b for more complex domains can be very difficult. We discussed another way of implementing the sponge that should be attempted in the future. The two layer model with an incoming S-wave also needs attention, as this does not work with the current implementation. A finite element analysis should be made in FEniCS to better understand the behaviour of the discon-tinuities in the solid-liquid boundary, so the problems can be handled. After such an analysis is made, the problem in section 7 should be implemented and tested. Further research can be made by inverstigating the P-SV wave system in more complex domains. A reasonable goal is then the earthquake model in figure 21, examining the propagation of seismic waves in a realistic problem, and
investigating the full tsunami story that follows. In the end, we would like to remark that although the methods in this work are directed toward seismology, the general theory of the multilayer approach can also be implemented in other aspects of science.
9 Appendix
Below are listings of the codes used in the thesis. The codes are written in python version 2.7.6, and the FEniCS version 1.3. In total, 4 codes have been used. The wave project, the two seismic test solution projects and the project with two layers.
9.1 Code for the sponge layer project
1 f r o m d o l f i n i m p o r t ∗ 2 i m p o r t math a s mt 3
4 d e f s o l v e r ( L , h , x e l , y e l , xs , dt , T , omega , v e l , k , damp , v i z , s a v e ) :
5 ” ” ”
6 F u n c t i o n f o r s o l v i n g t h e s c a l a r wave e q u a t i o n i n a 7 r e c t a n g u l a r domain w i t h g i v e n b o u n d a r y and i n i t i a l 8 c o n d i t i o n s
9 −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
10 INPUT :
11 L : L e n g t h o f domain
12 h : H e i g h t o f domain
13 x e l : Number o f e l e m e n t s p e r u n i t l e n g t h i n t h e x−d i r e c t i o n 14 y e l : Number o f e l e m e n t s p e r u n i t l e n g t h i n t h e y−d i r e c t i o n 15 x s : C o o r d i n a t e o f t h e v e r t i c a l l i n e s e p e r a t i n g f l u i d and s p o n g e
16 d t : Time s t e p
17 T : T o t a l s i m u l a t i o n t i m e 18 omega : A n g u l a r f r e q u e n c y 19 v e l : V e l o c i t y o f waves 20 k : C o n s t a n t d e t e r m i n g i n g t h e
21 damp : l i n f o r l i n e a r , and quad f o r q u a d r a t i c damping 22 v i z : True f o r s h o w i n g s i m u l a t i o n p l o t
23 s a v e : True f o r s a v i n g e r r o r s a t t i m e T
24 −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
25 OUTPUT:
26 R e t u r n s t h e e r r o r b e t w e e n a n a l y t i c and e x a c t 27 s o l u t i o n i n t h e f l u i d l a y e r
28 S a v e s p l o t s o f t h e component e r r o r s i f s a v e=True
29 ” ” ”
30 # S t a r t i n g t i m e
31 t = 0
32
33 # e l e m e n t s p e r l e n g t h
34 l = L∗x e l
35 m = h∗y e l
36
37 # D e f i n e f u n c t i o n s p a c e
38 mesh = R e c t a n g l e M e s h ( 0 , 0 , L , h , l ,m) 39 V = F u n c t i o n S p a c e ( mesh , ”CG”, 1 ) 40 u = T r i a l F u n c t i o n (V)
41 v = T e s t F u n c t i o n (V) 42
43 # D e f i n e s u b d o m a i n s 44 c l a s s F l u i d ( SubDomain ) :
45 d e f i n s i d e ( s e l f , x , o n b o u n d a r y ) : 46 r e t u r n ( b e t w e e n ( x [ 0 ] , ( 0 , x s ) ) ) 47
48 c l a s s Sponge ( SubDomain ) :
49 d e f i n s i d e ( s e l f , x , o n b o u n d a r y ) : 50 r e t u r n ( b e t w e e n ( x [ 0 ] , ( xs , L ) ) ) 51
52 f l u i d = F l u i d ( ) 53 s p o n g e = Sponge ( )
54 d o m a i n s = C e l l F u n c t i o n (” s i z e t ”, mesh ) 55 d o m a i n s . s e t a l l ( 0 )
56 f l u i d . mark ( domains , 0 ) 57 s p o n g e . mark ( domains , 1 ) 58
59 # C r e a t e submesh f r o m f l u i d domain 60 submesh = SubMesh ( mesh , f l u i d ) 61 Vf = F u n c t i o n S p a c e ( submesh , ”CG”, 1 )
62
63 # V a r i a b l e e x p r e s s i o n s 64 c e = C o n s t a n t ( v e l ) 65
66 # S e t t h e damping t o l i n o r quad 67 i f damp==” l i n ”:
68 be = E x p r e s s i o n (” x [ 0 ] < x s ? 0 : 1 0∗( x [ 0 ]−x s ) ”, x s=x s )
69 e l i f damp==” quad ”:
70 be = E x p r e s s i o n (” x [ 0 ] < x s ? 0 : 1 0∗( pow ( x [ 0 ] , 2 )−2∗x s∗x [ 0 ] + pow ( xs , 2 ) ) ”,
71 x s=x s )
72 e l s e:
73 p r i n t ” I n s e r t l i n o r quad ”
74 e x i t ( )
75
76 # D e f i n e i m p o r t a n t c o n s t a n t s 77 s t e p 2 = C o n s t a n t ( 1 / d t∗ ∗2 ) 78 s t e p 3 = C o n s t a n t ( 1 / ( 2∗d t ) ) 79
80 # I n i t i a l c o n d i t i o n s 81 I x y = C o n s t a n t ( 0 ) 82 Vxy = C o n s t a n t ( 0 ) 83
84 # E s s e n t i a l b o u n d a r y c o n d i t i o n s
85 i n f l o w = E x p r e s s i o n (” s i n ( omega∗t )∗c o s ( p i∗x [ 1 ] / ( 2∗h )∗( 1 + k ) ) ”,
86 omega=omega , h=h , k=k , t=t )
87 f r e e = C o n s t a n t ( 0 )
88 d e f s u r f a c e ( x , ob ) : r e t u r n ob and a b s( x [ 1 ]−h )< DOLFIN EPS 89 d e f l e f t f u n ( x , ob ) : r e t u r n ob and a b s( x [ 0 ] ) < DOLFIN EPS 90
91 l e f t = D i r i c h l e t B C (V, i n f l o w , l e f t f u n ) 92 t o p p = D i r i c h l e t B C (V, f r e e , s u r f a c e ) 93 b c s = [ l e f t , t o p p ]
94
95 # S e t a l l f u n c t i o n s i n t o domain 96 c = i n t e r p o l a t e ( c e , V)
97 b = i n t e r p o l a t e ( be , V) 98 u1 = i n t e r p o l a t e ( I x y , V) 99 u2 = i n t e r p o l a t e ( Vxy , V) 100
101 # V a r i a t i o n a l f o r m s
102 F = s t e p 2∗i n n e r ( u , v )∗dx− 2∗s t e p 2∗i n n e r ( u1 , v )∗dx + s t e p 2∗i n n e r ( u2 , v )∗dx +\
103 b∗s t e p 3∗i n n e r ( u , v )∗dx− b∗s t e p 3∗i n n e r ( u2 , v )∗dx +\
104 c∗i n n e r ( n a b l a g r a d ( u1 ) , n a b l a g r a d ( v ) )∗dx 105
106 A = a s s e m b l e ( l h s ( F ) ) 107 u = F u n c t i o n (V)
108 t = 2∗d t
109
110 w h i l e t<= T + 2∗d t + DOLFIN EPS :
111 # P l o t i f v i z=True
112 i f v i z==True :
113 p l o t ( u2 , r a n g e m a x = 1 . 0 , r a n g e m i n =−1.0 , t i t l e =” N u m e r i c a l s o l u t i o n ”)
114 i n f l o w . t = t
115 b e g i n (” Computing a t t i m e l e v e l t = %g ”%t ) 116 LL = a s s e m b l e ( r h s ( F ) )
117 [ bc .a p p l y(A, LL ) f o r bc i n b c s ] 118 s o l v e (A, u . v e c t o r ( ) , LL )
119 end ( )
120
121 u2 . a s s i g n ( u1 ) 122 u1 . a s s i g n ( u ) 123
124 t += d t
125
126 # E x a c t s o l u t i o n
127 l k = mt . p i∗( 1 + k ) / ( 2 .∗h )
128 kk = mt . s q r t ( omega∗∗2/ v e l∗∗2 − l k∗ ∗2 )
129 ue = E x p r e s s i o n (” s i n ( omega∗t − kk∗x [ 0 ] )∗c o s ( l k∗x [ 1 ] ) ”,
130 omega=omega , kk=kk , l k=l k , t=t−2∗d t )
131
132 # I n t e r p o l a t e i n t o f l u i d domain 133 u 2 e = i n t e r p o l a t e ( ue , Vf ) 134 u 2 s = i n t e r p o l a t e ( u2 , Vf ) 135
136 d i f f = T r i a l F u n c t i o n ( Vf ) 137 v f = T e s t F u n c t i o n ( Vf ) 138 l e f t = i n n e r ( d i f f , v f )∗dx
139 r i g h = i n n e r ( u2e , v f )∗dx− i n n e r ( u2s , v f )∗dx 140 l a s s = a s s e m b l e ( l e f t )
141 r a s s = a s s e m b l e ( r i g h ) 142 d = F u n c t i o n ( Vf )
143 s o l v e ( l a s s , d . v e c t o r ( ) , r a s s ) 144
145 # S a v e e r r o r t o f i l e 146 i f s a v e==True :
147 f i l e 1 = F i l e (” e−d−%s−L−%s−h−%s−x e l−%s−y e l−%s−xs−%s−dt−%s−T−%s . pvd ” \
148 % ( damp , L , h , x e l , y e l , xs , dt , T ) )
149
150 f i l e 1 << d 151
152 # R e t u r n t h e a b s o l u t e v a l u e o f t h e e r r o r
153 e r r o r = a b s( d . v e c t o r ( ) . a r r a y ( ) ) 154 r e t u r n e r r o r
155 156
157 d e f r u n s i m u l a t i o n ( ) :
158 ” ” ”
159 T e s t program f o r r u n n i n g an e x p e r i m e n t s h o w i n g t h e 160 p l o t on s c r e e n w i t h g i v e n v a l u e s and r e t u r n i n g t h e 161 e r r o r . The maximum and L2 norm e r r o r s
162 a r e p r i n t e d a t t e r m i n a l
163 ” ” ”
164 L = 3
165 h = 1
166 x e l = 24
167 y e l = 24
168 x s = 1
169 d t = 0 . 0 1
170 T = 10
171 omega = 1 0 .
172 v e l = 1 .
173 k = 0
174 damp=” l i n ”
175 v i z = True
176 s a v e = F a l s e
177 e r r o r = s o l v e r ( L , h , x e l , y e l , xs , dt , T , omega , v e l , k , damp , v i z , s a v e ) 178 e r r o r m a x = e r r o r .max( )
179 e r r o r l 2 n = mt . s q r t (sum( e r r o r∗∗2/l e n( e r r o r ) ) ) 180
181 p r i n t ”Maximum e r r o r : ”, e r r o r m a x 182 p r i n t ” L2 norm e r r o r : ”, e r r o r l 2 n 183
184 d e f t e s t c o n v e r g e n c e ( ) :
185 ” ” ”
186 Program f o r r u n n i n g a c o n v e r g e n c e t e s t w i t h g i v e n p h y s i c a l 187 v a l u e s . The t i m e and s p a t i a l s t e p s a r e h a l v e d t o t e s t t h a t 188 c o n v e r g e n c e i s r e a c h e d . Component e r r o r s a r e t h e n s a v e d t o VTK 189 f i l e s .
190 ” ” ”
191 L = 3
192 h = 1
193 x s = 1
194 T = 10
195 v e l = 1 .
196 omega = 1 0 .
197 k = 0
198 damp = ” quad ”
199 v i z=F a l s e
200 s a v e=True
201
202 # L i s t s t o s t o r e e r r o r v a l u e s
203 E max = [ ]
204 E l 2 n = [ ]
205
206 # L i s t s w i t h dt , dx and dy v a l u e s 207 t i m e s t e p = [ 0 . 0 1 , 0 . 0 0 5 , 0 . 0 0 2 5 ] 208 x e l e m e n t = [ 2 4 , 4 8 , 9 6 ]
209 y e l e m e n t = [ 2 4 , 4 8 , 9 6 ] 210
211 f o r i i n r a n g e(l e n( t i m e s t e p ) ) : 212 d t = t i m e s t e p [ i ]
213 x e l = x e l e m e n t [ i ] 214 y e l = y e l e m e n t [ i ]
215 e r r o r = s o l v e r ( L , h , x e l , y e l , xs , dt , T , omega , v e l , k , damp , v i z , s a v e ) 216
217 e r r o r m a x = e r r o r .max( )
218 e r r o r l 2 n = mt . s q r t (sum( e r r o r∗∗2/l e n( e r r o r ) ) ) 219
220 E max . append ( e r r o r m a x ) 221 E l 2 n . append ( e r r o r l 2 n ) 222
223
224 # Check c o n v e r g e n c e
225 C max = [ ]
226 C l 2 n = [ ]
227 f o r i i n r a n g e(l e n( E max )−1 ) : 228 C max . append ( E max [ i +1]/ E max [ i ] ) 229 C l 2 n . append ( E l 2 n [ i +1]/ E l 2 n [ i ] ) 230
231 p r i n t 40∗’−−’
232 p r i n t ’MAXIMUM ERROR ’
233 p r i n t E max
234 p r i n t 40∗’−−’
235 p r i n t ’ L2 NORM’
236 p r i n t E l 2 n
237 p r i n t 40∗’−−’
238 p r i n t ’CONVERGENCE MAXIMUM ERROR ’
239 p r i n t C max
240 p r i n t 40∗’−−’
241 p r i n t ’CONVERGENCE L2 NORM’
242 p r i n t C l 2 n
243 p r i n t 40∗’−−’
244 245 246 247 248
249 d e f main ( ) :
250 r u n s i m u l a t i o n ( ) 251 #t e s t c o n v e r g e n c e ( ) 252
253 i f n a m e ==’ m a i n ’:
254 main ( )