Stability Analysis of Finite Difference Method in the Propagation of Surface Waves due to Earthquake

Stability Analysis of Finite Difference Method in the Propagation of Surface Waves due to Earthquake

Rekha Tiwarya

aDepartment of Mathematics, RVSCET, Jamshedpur,

India

Anjana P. Ghoraib bDepartment of Mathematics, BIT Mesra, Ranchi,

India

Abstract-In this present context, mathematical modeling of the propagation of Surface waves (Love waves) in a fluid saturated poro-elastic medium has been considered using time dependent higher order finite difference method (FDM). It has been shown that the dispersion curves of Love waves are less dispersed for higher order FDM than of lower order FDM. Stability analysis has been done following conventional Eigen- value method. The variation of stability factor has been derived for r and M (for different order finite difference schemes). It has been shown that the method is always stable as we have considered r 1. It is also shown graphically that it is more stable for lower values of r and for higher order finite difference schemes.

Keywords- Love waves; Fluid saturated porous layer; Time- Space domain; Finite Difference scheme; Accuracy; Dispersion analysis; Phase velocity.

1. INTRODUCTION

   The simulation of surface waves propagating in a fluid saturated poro-elastic media is of great importance to seismologists due to its possible applications in geophysical prospecting, survey techniques and reservoir engineering for understanding the cause and estimation of damage due to natural and manmade hazards.Since poroelastic theory was developed by Biot [1,2], many efforts have been made in using experimental and numerical methods to characterize elastic wave propagation in porous, liquid-saturated solids. The finite difference method is an important tool for numerical simulations of partial differential equations and has been used widely in simulating elastic waves.

   This method depends on the approximation of the temporal and spatial derivatives present in these equations. Also the finite difference method has greater advantage over the other purely numerical methods for using relatively small memory and computation time. It specifies the model as a series of grid points and approximates the spatial and temporal derivatives by using the model values at nearby grid points. To approximate temporal derivatives, a second order finite difference scheme is generally used to limit the accuracy of modeling. A smaller time step or grid size may increase the modeling accuracy but also increase the computation time. But the availability of the modern digital computers has overcome this difficulty and the

   efficient and more general finite difference method has been smoothly developed. Therefore higher order finite difference on the space derivative appears to be a popular method to increase modeling accuracy.

   Since poro-elastic theory was developed by Biot[1,2], many efforts have been made in using experimental and numerical methods to characterize elastic wave propagation in porous, liquid-saturated solids. The finite difference method (FDM) is an important tool for numerical simulations of partial differential equations and has been used widely in simulating elastic waves. Quite a good amount of information about the numerical modeling and propagation of seismic waves using FDM is available in the literature of many authors, namely, Boore [3], Alford [4]. Virieux [5] have used velocity-stress finite difference method for the propagation of SH wave in heterogeneous media. To improve the accuracy and stability of FDM many authors has used and developed different types of difference schemes. Levander [6] applied 4th order approximation in space to the P-SV scheme. Kristek et al.

   [7] considered seismic wave propagation in visco-elastic media using 3D forth order staggered-grid finite difference scheme. Kristek et al. [8] considered 1D elastic problem on the accuracy of the finite-difference schemes. Tessmer [9] discussed Seismic finite difference modeling with spatially variable time steps. Finkelstein et al. [10] developed finite difference time domain dispersion reduction schemes. Y. Liu et al. [11, 12] considered an implicit finite difference scheme for seismic modeling and also considered a new time-space domain high-order difference method for the acoustic wave equation. Lie et al. [13] discussed finite difference numerical modeling in two phase anisotropic media with even order accuracy. Zhu et al. [14] developed finite difference modeling of the seismic response of fluid saturated, porous, elastic solid using Biot theory.

2. FORMULATION OF THE PROBLEM

   We consider a model consisting of water saturated anisotropic poro-elastic layer of finite thickness H ; the z axis are taken vertically downward. The x axis are chosen

   parallel to the layer in the direction of propagation of surface wave (Fig.1).

   The equations of motion for the fluid-saturated anisotropic porous layer without body force and neglecting the viscosity of the fluid are

   2

   a

   where, d

   1/ d

   2 / and

   (5)

   N / ' ,

   u U

   11 12 22 a

   ij, j

   t 2

   1. i

      12 i

      (1a)

      the velocity of the shear wave in the corresponding non- porous anisotropic elastic medium along the direction of x.

      2

      u U

      Also,

      11 11 / ,

      12 12 / ,

      22 22 / .

      (6)

      , i t 2

   2. i

   22 i

   (1b)

   are the non-dimensional parameters for the material of the porous layer .

   where

   i j

   are the components of stress tensor in the

   To improve the accuracy, we have considered the higher order finite difference scheme for spatial derivatives

   solid skeleton, ( f p) is the reduced pressure of the as

   fluid ( p is the pressure in the fluid , and f is the porosity

   2 v 1

   0

   0

   a v

   • a

     (v0

     v0 )

     of the porous layer),

     ui are the components of the

     x2

     p

     0 0,0

     m

     M

     M

     m1

     m,0

     m,0

     displacement vector of the solid and Ui are those of fluid.

     (6a)

     The dynamic coefficients,

     11 ,

     12 ,

     22

     take into

     2 v 1

     0

     0

     a v

   • a

     (v0

     v0 )

     account of the inertia effects of the moving fluid and are

     z 2

     p

     0 0,0

     m

     M

     M

     m1

     0, m

     0,m

     related to the densities of the solid s , the fluid f

     the layer by .

     and

     (6b)

     As generally higher order finite difference scheme on temporal derivatives requires large space in the

     Using the conventional Love wave conditions,

     computer memory and usually unstable, 2nd order finite

     ie.

     u (0, u

     ,0)

     and

     (0,U

     ,0)

     difference scheme is used for temporal derivatives as:

     y U y

     2v

     1 0

     1 1

     where

     ,

     uy v (x, z,t)

     and

     U y V (x, z, t)

     t 2 2

     2v0,0 (v0,0 v0,0 )

     v

     v

     (6c)

     and solving we have

     2

     2

     N v

     2 v

     2 v

     where

     n m, j

     v (x mh, z jh , t n ) ,

     x2 L z 2

     d t 2

     h is the grid size and is the time step.

     Using equations (6a), (6b), (6c) into equation (4) we have

     (2) a (1 )v0 M a (v0

   • v0

   • v0

     )

     where

     d 2 /

     0 0,0 m

     m,0

     m,0

     0, m

     0,m

     11 12 22

     (3)

     p

     m 1

     0 1

     1

     where N , L correspond to the familiar Lames

     2 2

     2v0,0

     (v0,0

     v0,0 )

     constants, From the equation (2) it can be seen that the velocities of shear waves in the porous medium in x and

     (7)

     z directions are

     N / d

     and

     L / d

     respectively .

     Using the plane wave theory, let us consider

     v

     v

3. FINITE DIFFERENCE APPROXIMATION

   n m, j

   e i[kx ( xmh)kz ( z jh) (tn )]

   (8)

   Equation (2) can be written as

   Substituting equation (8) into equation (7) and simplifying, we have

   2 v 2 v 2 v

   1 M

   x2

   z 2

   2 t 2

   a0 (1 ) am cos(mkhcos ) cos(mkhsin )

   2 m 1

   Where,

   (4)

   N / L , the anisotropic parameter and

   r 2 1 cos( )

   (9)

   N / d '

   , the velocity of shear wave in the

   where,

   r

   and

   porous medium in x direction.

   The shear wave velocity in the x-direction may be expressed as

   h

   kx k cos , kz k sin , and , being the propagation direction angle of the plane wave.

   Using the Taylor series expansion for cosine functions, we have from the equation (9)

   the eigen-value equation

   2 g 1 0 will be less

   1 M

   [ cos 2 j sin 2 j ](mkh )2 j

   than or equal to 1 if g

   2.

   j

   j

   a0 (1 ) am (1 ) (1)

   Since the error generally increases with the increase of the

   2 m 1

   j 1

   (2 j) !

   wave number, lets consider the maximum wave number(

   r (kh)

   r (kh)

   2 j 2 2 j

   (1) j

   Nyquist frequency) as

   j 1

   (2 j)!

   (10)

   kx kz

   .

   h

   Comparing the coefficients of k 2 j , we get,

   M

   M

   a0 2am 0

   m1

   (17)

   Using the equation (17) into the equation (16), we have

   4r 2 (1 ) M1

   M

   M

   a ( cos2 j sin2 j )(m)2 j r2 j 2 ,

   (11)

   j 1,2…….M

   g 2

   a2m1

   m1

   (18)

   m

   m 1

   Where M1 int[(M 1) / 2], int[] is a function to get

   This equation indicates that the coefficients

   (12)

   am are the

   the integer part of a value. Therefore the stability condition is

   function of . To obtain a single set of coefficients, an

   4r 2 (1 ) M1

   optimal angle has to be chosen. We solve the equation (12) 2

   a2m1 m1

   2,

   to get

   am by using

   / 4 and then

   a0 can be

   (19)

   obtained from equation (11).

   (1 ) M1 2

4. STABILITY ANALYSIS

   i.e.,

   r

   a2m1

   The recursion equation of finite difference scheme can be obtained from the equation (7) as follows:

   m1

   1

   1

   (20)

   As a particular case if we take 1(the case of isotropic

   v

   v

   1

   1

   1

   0,0

   medium)), the stability condition is reduced to

   (1 ) r 2

   a

   2v0

   • r 2

     a (v0

     m

     m

   • v0

   ) (v0

   v0

   ) v1

   M1 2

   0

   0

   0,0

   M

   M

   m 1

   m,0

   m,0

   0, m

   0, m

   (13)

   0,0

   r 2a

   2m1

   Using the conventional eigen-value method of stability

   m1

   (21)

   analysis, lets consider:

   which is discussed by Liu and Sen [12].

   To calculate and analyze the stability of the finite

   p

   p

   0

   m,m

   U

   U

   0

   m,m

   0

   v

   v

   p

   p

   ;

   ;

   m,m

   0

   m,m

   0

   q

   q

   m,m

   , q

   , q

   0

   m,m

   1

   v

   v

   ;

   ;

   m,m

   T W 0 ei(kx m h kz m h) ;

   difference scheme, we define the stability factor s

   according to the equation (20) as follows:

   1

   U 1 p1

   , q1

   T W 1ei(kx m h kz m h)

   (1 ) M1 2

   m,m

   m,m

   m,m

   s

   a2m1

   (14)

   m1

   (22)

   Using equation (14) in equation (13), we obtain

   We calculate the variation of s with r and M .

   1

   1

   W 1 GW 0 g

   1W 0

   0

   0

   (15)

5. DISPERSION ANALYSIS

   Let us define a parameter to describe the dispersion of the Finite difference scheme by using

   Where G is the transition matrix and

   equation (9) as follows:

   g 2

   2r 2

   M

   M

   am [ cos(kx mh) cos(kz mh) ( 1)]

   m1

   (16)

   vFD

   a

   2 sin 1

   2 M

   r

   r

   a

   2 mkhcos

   2 mkhsin

   The recursion relations of finite difference scheme will be stable if the absolute values of the eigen-values of the transition matrix are less than or equal to 1. The roots of

   rkh d

   m1

   m sin

   sin

   2 2

   (23)

   If is equal to 1, then there is no dispersion. However, if

   is far from 1, a large dispersion will occur. Also kh is equal to at the Nyquist frequency, so in calculating , kh only ranges from 0 to and the variation of is

   To analyze the stability of the finite difference scheme, we calculate the variation of stability factor s from the equation (22) for different values of r 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 and M 6,8,10,12.

   from 0 to

   / 4 .

   Figs. 2-4 display the dispersion curves of Love waves at different values of M in a homogeneous non

6. NUMERICAL CALCULATION AND DISCUSSIONS WITH RESPECT TO GRAPHICAL

   REPRESENTATIONS

   The numerical calculation of the equation (23) has been done for different values of the parameters , d and

   by taking, 3000 . The phase velocity vFD / a of

   porous elastic solid, porous isotropic and anisotropic layer respectively. It is found that dispersion is more for the lower values of M and decreases for higher values of M. Here it is also observed that the increase in porosity leads to the decrease in the magnitude of the phase velocity of Love waves and an increase in anisotropy leads to increase in the phase velocity of Love waves

   Fig. 5 displays that the area for stable recursion

   Love wave from the equation (23 ) verses kh has been computed for different values of

   decreases with the increase of M.

   Figs. 6 and 7 displays the variation of stability factor s

   for

   d 1, 0.9, 0.8, 0.7 , 1, 2, 3, propagation

   anangglele 00,, //1616,, / 8, 3 /16, / 4 and

   different values of r and M . It has been shown that the

   method is always stable and more stable for lower values

   of r and higher values of M .

   t 0.0005, 0.001, 0.0015, 0.002, 0.0025..

   Rigid layer

   z = 0 X

   fluid saturated anisotropic poro-elastic medium

   Z

   Fig. 1: Geometry of the problem.

7. CONCLUSION

It is observed that the higher order time dependent Finite Difference Method plays an important role in the propagation of Love wave in a porous layer. Graphically, we have shown that the dispersion curves of Love waves are less dispersed for higher order Finite Difference Method. The results are copared for different order finite difference scheme (M= 6, 8, 10, 12). Here it is also observed that the increase in porosity leads to the decrease in the magnitude of the phase velocity of Love waves and

an increase in anisotropy leads to increase in the phase velocity of Love waves

Stability analysis has been done following conventional Eigen- value method. The variation of stability factor has been derived for r and M (for different order finite difference schemes). It has been shown that the method is always stable as we have considered r 1 . It is also shown graphically that it is more stable for lower values of r and for higher order finite difference schemes.

M=4

M=6

M=4

M=6

1

M=8 M=10 M=12

M=8 M=10 M=12

0.98

0.96

0.94

vFD

a

0.92

0.9

0.88

0 0.5 1 1.5 2 2.5 3 3.5

Fig.2 Dispersion curves for Love waves for different values of M when d = 1 and = 1.

vFD

a

1.14

M=4

M=4

M=6 M=8 M=10 M=12

M=6 M=8 M=10 M=12

1.12

1.1

1.08

1.06

1.04

1.02

1

0.98

0 0.5 1 1.5 2 2.5 3 3.5

Fig.3 Dispersion curves for Love waves at different values of M when d = 0.8 and = 1.

M=4

M=4

1.26

M=6 M=8 M=10 M=12

M=6 M=8 M=10 M=12

1.24

1.22

1.2

1.18

vFD

a

1.16

1.14

1.12

1.1

1.08

0 0.5 1 1.5 2 2.5 3 3.5

kh

Fig.4 Dispersion curves for Love waves for different values of M when d = 0.8 and = 2.

0.6

0.59

0.58

s for r 0.1

s for r 0.2

s for r 0.3

0.57

0.56

0.55

0.54

0.53

0.52

0.51

s

s for r 0.4

0.54 5 6 7 8 9 10 11 12

M

Fig. 5 The variation of stability factor s verses M for different values of r .

0.66

0.64

0.62

0.6

0.58

0.56

0.54

s for M=6

s for M=8

s for M=10 s for M=12

0.52

s

0.5

0.1 0.2 0.3 0.4 0.5 0.6

r

Fig 6. The variation of stability factor s verses r for different values of M .

0.8

0.7

s = r

s for M=12

0.6

0.5

0.4

0.3

0.2

s

0.1

0.1 0.2 0.3 0.4 0.5 0.6

r

Fig 7 The variation of stability factor s verses r for different values of M .

REFERENCES

1. Biot M. A. (1956a) Theory of propagation of elastic solid in a fluid-saturated porous solid. 1. Low-frequence range, J. Acoust. Soc. Am.28,168-178.

2. Biot M. A. (1956b) Theory of propagation of elastic solid in a fluid-saturated porous solid. II. Low-frequence range, J. Acoust. Soc. Am.28,179-191.

3. Boore, D.M. (1972) Finite difference methods for seismic wave propagation in heterogeneous materials, Methods in Computational Physics, Vol. II, Academic Press, New York.

4. Alford R. M., Kelly K. R. and Boore D. M. (1974), Accuracy of finite difference modeling of the acoustic wave equation, Geophysics, 39(6) 834-842.

5. Virieux J. (1984), SH-wave propagation in heterogeneous media: velocity stress finite difference method, Geophysics 49, 1933-1957.

6. Levander, A. (1998) Fourth-order finite difference P-SV seismograms: Geophysics, 53,1425-1436.

7. Kristek J and Moczo P.(2003) Seismic wave propagation in visco- elastic media with material discontinuities: A 3D forth-order

   staggered-grid finite difference modeling, Bull. Seism. Soc. Am. 93(5) 2273-2280.

8. Kristek J and Moczo P.(2006) On the Accuracy of the Finite- Difference Schemes: The 1D elastic problem, Bull. Seism. Soc. Am.96: 2398-2414.

9. Tessmer E.(2000) Seismic finite difference modelingwith spatially variable time steps,

   Geophysics, 65, 1290-1293.

10. Finkelstein B. and Kastner R.(2007) Finite difference time domain dispersion reduction

    schemes, J. Comput. Phys., 221, 422-438.

11. Liu Y. and Sen M. K.(2009a) A practical implicit finite difference method: examples from seismic modeling:, J. Geophys. Eng, 6 231- 249.

12. Liu Y. and Sen M. K.(2009b) A new time- space domain high- order finite difference method for the acoustic wave equation., J. Comput. Phys.,228, 8779-8806.

13. Liu Y. and Xiucheng W.(2008) Finite difference numerical modelingwith even order accuracy in two phase anisotropic media, App. Geophysics, 5, 107-114.

14. Zhu X. and McMechan G. A., (1991) Finite difference modeling of the seismic response of fluid saturated, porous,elastic solid using Biot theory, Geophysics 56,424-435.