S T A B I L I T Y PROPERTIES OF T H E DIRECT B O U N D A R Y E L E M E N T M E T H O D A P P L I E D TO T H E E L A S T O D Y N A M I C E Q U A T I O N S by S H E L L Y D E L E E N PINDER B . S c , University of Lethbridge, 1994 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T OF M A T H E M A T I C S We accept this thesis as conforming to the)required standard U N I V E R S I T Y OF BRITISH C O L U M B I A August 2000 ® Shelly Deleen Pinder, 2000 In presenting t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I f u r t h e r agree that permission f o r extensive copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the head of my department or by h i s or her r e p r e s e n t a t i v e s . I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n permission. Department of fY]r;|-f~rT€ lYKCA\~l C - S The U n i v e r s i t y of B r i t i s h Columbia Vancouver, Canada Date Afft. 1 ? JOnO Abstract The elastodynamic equations are frequently encountered in the geosciences for assessing the stability of mining excavations or performing seismological analyses of rock bursts and earthquakes. However, the geometric complexity of the problems encountered usually prevents exact solutions from being determined. Thus it is necessary to resort to numerical approximations in order to obtain solutions to elastodynamic problems with general geometries. In previous work, using boundary element methods for approximating solutions of the elastodynamic equations has been limited due to numerical instabilities that appear sporadically when applying this technique. Until recently, finite difference and finite element methods have been used almost exclusively for these approximations. Besides being computationally more expensive since the entire domain needs to be discretized, finite difference and finite element methods also have problems with numerical dispersion. Examining boundary element methods could lead to the development of techniques that are less intense computationally and avoid any numerical dispersion problems. For boundary element methods to become a standard tool for approximating these solutions, there must be some type of criterion established for choosing meshing parameters to ensure stability. A complete analysis of a one-dimensional model problem is performed via the z - transform. For this model problem, the validity of the stability analysis is confirmed from a comparison of the analytic results with numerical experiments. This allows some guidelines to be made to ensure that a particular numerical approximation is stable when applied to the model problem. i i Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgement viii Chapter 1 Introduction 1 Chapter 2 Formulation of the Problem 3 2.1 General Boundary Integral Formulations 3 2.2 The Model Problem 5 Chapter 3 Discretization Schemes 9 3.1 The Forward Euler Method 10 3.2 The Backward Euler Method 11 3.3 The Midpoint Rule 13 3.4 The Trapezoidal Rule 15 3.5 Simpson's Rule 16 3.6 The e -Scheme 18 3.7 The Half-Step Scheme 20 Chapter 4 Stability of the Model Problem 24 4.1 Stability Analysis Using the z-transform 24 4.2 The Forward Euler Method 25 4.3 The Backward Euler Method 27 4.4 The Midpoint Rule 28 4.5 The Trapezoidal Rule 28 4.6 Simpson's Rule 30 4.7 The €-Scheme 31 4.8 The Half-Step Scheme 33 4.9 Summary 36 Chapter 5 Numerical Results 37 5.1 Intermittent Instabilities 37 5.2 Numerical Experiments 38 5.3 The Forward Euler Method 39 5.4 The Backward Euler Method 40 i i i 5.5 The Midpoint Rule 4 1 5.6 The Trapezoidal Rule 4 3 5.7 Simpson's Rule 4 4 5.8 The e-Scheme 4 6 5.9 The Half-Step Scheme 4 8 5.10 Summary 5 0 Chapter 6 Conclusions 51 Bibliography 5 3 iv List of Tables 1 Analytic stability classification 36 2 Evidence of intermittent instabilities 37 3 Experimental stability classification 50 v List of Figures 1 The model problem 6 2 Exact solutions (u= l , p = l , X=0) 8 3 The Forward Euler method 10 ,T - | -4 Approximation of J o p(0,T)dT using the Forward Euler method 11 5 The Backward Euler method 12 ,T - | -6 Approximation of J0 p(0,T)dT using the Backward Euler method 12 7 The Midpoint rule 13 8 Approximation of Jo p(0,T)dT using the Midpoint rule 14 9 The Trapezoidal rule 15 10 Approximation of Jo p(0,T)dT using the Trapezoidal rule 15 11 Simpson's rule 16 T-2. o p(0,T)dT using Simpson's rule 18 13 Thee-scheme 19 14 Approximation of J o p(0,T)dT using the 6-scheme 20 15 The half-step method (n even) 20 16 Approximation of Jo p(0,T)dT using the half-step method (N and k even) . . 21 17 Poles of T F(z) (At = 0.06 => k=23,Y=0) 26 18 Poles of TT(z) (At = 0.06=> k=23,Y=0) 29 19 Poles of T s(z) (At = 0.05 => k=28,Y=0) 31 20 Poles of T e (z) (€ = 0.2, At = 0.05 => k=28,Y=0) 32 21 Poles of T e(z) (€ = 0.35, At = 0.05 = • k=28, Y=0) 32 v i 22 Poles of T H (z) (At = 0.06 => k=23, Y =0) 35 23 Approximation using the Forward Euler method (At = 0.01 ==> Y^O) . . . . 39 24 Approximation using the Forward Euler method (At =*• Y=0) 39 25 Approximation using the Backward Euler method (At = 0.01 =*• Y^ 0) • • • • 40 26 Approximation using the Backward Euler method (At = T ^ =* Y=0) . . . . 41 27 Approximation using the Midpoint rule (At = 0.01 => Y# 0) 42 28 Approximation using the Midpoint rule (At => Y=0) 42 29 Approximation using the Midpoint rule (At =* Y= 0, T=50) 43 30 Approximation using the Trapezoidal rule (At = 0.01 ==> Y^ 0) 43 31 Approximation using the Trapezoidal rule (At =>• Y=0) 44 32 Approximation using Simpson's rule (At = 0.01 = • Y# 0) 45 33 Approximation using Simpson's rule (At =* Y=0) 45 34 Approximation using the 6 -scheme (e= 0.2, At = 0.01 => Y^ 0) 46 35 Approximation using the € -scheme (6= 0.35, At = 0.01 =>• Y^0) 46 36 Approximation using the € -scheme (6= 0.2, At = ^ Y=0) 47 37 Approximation using the 6 -scheme (€= 0.35, At = y ^ : => Y=0) 47 38 Approximation using the half-step scheme (At = 0.01 =>• Y# 0) 49 39 Approximation using the half-step scheme (At =T^" ^ Y=0) 49 v i i Acknowledgement I would like to thank my supervisor, Dr. Anthony Peirce, for his guidance over the last couple of years. His assistance in suggesting a topic and helping me work through the research process is greatly appreciated. As always, my parents have been relentless in their support and encouragement. Thank you for being there for me. The last people who I would like to recognize have been invaluable to me. M y good friends David Ferguson and Seema A l i have been instrumental in my accomplishment. Their support has been constant, whether it be reading over my latest work or just relaxing. Thanks guys! vm Chapter 1 Introduction In the geosciences, the solution of the elastodynamic equations is important for a number of applications such as fracture mechanics, ice mechanics, seismology, soil mechanics, and structural mechanics. Geological prospecting and assessing the dynamic effects of stress waves on surface structures and excavations, excavation support, and fractures are just a few examples where the elastodynamic equations are essential in the modelling process. These problems are almost always too difficult to solve analytically so a numerical solution must be found. Finite element and finite difference methods have been used for the most part in the past, but the boundary element method offers an attractive alternative for numerical approximations. As opposed to finite element or finite difference methods in which the entire domain of the problem is discretized, only the boundary is discretized in the boundary element method. This reduction of the dimension of the problem brings about a substantial reduction in computational costs. In addition, when modelling such complex phenomena as fracture propagation, computational costs decrease since there is no need to re-mesh the domain at each growth increment in fracture growth models. The boundary element method is based on two different types of integral equation formulations which are described in detail in Chapter 2. The direct boundary element method is formulated in terms of physical quantities such as displacements and tractions while the indirect boundary element method involves fictitious quantities which have no physical significance. Stability analysis has thus far been done only on indirect boundary element equations [11], so the object of this thesis is to analyse the stability properties of the direct boundary element method. Boundary element methods have not been used extensively for numerical solutions in elastodynamics since they seem to work well for some problems but have exhibited instabilities for others. In the literature, numerical results have often been presented for a given time-step At with no description of the consequences of changing this time-step. These results can be misleading since intermittent instabilities may occur i f the time-step is altered slightly. One would think that decreasing At can only improve the stability of a given numerical approximation. This is not always the case. A problem has intermittent instabilities i f it is unstable for a given At, becomes stable as At increases, and goes unstable once more as At continues to grow. The existence of intermittent instabilities has long been debated and even ignored, but growing evidence of instabilities in boundary integral elastodynamic models [1, 11] provides motivation to investigate the cause of this phenomenon. Understanding why these instabilities occur is an initial step in designing discretization procedures with improved stability characteristics. With instabilities occurring in a seemingly random manner, boundary element methods have not been a feasible alternative for numerical approximation of the elastodynamic equations. The objective of this thesis is to use the model problem as given by Frangi and 1 Novati [6] to study the stability of various discretization schemes for the direct boundary element method. In Chapter 2, the general boundary element formulation is described and the one-dimensional model problem is presented. In Chapter 3, various discretization schemes are applied to the model problem. In Chapter 4, results obtained from stability analysis performed via the z -transform on the model problem are presented. The numerical results obtained by implementing the discretization schemes for approximating the solution of the model problem are given in Chapter 5. In Chapter 6 we present some concluding remarks. Even for this model problem with a simple geometry, intermittent instabilities appear. The model problem is useful because the exact solution can be solved for analytically and compared with the numerical approximations. Real-life problems are often too complicated to diagnose the numerical difficulties. If a problem is too complex, it is difficult, i f not impossible, to isolate the cause of the instabilities which limits the applicability of the results. Stability analysis is performed to try and predict when an instability wi l l occur. The ultimate goal is to discover what causes the instabilities and to develop some criteria for choosing a stable discretization scheme and appropriate time-step for the model problem. Analysis of a model problem, such as the one examined here, provides a framework to design discretization schemes with improved stability characteristics. Future work on problems with more complicated geometries can then benefit from the insight gained from the model problem. 2 Chapter 2 Formulation of the Problem 2.1 General Boundary Integral Formulations The general equations of linear elastodynamics are the equation of motion, the kinematic equation, and the constitutive equation [4, 9]. These are given respectively as Ojjj + pbj = piij (2.1a) S,y = y ( " y + Uji) (2.1b) a ij = Xzkk5 ij + 2fie ij (2.1c) where Uj(x, t) is the displacement vector at point x and time t, ay and ey are the stress and strain tensors, and bj is the body force per unit mass. The Lame constants are X and p, the mass density is given by p, and Sy is the Kronecker delta function. Spatial differentiation is indicated by commas and dots represent time derivatives. In three-dimensions, ij = 1,2,3. The equations for the one-dimensional model problem can be reduced from the general three-dimensional boundary integral equations. These general equations are based on fundamental solutions and the dynamic generalization of Betti-Rayleigh's reciprocal theorem. In dynamics, the reciprocal theorem specifies a relationship between a pair of elastodynamic states. With no body forces present and given zero initial conditions, the reciprocal theorem is given by [ tt(x,t-s) * u'i(x,s)dS(x) = [ t'Xx,t-s) * Ui(x,s)dS(x) (2.2) J s J s where and w,- are tractions and displacements in one elastodynamic state and t\ and u\ are the tractions and displacements in the other elastodynamic state. One state is chosen to be the fundamental solution to a unit impulsive loading and the other unknown state is to be found from the boundary integral equations. The symbol * is the time convolution operator. The fundamental solution is given by Stokes' displacement solution to a concentrated unit pulse applied at time t = 0. That is, 3 where and 8jj is the Kronecker delta function. The tensor Uy represents the z'th component of displacement at point x and time t due to a force at point | and time s = 0, parallel to the j'th axis, and having a time dependent magnitude off[t). B y substituting the fundamental displacements defined by equation (2.3) into Hooke's law (2.1c), the fundamental stresses are obtained and are given by Combining the dynamic reciprocal theorem with the appropriate fundamental solutions, the direct boundary element method equations are obtained. Without loss of generality, assume zero initial conditions. In the absence of body forces, the direct boundary integral equations are There are two types of indirect boundary element methods : the fictitious stress method and the displacement discontinuity method. Both of these can be derived from the direct boundary integral method by adding an interior to an exterior domain problem and subtracting the interior region equations from those for the exterior region. Assuming that the displacement jumps across the boundary between interior and exterior domains are zero, the fictitious stress method is obtained. The displacements and stresses respectively are (2.4) " * ( £ 0 = j [Uik(x,t;£,0) * ti{x,t) - TugU,t,£,0)nj * Ui(x,t)]dS(x) (2.5) 4 « * ( * , / ) = f Uu(x.,t;£,0) * F,($,t)dS(& J s ~ ~ okj(x, t) = f Tkji(x, t; £ 0) * FM, t)njdS^) (2.6) J s where F, = ft - tj are the traction jumps across the fictitious stress surface S. If the tractions are required to be continuous across the boundary, then the displacement discontinuity method is obtained. The displacements and stresses here are given by uk(x,t) = J J * * , ? 0 ) w ; * Di(£,t)dS(§) °kj(x, t) = f Sklij(x, t; & 0)nj * Dtf, t)dS($) (2.7) * s ~ ~ where £), = u\ -uJ are the displacement jumps across the displacement discontinuity surface S, defined by The direct boundary integral method is defined in terms of the displacements and stresses given by (2.3) and (2.4). For the indirect boundary integral method, the equations are given in terms of the traction jumps Ft or the displacement jumps Ti. Both of these variables have no physical significance and are called fictitious quantities. Using a model problem, the stability of the direct method wil l be studied in this thesis. 2.2 The Model Problem The model problem of interest [6] has a relatively simple geometry and can be solved analytically, so the results obtained by using various numerical approximation schemes can be compared to the exact solution. This gives insight to the strengths and weaknesses of the different approximation methods which wil l be discussed in Chapter 3. Consider an infinite strip of height L = 1 as in Figure 1. A uniform traction of P = H(t) is applied to the upper surface, where H is the Heaviside function. Let u be the displacement in the X2 - direction. The displacement u and the traction ux are independent of x i and depend only on X2 so the problem is one-dimensional in space. With this being the case, let X2 = x to simplify notation. Then the governing equation is (A + 2 / z ) ^ 0 _ = Q 5 with the boundary conditions p(l,t) = (l + 2u)^^- = 1, M(0,0 = 0 (2.10) where X and p. are the Lame constants and p is the material density. Equation (2.9) can be rewritten as Uxx - \ u u = 0 where The quantity c can be interpreted as the propagation velocity of pressure waves. (2.11) (2.12) mmtmmmm P = H(t) L = 1 L. ///////////////////////////////////// Figure 1 : The model problem The Green's function for the model problem can be determined by taking the Laplace transform of - 7T G « - Y+~TU (2-13) and solving for the constants using the continuity and jump conditions. The displacement atx at time t due to a concentrated impulsive force acting at x at time r is given by G(x,x, t-x)= 1 { X C + 2 ^ H[c(t -r)-r] (2.14) where 6 r = |x - x|, (2.15) c is the longitudinal wave velocity as given in (2.12) and / / i s the Heaviside function. Multiplying (2.9) by the Green's function and integrating over the space-time interval [0,1] x [0,r], the following identity is obtained : 0 [(A + 2n)uxx - pult]G(x,x,t - r)dxdr „ o J : 0, - j " ^ u(x,t)8(x - x)8(t - r)dx + p(l,r)G(x, l,t - r) +p(0, r )G(x, 0, t - r) - (A + 2u)u( \,T)Gx(X,\J - z) dr (2.16) wherep(0,t) = -(X + 2jj)ux(0,t). Then ' 0 J 0 O u(x, r)S(t - r)r5(x - x)dxdt „ o f' [p( 1, T)G(X, \,t-r)+ p(0,r)G(x,0,t-r)-(X + 2u)u{ 1, r)Gx(x, 1, / - T ) 1 * . (2.17) From this equation, u(x, t) can be expressed as u(x,t) = \u(\,t-^-) + 2(A + 2|/) where ro and r\ are the distances of x from x = 0 and x = 1 respectively. [ ri-— (t-l±)H(ct-ri) + H(ct-ro)\ ' p(0,r)dt (2.18) Using the direct boundary integral method, a system of two integral equations with two unknowns, u(\,t) andp(0,t), is obtained by evaluating (2.18) at x = 0 and x = 1 : r 1 ' ( 2 - 1 9 ) " (1 ,0 = Tf2j[tH(ct)+H(ct-l)jo cp(0,T)drj This system of equations can be solved for exactly by looking at intervals of size At = Y . In any interval n \ 1 ) where n = 0,1,2, • • •, p(0, t) can be solved for using the first equation of (2.19) and then u(\,t) can be found from the second equation of (2.19). At this point,p(0,t) and u(l,t) are known for every value of t in the interval /, so this information can be used to solve forp(0, t) and «(1 , t) in the next time interval, / + At. This process leads to the exact solution for all time, given as 7 u(\,t) = H{ct) X + 2u < t i f t < f Am. _ t ^ Am_z2. < t < Am. t - Am -r- Am s t < Am+ 2 1 c c — c (2.20) , . 0 i f r < -TT or 4m - 1 < ? < 4w + 1 p(0,t) = <{ c -2 i f Am - 3 < t < Am - 1 (2.21) where w = 1,2,3, • • •. In Figure 2, the exact solutions to the model problem are shown when u = 1, p = 1, and A = 0. Figure 2 : Exact solutions (/i = 1, p = 1, A = 0) Various numerical methods such as the Forward and Backward Euler methods, Trapezoidal rule, Midpoint rule, and Simpson's rule are used to approximate the solution of this system. Also being tested on the model problem are the e and half-step schemes [11] which are both methods obtained by modifying the Trapezoidal rule to achieve enhanced stability characteristics. 8 Chapter 3 Discretization Schemes In the system of boundary integral equations (2.19), the integral ^p(0,T)dr appears. When solving for the exact solution, the Fundamental Theorem of Calculus is applied and this integral is dealt with easily. Numerical integration must be used when approximating the solution. There are many options available when choosing a method for numerical integration. The various schemes tested on the model problem are outlined in this chapter. In practice, the same shape functions and collocation points are used to approximate (2.5) in the general formulation, resulting in the same set of approximations that are obtained for the model problem. The model problem is simplified by equating the two equations of (2.19) evaluated at time t — -jr. Through this process, u is eliminated altogether and the system becomes 2[t-^H(ct-\)+H(ct-2)^ ~p(0,T)dr + H(ct) j " p(0,r)dr = 0. (3.1) If t > then (3.1) can be rewritten as $'op(0,x)dx + \'~Tp(0,r)dv = 2 (- J - - / ) . (3.2) For the remainder of this thesis, it is assumed that t > when dealing with the model problem. Let y = f(x) be a continuous function on some interval [a, b]. Consider a partition of [a,b] into n subintervals, each of width Ax = ^ ~ a . Let/j = / (x ( ) . For the model problem, assume that [0,7] is divided into A 7 intervals of size At = where T > -jr. Let /?, « p(0JAt), where i is a positive integer. To approximate j c p(0, r)dz, assume that A: is a positive integer chosen so that kAt < ± (k+l)At>^. (3.3) Define y as r = k~i5t W and let 9 t = T+At. 3.1 The Forward Euler Method (3-5) y=f(x) Figure 3 : The Forward Euler method In this case,/ is assumed to be constant on each subinterval. Sampled at the left endpoint of each subinterval, /r-i f f(x)dx * A x V / + 0(Ax) as depicted in Figure 3. (3.6) /V-l i=0 For the model problem, \Tp(0,r)dr * A / V V and from Figure 4, it can be determined that [ c p(0,r)dr « A H V + ypN-k-\ Then /?(0, T) « is given by the formula '"A'-A-1 (3.7) (3.8) 10 (3.9) Figure 4 : Approximation of j c p(Q, z)dr using the Forward Euler method 3.2 The Backward Euler Method For this method, / is approximated by a constant function on each subinterval and is collocated at the right endpoint of the subinterval, so b " f f(x)dx * Ax + O(Ax) (3.10) ,=i as illustrated in Figure 5. 11 y=f(x) Figure 5 : The Backward Euler method Figure 6 : Approximation of J c p(0, z)dz using the Backward Euler method Considering the model problem, \Tp(0,T)dT * A(Vp,-J o and from Figure 6, (3.11) N-k [ c p(0,z)dz « A M V p , - + 7/>/v-t 0 I M (3.12) 12 The value p(0, T) is approximated byp N . At time t, the function is approximated by p^+\ which is given by PN+\ = -N 1=1 N i f 0 < 7 < 4-i f - i - < 7 < - 2-1 1 c - c (3.13) ;=1 1=1 XPAf-t+i i f t > -jt 3.3 The Midpoint Rule Figure 7 : The Midpoint rule Once again, / is approximated by a constant function on each subinterval but here it is sampled at the midpoint m, of the subinterval, where rrij =f^ X>-X+Xi ^ _ Then b " n — ' | f(x)dx ~ Ax^m,• + 0(Ax2) = Ax^fi+± + 0(Ax2), (3.14) " (=1 i=0 as shown in Figure 7. 13 For the model problem, A M f p(0,r)dr « AtY\p, (3.15) and from Figure 8, T-2. (N-k-l | " P(0,T)dT ~ At< V p / +_L + 7 / ? ^ . (3.16) Ifp(o, 7 - -4/-') ~ PN—L> t n e n a t time T+ -4/-, is approximated by pN+±. This is given by the formula r PN+A N-l Y,Pi+Ar 1=0 I X + + £(7-4-) (=0 N-1 iV—t ,-=o " l ' ,=0 2 i f 0 < / < i f — < 7 < — 1 1 c — c YPN-k+i. *f * - c (3.17) 14 3.4 The Trapezoidal Rule On each subinterval, / is assumed to be linear with the Trapezoidal rule. Figure 9 shows that \bMdx = ^ { l > - y ( / o + / » ) } + 0(Ax2). (3.18) The standard basis functions for the Trapezoidal rule are depicted at the bottom of Figure 9. yf Figure 10 : Approximation of j c p(0,r)dr using the Trapezoidal rule 15 The integrals appearing in (3.2) are approximated as follows : r _ J L C N-k j ^ c p{0,r)dx ~ Af- < 2 ^ jTjPi-Po + YPN-k-\ + (y - \)PN-k (3.19) (3.20) Equation (3.20) is illustrated in Figure 10. The function p(0, T) for the model problem is approximated byPN, SOp(0,l) is approximated byp^+\ which is given by the formula r PN+\ = -<, 22> (=0 -Po i f 0 < / < - i -N /=0 i f T * 1 < T N -2po + ^ { l - i ] N-k+l 2zZPi i=0 ) + 2 X p , + YPN-k + (y - i)pN-k+\ 1=0 i f 7 > | -(3.21) 3.5 Simpson's Rule Assume there are an even number of subintervals, so n = 2m where m is a positive integer. In this case, / is approximated by a quadratic function on each subinterval. That is, 16 b f " a " / 2 _ 1 I \ f[x)dx * 2 J > , + 4 2 / ' + ' - (/"o + / » ) r + 0 ( A x 4 )- (3"22) 0 I i=0 1=0 J This method is shown in Figure 11. The basis functions for Simpson's rule are at the bottom of the picture. Assume there is an even number of subintervals, so N = 2m where m is a positive integer. Each subinterval of size 2At is approximated by a quadratic function. For the model problem, \Tp(0,T)dT * Jo 3 N/2 N/2-\ 2 2 ^ 2 / + 4 ^ />2H-I - (po + PAO i=0 1=0 (3.23) and from Figure 12, it can be determined that r l~Cp(0,T)dT~ ^ (N-k-2)/2 (N-k-2)/2-\ 2 ^ / ?2<+4 ^ P2M-P0 i=0 1=0 ^ (3.24) +(1 + y)pN-k-2 + (2 + 7)(4pyv_i_i +pN-k) where both A 7 and A: are positive even integers. The function />(0, T) is approximated by At time 7, « p(0,l) and is given by N/2 N/2-l 2 X)/>2i + 4 X) P2i+\ ~(po+ PN) i=0 1=0 N/2 N/2-l i f 0 < t < 1 2 $ > 2 , + 4 X P21+1 ~(PO+PN) + 4r(t - ^) i f - J - < f < i=0 i=0 A f where JV/2 M2-1 2HP2i + 4 X /72/+1 - 2^0 -pN-k-2 ~PN i=0 1=0 (N-k)/2 {N-k)/2-\ + l r ( 7 - i ) + 2 2>< + 4 S " i=0 1=0 +(2 + y){pN-k + 4/>AHt+l + 2/+1 I f t > f (3.25) 4- i f ? is odd 4 •y i f / is even • 1 i f i = 0,N (3.26) 17 Notice that (3.25) is dependent on N, unlike any of the other methods that have been looked at so far. Figure 12 : Approximation of j c p(0,T)dr using Simpson's rule 3.6 The £ -Scheme From previous work [11], the Trapezoidal rule is not always stable. However, using linear functions on each subinterval to approximate the actual function is desirable because then discontinuities introduced by using constant approximation schemes are avoided. Maintaining this linear quality, the e -scheme was designed as a perturbation of the Trapezoidal rule with improved stability characteristics [11]. The e -scheme is simply the Trapezoidal rule with more weight placed on the last unknown. The function/ is assumed to be linear on each subinterval, so \"f(x)dx = ^ j | > - y/o + jf*«\ + 0(A*2) where f„. 1 + 6 1 - e -/„. This is depicted in Figure 13. (3.27) Assume that l - e 1 + 6 Then (3.28) 18 yv-i f p(0,z)dr * A M VV- - \p0 + J o VPS 2 1 + 6 2 ( 1 - 6 ) PN (3.29) and from Figure 14, it can be determined that N-k-\ J O C p(0,T)dr * - y - J 2 />,• - / > 0 + XPJV-JW + (1 + r)jB/>AHt | . (3.30) At time 71, /?(0,7) « pn. At the next time-step,p(0, t) « /?AM which is given by the formula PN+\ = -B< I=0 1 e 2 g / ? / - J P o + T ^ 7 ^ + ^ - ( 7 - - J - ) i f 0 < 7 < - J -i f - i - < 7 < -2-1 1 c — c i=0 N-k i f 7 > +2 + TP N-k + ( l + y)PpN-k+i 1=0 (3.31) 3.7 The Half-Step Scheme y f Figure 15 : The half-step method {n even) To attempt to increase the accuracy and stability of the Trapezoidal rule, the half-step method was developed [11]. For the indirect boundary element method, the half-step scheme has been demonstrated to be more stable than other numerical methods. The solution is advanced in two half-steps of size -y- while the time-step for the convolution remains to be of size At. Since the convolution is the most expensive operation computationally, the half-step 20 method is a viable numerical scheme because the time-step used for the convolution is the same as that used for the other methods discussed in this chapter. Computer run times do increase by about 50%, but the slower running time is a sacrifice made in the hope that the half-step scheme is more accurate and more stable than the other methods. Similar to the Trapezoidal rule, the solution for the half-step scheme is assumed to be linear on each subinterval. Assume that there are n subintervals, each of size -4^ -. Then r nil J f(x)dx = A x J X hi - 9"(/o + /«) i f " is even («-D/2 H hi - 4-/o - \ U - \ + -jh i f n is odd i=0 ^ + 0(Ax2) (3.32) as shown in Figure 15. Also shown in Figure 15 are the standard basis functions for the half-step scheme. - i — — i - \ — — i - - i — — i -T 4 H 1 1 1 • t = T N^-k-2^ N-k-l LN-k t 0 = 0 Figure 16 : Approximation of J c p(Q,z)dT using the half-step method (A7 and k even) Assume there are N subintervals, each of size -4/-. For the half-step scheme, the following integrals are defined as j o 2 p(0,r)dT * At-l P 2 i - ±.(po+Ptf-\)> i f A^is odd, (3.33) J p(0,T)dr * At4^P2i--j(po+PN) M f A ' i s e v e n . (3.34) i=0 21 The value {N - k) must be even to approximate the integral J c p(0, T)dr, so there are two possibilities to consider : 1) A 7 and k are both even or 2) A 7 and k are both odd. In both cases, k is a positive integer chosen using the following criteria : JcM. < 2. k 2 - c (k + 2)-f > (3.35) Let 7 = r+4L. (3.36) Then (N-k-2)/2 c p(0, x)dx ~ " j " S 4 ^ P2i-po + (2y- k)pN-k-2 + (1 + 2 / - A)/?*-* 1=0 as shown in Figure 16 and for the next time-step t, r {N-k-iyi 4 ^ P2i -2po + (2y -k+l)pN-k-2 1=0 + (2/ - £ + 3 )/>#_* (i\Mt)/2 4 ^ /72i - 2/70 + (2y - k - 1 )pN.k ;=o + (2/ - k+ \)pN-k+2 (3.37) i f {k+l)-^f < J -i f ( * + l ) ^ - > (3.38) The function/>(0, J) is approximated bypn and at the next time-step,p(0, t) « • For the case where (£ + l)-4p < /7/v+i is given by the formula 22 (AM)/2 2 E P2i~P0 1=0 2 E ^ - / > o + £ ( / - - £ - ) (=0 " f (N-W2 (N-k-2)/2 2 P2i-2po + - ^ ( t - ^ ) + 2 2 P2i '=° (=0 {(2/ -k+1 )pN-k-2 + (2y-k + 3)pN-k} i f A 7 is odd and i f A 7 is even,pN+\ is given by r PN+\ = -< N/2 ^LP2i - 2p0 -pN (=0 N/2 4Y,P2i-2Po-PN+4r(l-f) i=0 " f N/2 {N-k-2)/2 i f 0 < t < -jr i f - i - < 7 < -2-1 1 c - c i f 7 > (3.39) 4TP2i-4po-pN+ -£r(t - ± ) +4 ^21 <=° 1=0 i f 0 < t < - J -i f < 7 < • 2-1 1 c - c i f 7 > f +(2y - * + 1 )p w-i_ 2 + (2 / - k + 3)p N-k (3.40) If ( k + l)-4f- > the formulas for/?jv+i remain unchanged for t < -2-. For / > -2-, the formula becomes {N-k)/2 PN+\ = -< 2 E ^ 2 / - 2 p 0 + -^ -(/-4") +2 2 ^2< <=° i=0 {(2/ - A: - 1 )/?Ar^t + (2y - k + 1 )pN-k+2} N/2 ( N ~ k V 2 4TP2i-4po-pN + -£r(j-Y)+4 X) P2i < = 0 1=0 +(2/ - k - 1 + (2/ - A: + l)pN-ic+2 i f A 7 is odd i f A 7 is even (3.41) 23 Chapter 4 Stability of the Model Problem 4.1 Stability Analysis Using the z -transform To determine the stability properties of the various discretization schemes applied to the model problem, the z -transform is used. As the discrete analogue of the Laplace transform, the z -transform is used for analyzing discrete time systems and is a technique for solving difference equations with constant coefficients. The z -transform of a sequence of numbers {f„} is defined as oo Ztf„} = F(z) = X > K " (4-1) «=o while if j{t) is a continuous function with the sampling period At, the z -transform is defined to be co Z{f{kAt)} =F(z) = Y 1, then an exponential numerical instability results and the numerical approximation is unstable. If there is a pole of the function T(z)F(z) of multiplicity greater than one occurring on the unit disk |z| = 1, then a resonant numerical instability may occur and the numerical method is called marginally stable. If all of the poles of T(z) fall inside the unit disk |z| < 1, then the discretization scheme is stable. The roots of the transfer function are sometimes difficult to determine analytically. In these cases, M A T L A B was used to calculate the roots for various values of k. Once known, the roots can be plotted to easily display where they are located with respect to the unit circle. For each approximation method, there is a unique transfer function which characterizes the stability of the method. The transfer functions for the methods outlined in Chapter 3 wi l l be derived in the remainder of this chapter. Assume that p(0, T) is approximated by px, so the interval [0,7] is divided into A 7 intervals of size At = where T > -7-. Let pt = p(0,iAt), where / is a positive integer. Assume that A; is a positive integer chosen such that kAt < 1 (k+\)At>^. (4.7) For the remainder of the chapter, define y as 4.2 The Forward Euler Method B y approximating (4.3) at time Tand at the following time-step T+ At, two equations are formed : 25 {A M N-k-\ ~\ + + YPN-k-x \ = 2(4- - tN) = 2Atf-^-N) i=0 ;=o J L i'=0 i=0 J Subtracting (4.9) from (4.10) results in p N + ( \ + y)pN-k - YPN-k-\ = -2. Taking the z -transform of (4.11), the following is obtained : This can be rewritten as P(z) = z— . ' ^ 2 z _ z A + 1 +(1 + y)z-y z - 1 ' so the transfer function for the Forward Euler method is ~k+l TF{z) = zk+[ + (1 + y)z - y Notice that the transfer function is dependent on k and independent of N. (4.9) (4.10) (4.11) (4.12) (4.13) Figure 17 : Poles of 7>(z) (At = 0.06 => k = 23, y * 0) Now that the transfer function is known, it is necessary to find the poles of 7>(z). This 26 can be broken down into two cases : k = (7 = 0) and k * ( 7 * 0 ) . If 7 i= 0, then all of the poles of 7>(z) fall inside the unit circle implying that the Forward Euler method is stable (see Figure 17). If 7 = 0 => At = where A: is a positive integer, then (4.14) reduces to z* + 1 In this case, all k poles of 7>(z) fall on the unit circle so the Forward Euler method is marginally stable. 4.3 The Backward Euler Method The two equations found by approximating Equation (4.3) at time Tand T+ At are f N N-k ~\ ti\^Pi + Y^Pi + YPN-k> = l ( ± - t N ) (4.16) fN+\ N-k+\ When (4.16) is subtracted from (4.17), the resulting difference equation is pN+\ + (1 + y)pN-k+\ - ypN-k = - 2 . (4.18) Taking the z -transform of equation (4.18), P(z) is found to be p ( z ) = - H i ^ (4-19) w zk+i +(1 + 7 ) z - 7 z - 1 v ' and the transfer function for the Backward Euler method is z* + l + (1 + y)z - y TBV = . „ Z • (4-20) This is similar to the transfer function for the Forward Euler method and in fact, the denominators of 7>(z) and TB(Z) are identical. Since the transfer function for the Backward Euler method has the same poles as the transfer function for the Forward Euler method, these two approximation schemes have the same stability properties. If 7 * 0, all of the poles of TB(Z) fall inside the unit circle so the Backward Euler method is stable. When 7 = 0, the poles of TB(z) all fall on the unit circle so in this case the Backward Euler method is classified as marginally stable. 27 4.4 The Midpoint Rule Equation (4.3) approximated by the Midpoint rule at time Tand T+ At becomes C N-\ N-k-\ AU ^ p i + x + X) Ptr-L + rPN-k-± \ =2(T-^) (4-21) . /=0 <=0 N N-k At-l J > i + x + ! > < H + rPN-k+± l = 2 ( i - tN+l ) . (4.22) I i=0 i=0 J Subtracting (4.21) from (4.22), the result is pN+± + (1 + Y)pN-k+± - 7PN-k-± = -2- (4.23) The Midpoint rule approximates p(0, ti) by a constant function on each subinterval similar to the Forward and Backward Euler methods, so the results are closely related to the previous two methods. The equation P(z) = —= ^ . (4.24) w zk+i +(1 +y)z-y z- 1 v ' is found by taking the z -transform of Equation (4.23) and the transfer function for the Midpoint rule is = M ^ • <4'25) zk+' + (1 + y)z - y The poles of TM(z) are the same as the poles of 7>(z) and 7g(z) since the denominators of these three functions are identical. A l l of the poles of TM{Z) are contained in the unit circle i f 7 * 0 , implying that the Midpoint rule is stable. In the special case where 7 = 0, the Midpoint rule is marginally stable since the poles of TM(Z) all fall on the unit circle. Since the Midpoint rule is more accurate than the previous two methods but shares the same stability characteristics, it is the most desirable scheme thus far. 4.5 The Trapezoidal Rule At time Tand at the next time step T+ At, Equation (4.3) is approximated by the following equations : 28 j N N-k 1 yM 2 + 2 2^'' ~ 2 / ? 0 + + - -PN > =2f^-tN} L /=0 1=0 J (4.26) r M - l N-k+\ i •y-S 2 ^ p , + 2 ^ /?/ - 2/?0 + 7/?^-* + (7 - 1 )/>#-*+1 -/'A'+I f = - tN+l ). I ;=0 i=0 J (4.27) B y subtracting (4.26) from (4.27) the resulting equation is PN+\ +PN + (1 + j)pN-k+\ + PN-k - 7PN-k-\ = -4 . (4.28) The function ( Z ) = z i + 2 + z t + 1 + ( l + 7 ) z 2 + z - 7 " ( 4 - 2 9 ) is found by taking the z -transform of Equation (4.28), thus the transfer function for the Trapezoidal rule is T ^ Z ) = — 2 M f l " ' ^ 2 • (4'3°) Figure 18 : Poles of 7>(z) (At = 0.06 => A = 23, 7 * 0) When 7 * 0, z = -1 is a pole of 7>(z). This one pole falls on the unit circle while the rest of the poles fall inside the unit circle. Here, the Trapezoidal rule is marginally stable since the function F(z) = ~ 4 z also has a pole falling on the unit circle. If 7 = 0, then the transfer 29 function reduces to 7Y(z) Tk+\ + zk + z + i (4.31) The (k + 1) poles of this modified transfer function (4.31) all fall on the unit circle. Therefore, the Trapezoidal rule is marginally stable for any choice of At. Although the trapezoidal scheme has the same order of accuracy as the midpoint scheme, it has less desirable stability characteristics. 4.6 Simpson's Rule Approximating (4.3) at time Tand T+ 2At, the following equations are obtained : N/2 N/2-\ (N-k-2)/2 (N-k-2)/2-\ ALJ 2 2^ 2'+4 2 P I M + 2 X P 2 I + 4 2 PIM r <=0 ;=0 i=0 (=0 -2p0 + (1 + Y)PN-k-2 + (2 + j)(ApN-k-\ + PN-k) -PN (N+2)/2 (N+2)/2-\ (N-k)/2 {N-k)/2-\ AL) 2 2 P2'+4 2 PIM+2 2 P2T+4 X P2M (=0 1=0 1=0 1=0 -2p0 + (1 + r)PN-k + (2 + y)(4pN-k+\ +pN-k+2) -PN+2 r = 2 ( ± - t N ) (4.32) (4.33) The following difference equation results from subtracting (4.32) from (4.33) : PN+2 + 4pN+\ +pN + (2 + y)iPN-k+2 + 4pN-k+\) + pN-k - (1 + y)(4pN-k-i + PN-k-2) = -12. (4.34) The function P(z) is found to be P(z) ,k+2 . -12z zk+4 + 4ZM + zk+2 + ( 2 + 7 ) Z 4 + 4(2 + 7 ) z 3 + z 2 - 4 ( 1 + 7 ) 2 - ( 1 + 7) z-l (4.35) by taking the z -transform of Equation (4.34), so the transfer function for Simpson's rule is Tk+2 Ts(z) = zk+4 + 4zk+3 + zk+2 + (2 + yy> + 4(2 + r ) z i + z2 - 4(1 + r)z - (1 + 7) 3 , ,2 (4.36) Even though it is more accurate, Simpson's rule is classified as unstable for any At since there are some poles of the transfer function Ts(z) that fall outside of the unit circle when 7 * 0 30 and when 7 = 0 (see Figure 19). 1.5 0.5 -0.5 -1.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Figure 19 : Poles of Ts(z) (At = 0.05 => k = 28, 7 * 0) 4.7 The £ -Scheme The € -scheme has been developed [11] to improve the stability characteristics of indirect boundary element approximations to the elastodynamic equations. In this section we analyse the stability properties of the e -scheme when applied to the direct boundary element model problem. Using the e -scheme, Equation (4.3) at time Tand T+ At becomes f" A M N-k-\ "\ J^Pi + 2 2 P>- 2P° + rPN-k-\ + (1 + Y)PPN-k + PPN \ = 2(-J- - tN) 1=0 1=0 N-k (4.37) 1 2 : + 2 2^ ~ 2 p o + 7 P N - k + (1 + riBPN~M + PP™ f = 2 (c" ~ / j v + 1 ) • I i=0 (=0 J (4.38) Equation (4.37) is subtracted from (4.38) to obtain P>AM + (2 - P)pN + (1 + y)PpN-k+\ + [ 2 - / 3 + 7 ( 1 - p)]pN-k - ypN-k-x = - 4 . (4.39) Taking the z -transform of (4.39), P(z) is found to be 31 P(z\ = - • ——- (4 40) K J Bzk+2 + (2-P)zk+]+(l+y)pz2 + [2-B + y(l-B)]z-y z - 1 v ' ' hence the transfer function for the e -scheme is j (z) = z (4 4\) e w Pzk+2 + (2-B)zk+]+(\+y)Pz2 + [2-B + y(l-6)]z-y' v ' 7 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Figure 20 : Poles of Te(z) (e = 0.2, At = 0.05 => k = 28,y * 0) - i Figure 21 : Poles of Te(z) (e = 0.35, At = 0.05 => k = 28 , / * 0) 32 For any value of y * 0, all of the poles of Te(z) fall inside the unit circle and the e -scheme is said to be stable. If y = 0, then the transfer function reduces to p V + 1 + (2 - B)zk + Bz + (2-B) (4.42) There are some poles of T€ (z) that fall on the unit circle implying that the e -scheme is marginally stable in this case. Figures 20 and 21 show the poles of Te(z) when At = 0.05 for € = 0.2, 0.35. 4.8 The Half-Step Scheme First, consider when A 7 and k are both odd. Equation (4.3) at time 7 is At {N-\)I2 (N-k-2)/2 4 2 P 2 i + 4 2 P 2 i ~ 4P° + ~ k)PN~k-2 i=0 i=0 +(2y -k + 2)pN-k -pN-\ + PN J (4.43) Now there are two cases to examine : (A: + l ) - y - < or + 1) > If (*+ 1 ) ^ - < then at time T+ (4.3) becomes r At (N+\)/2 (N-k-2)/2 4 2 P2i + 4 2 P2i ~ 4P° + (2T-k + 1)/>AHM 1=0 i=0 +(2y -k + 2)pN-k - 2pN+\ and the difference equation obtained by subtracting (4.43) from (4.44) is 2/7/V+l ~PN+ PN-l + PN-k + PN-k-2 = -4 . Taking the z -transform of (4.45), P(z) is found to be P(z) = ,k+2 2zk+3 _ zk+2 +zk+\ + z 2 + l z _ 1 (4.44) (4.45) (4.46) and the transfer function is 33 1 2zk+3 - zk+2 + zk+] + z2 + 1 (4.47) In the second case where (k+ l)-4r > 7 k Equation (4.3) becomes A/ < {N+\)/2 {N-k)l2 4 X] + 4 _ ) ^ 2/ - 4/»0 + (2y - * - l ^ w (=0 1=0 +(2y - A: + l)/?jv-i+2 - 2pN+\ 1 = 2 ( - J - - ^ i ) (4.48) at time T + -y- . Subtracting (4.43) from (4.48) results in the following difference equation : 2pN+\ ~PN+ PN-\ + (2/ - k + \)pN-k+2 + PN-k - (2y - k)pN-k-2 = -4 . (4.49) Taking the z -transform of (4.49), . -4z P(z) = Tk+2 2z*+3 _ z*+2 + zk+\ + (2/ - A: + 1 )z 4 + z 2 - (2 / - k) z - 1 resulting in the transfer function 7W2 _ 7k+2 2z* + 3 - zk+2 + zM + (2y - k + l ) z 4 + z 2 - (2/ - A:) ' (4.50) (4.51) Consider the case where N and k are both even. Equation (4.3) approximated at time T is r At N/2 {N-k-2)/2 4 _ ^ 2 , ' + 4 2 P2i - 4po + (2y - k)pN-k-2 1=0 1=0 +(2/ -k + 2)pN-k - 2pN l = 2 ( - i - - ^ ) . (4.52) Once again, this can be broken down into two cases : (A: + l)-4p < or (A: + < If (* + 1)4*- < - 2 r , (4.3) becomes M2 {N-k-2)/2 AtJ 4 _)^2' + 4 X ^ 2 / - 4/>o + (2y - A; + 1 )/7^-t-2 i=0 i=0 +(2y - A: + 3 )/?#-* + PN+\ l = 2 ^ - t N + i ) (4.53) 34 at time T+ -4p and the difference equation found by subtracting (4.52) from (4.53) is PN+I +PN+ PN-k + PN-k-i = -4 . (4.54) The z -transform of (4.54) is taken to find that „k+2 A„ and the transfer function is -k+2 TH, = -n—h—5 • ( 4 - 5 6 ) 3 zk+3 + zk+2+z2 + l V ' l .5r -0.5 -1.5 -1.5 -1 -0.5 0.5 1 1.5 Figure 22 : Poles of THL (At = 0.06 => k = 23, y * 0) \f(k+ 1 ) ^ - > | - , Equation (4.3) becomes r NI2 (N-k)/2 At J 4 ^ 4^P2i + 4 2 pv -4p0 + (2y-k-\)pN-k i=0 i=0 I =2(4--/AW). (4.57) +(2/ - A + \)pN-k+2 -PN+ PN+\ Equation (4.52) is subtracted from (4.57) to obtain the difference equation PN+I + PN + (2y - k + l)pN-k+2 + PN-k - (2y - k)pN-k-2 = -4 . Taking the z -transform of (4.58), the result is (4.58) 35 P(z) = —^ . ——- (4 59) v ; zk^+zk+2 + (2y-k+l)z4+z2-(2y-k) z-l K J y ) so the transfer function is zk+2 Tfi4 = zk+i + zk+2 + (2y-k+iy+z2- (2y - k) ' ( 4 " 6 0 ) For any value of y, there are poles that fall outside of the unit circle for each transfer function TH{, i = 1,2,3,4. This means that the half-step scheme is unstable for any choice of At. The roots of THL are plotted in Figure 22 for At = 0.06. 4.9 Summary of Theoretical Stability Properties In some cases, the stability of the time-stepping schemes is dependent on y = k -where k is a positive integer. A stable approximation method for values of / * 0 is only marginally stable i f y = 0 => kAt = This value is significant because after the initial 1 1 time interval of size the function p(0,t) is constant over time intervals of size It is interesting to note that the stability is dependent on the value k = -^j, so the stability is dependent on the chosen value of At. The analytic results summarized in Table 1 wi l l be compared with the experimental results discussed in the Chapter 5. Numerical Scheme Stability y * 0 y = 0 Forward Euler Method Stable Marginally Stable Backward Euler Method Stable Marginally Stable Midpoint Rule Stable Marginally Stable Trapezoidal Rule Marginally Stable Marginally Stable Simpson's Rule Unstable Unstable e -Scheme Stable Marginally Stable Half-Step Method Unstable Unstable Table 1 : Analytic stability classification 36 Chapter 5 Numerical Results 5.1 Intermittent Instabilities At Stable? 100 JL 50 Jl_ 46 K 100 J2_ 38 ft_ 30 0.02 = 0.0271 ••• = 0.0282-'•• 0.03 = 0.0307-•• = 0.0314--0.035 = 0.0372-•• 0.04 = 0.0471---0.05 yes yes no yes no yes yes no yes no yes Table 2 : Evidence of intermittent instabilities for the Forward Euler method Intermittent instabilities in boundary element approximations of the elastodynamic equations have been reported in the literature [1, 11]. Decreasing the time-step does not ensure that an instability wi l l be removed. In fact, i f a numerical approximation is stable for a particular At, it may become unstable as At increases and become stable again as At continues to increase. This type of intermittent instability is evident when numerical approximation methods are used to solve the boundary integral formulation of the elastodynamic equations. Stability analysis enables one to predict when an instability wi l l occur. In this chapter, numerical evidence is provided to illustrate the validity of the stability analysis that was performed in Chapter 4. 37 When the direct boundary integral method is employed, problems occur regardless of the numerical scheme being used when At is chosen such that kAt = -jr where k is a positive integer. For example, i f « = 1, p = 1, and X = 0, then c = Jl. As shown in Table 2, an intermittent instability is evident for a standard constant approximation method such as the Forward Euler method. 5.2 Numerical Experiments A l l of the programs used to test the various discretization schemes were written in M A T L A B . The equations necessary to approximate p(0,t) are " ( V - i ) = j-f^[(t-±)H(ct-l)+H(ct)l'Qp(P,T)dz] (5.1) u(\,t) = ^ tH(ct)+H(ct- ' p(Q,r)dx For the remainder of the chapter, let pt = p(0, iAt) where i is any non-negative integer and define where A: is a positive integer chosen such that kAt < £, (k+\)At>^. (5.3) Assume that \i = 1, p = 1, and X = 0, hence c = 41.. The time interval [0,15] is divided into N subintervals each of width At. For each approximation method, the numerical solution is plotted with the exact solution for two different values of At. In the first figure of each section, At = 0.01. This graph is representative of results obtained when choosing any At such that kAt + (y * 0). For the second figure of each section At = ^2c ~ ^-00995---. Typical results when the time-step is chosen such that kAt = -7- (y = 0) are illustrated in this figure. A value for k can be found from the conditions given by (5.3). For At = 0.01, k = 141 => 7 * 0 and for A/ = ^ , k = 142 => y = 0. 38 5.3 The Forward Euler Method -0.5 © " 1 -1.5 - exact solution numerical solution 10 - J — • 15 Figure 23 : Approximation using the Forward Euler method (At = 0.01 => y * 0) o _ Figure 24 : Approximation using the Forward Euler method {At 142c 7 = 0 ) 39 Using the Forward Euler method,PN is given by Equation (3.9). The numerical approximation forp(0,t) is plotted with the exact solution for At = 0.01 and for a slightly smaller At 142c = 0.00995- • • in Figures 23 and 24 respectively. The analytic results in Chapter 4 predict that the numerical solution is stable for At = 0.01 and marginally stable for At = 142c ' Coinciding with the analytic results, the Forward Euler method approximates the solution accurately when At - 0.01 and when At = 142c the approximation goes unstable. 5.4 The Backward Euler Method -0.5 o -1 -1-5F - e x a c t s o l u t i o n n u m e r i c a l s o l u t i o n 10 15 Figure 25 : Approximation using the Backward Euler method (At = 0.01 => y * 0) Givenp^, the approximation at the next time-step is defined by (3.13). In Figures 25 and 26, the numerical approximation is plotted against the exact solution for At = 0.01 and At = -r- 2— = 0.00995-••. In Section 4.3, the analytic results predict that the approximation is stable for At = 0.01 40 and marginally stable for At = . Figure 25 shows that the Backward Euler method is stable for At = 0.01 and Figure 26 shows that it is unstable for At = . ~. . 6 142c 5.5 The Midpoint Rule For the Midpoint rule approximation, pN+± is given by Equation (3.17). Figures 27 and 28 show the numerical approximation plotted against the exact solution for At = 0.01 and At = -r^- = 0.00995-••. The approximation is stable for At = 0.01 as predicted by the analytic results. When At = £ , the Midpoint rule approximation is predicted to be marginally stable and from Figure 28, it is evident that there is an instability present. Since the Midpoint rule approximation behaves in a different manner from the Forward and Backward Euler methods when 7 = 0, Figure 29 shows the Midpoint rule with At = ^ applied over a longer time period. It is observed that although the Midpoint rule does perform better than the Forward and Backward Euler methods when y = 0, it is still not stable and is not a good numerical approximation to use in this case. 41 exact solution numerical solution -0.5 o -1 a . -1.5 k -2 10 15 Figure 27 : Approximation using the Midpoint rule (At = 0.01 => y * -0.5 S " 1 -1.5 - exact solution numerical solution 10 15 Figure 28 : Approximation using the Midpoint mle (At 142c 7 = 42 o -1 in J I ! Ll : L - exact solution numerical solution Hi 10 15 20 25 30 35 40 45 50 Figure 29 : Approximation using the Midpoint rule (At = -^lc r = o, T = 5.6 The Trapezoidal Rule -0.5 a -1.5 -2 -2.5 • exact solution numerical solution 10 15 Figure 30 : Approximation using the Trapezoidal rule (At = 0.01 => y * 0) 43 The formula necessary to approximatepN +\ givenPN is defined by (3.21). The numerical approximation is plotted along with the exact solution for At = 0.01 and At = = 0.00995-•• in Figures 30 and 31 respectively. From Section 4.5, the analytic results predict that the Trapezoidal rule is marginally stable for any value of At. The numerical approximation is unstable in both Figures 30 and 31. 5.7 Simpson's Rule Using Simpson's rule, the approximation forPN+\ is defined by (3.25). Figures 32 and 33 show the numerical approximation using Simpson's rule plotted against the exact solution for At = 0.01 and A ; = - i - = 0.00995-•-. Simpson's rule is predicted to be unstable for any value of At in Section 4.6 and the numerical results provide evidence to support that prediction. 44 -5 exact solution numerical solution 5.8 The £ -Scheme 0.5 -0.5 -1.5 -2 -2.5 - exact solution numerical solution 10 — ' • 15 Figure 34 : Approximation using the e - scheme (e = 0.2, At = 0.01 => y * 0) -0.5 h S "1 -1.5 k • exact solution numerical solution 10 15 Figure 35 : Approximation using the e - scheme (e = 0.35, At = 0.01 => y 0) 46 o . exact solution numerical solution Figure 36 : Approximation using the e - scheme (e = 0.2, At 142c 7 = 0 ) o a Figure 37 : Approximation using the e - scheme (e = 0.35, At = 142c 7 = 0 ) Using the € -scheme,PN+I is given by Equation (3.31). The exact solution and 47 numerical approximation for = 0.01 and At = ^ - 0.00995- • • with e = 0.2 are plotted in Figures 34 and 36. These figures are compared with Figures 35 and 37 which are plotted for e = 0.35. It is apparent that the e -scheme is more stable in the case where e = 0.35. The analytic results show that the stability roots for e = 0.2 and e = 0.35 were fairly similar, with a couple of roots located closer to the center of the unit circle for e = 0.35. The corresponding solution in which the trapezoidal scheme is used is shown in Figures 30 and 31. In this case, the instability of the trapezoidal scheme (e = 0) can be clearly observed. For At = 0.01, the e -scheme is stable as predicted in Section 4.7. From the analytic results, the e -scheme is expected to be marginally stable for At = an<^ n u m e " c a ' approximation appears to be unstable in Figures 35 and 37. Further numerical experiments reveal that the value of the parameter e is an important factor in determining the behaviour of the approximation. The e -scheme is merely a perturbation of the Trapezoidal rule. In fact, i f e = 0 then the € -scheme reduces to the Trapezoidal rule which was determined to be unstable for all values of At in Section 5.6. If e > 1, the solution explodes exponentially and the e -scheme is not useful. The range of values for which the e -scheme is valid for is 0 < e < 1. Even within this range, the numerical solution behaves differently depending on what value of e is chosen. As e increases, stability of the approximation increases while accuracy decreases. If e is chosen to be in the interval [0.5, 1), the solution dies out as time increases. Optimal results or the best combination of stability and accuracy are achieved by choosing e such that 0.3 < e < 0.45. 5.9 The Half-Step Scheme For the half-step scheme, i f (k + 1)-^- < -jr, then the valuePN+\ givenPN is defined by Equation (3.39) i f N is odd and (3.40) i f N is even. If(k+l)4f > TT,PN+I is given by (3.41). In Figures 38 and 39, the exact solution and numerical approximation using the half-step scheme are plotted for At = 0.01 and At = -r^— = 0.00995-••. The analytic results predict the half-step approximation to be unstable. From Figures 38 and 39, it is apparent that this prediction is correct. The half-step approximation results in an exponential instability for any choice of At. This result is consistent with numerical evidence which indicates that the half-step scheme is unstable for the direct boundary element method [3]. 48 exact solution numerical solution H -9x 10 o _ j I -18x 10 Figure 38 : Approximation using the half-step scheme (At = 0.01 => y * 12 10 A 6 o 4 _ - exact solution numerical solution 16x 10 9x 10 c 10 15 Figure 39 : Approximation using the half-step scheme {At = => y = 49 5.10 Summary of Numerical Stability Results The numerical results provide evidence to support the idea that the stability of the discretization scheme is dependent on the value of y = k - ~^r- No approximation method is stable for this model problem when y = 0. Table 3 summarizes the experimental results. Numerical Scheme Stability y * 0 y = 0 Forward Euler Method Stable Unstable Backward Euler Method Stable Unstable Midpoint Rule Stable Unstable Trapezoidal Rule Unstable Unstable Simpson's Rule Unstable Unstable e -Scheme Stable Unstable Half-Step Method Unstable Unstable Table 3 : Experimental stability classification The results of the numerical experiments thus agree with the theoretical predictions of the stability analysis. B y using the z -transform, the stability of various discretization schemes applied to the model problem can be predicted accurately. To ensure a stable numerical approximation for the model problem, a constant approximation method or the 6 -scheme should be implemented and A/ should be chosen such that kAt * -fr. 50 Chapter 6 Conclusions The importance of the elastodynamic equations in the geosciences has provoked interest in using boundary element methods to obtain approximate solutions of the elastodynamic equations for problems with general geometries or fault and fracture interactions. Boundary element methods are an alternative to finite element and finite difference methods which have traditionally been used to solve these problems. Finite element and finite difference methods require that the entire domain be discretized and they also suffer from problems with numerical dispersion, so it would be advantageous to gain a better understanding of boundary element methods which could be useful when dealing with problems of this type. Although computational time is saved, there is evidence of intermittent instabilities when boundary element methods are applied. Even for the one-dimensional elastodynamic 9 model problem described in Chapter 2, instabilities appear when kAt = where k is any positive integer. A l l of the discretization schemes are unstable when kAt = The reason why this is the case may be that the reflected wave is causing a resonance, resulting in an instability. Analytically, it is easy to see why an instability occurs here. For any discretization scheme, at least one pole of the transfer function lands on the unit disk causing that numerical scheme to be marginally stable i f At is chosen such that kAt = -jr- Experimentally, all of the schemes classified as marginally stable in Section 4.9 have instabilities. This means that in practice, they would be considered unstable. The z -transform is used to perform a stability analysis on various discretizations of the model problem. A n equation of the form P(z) = T[z)F(z) is obtained, where T(z) is the transfer function that determines the stability of the given discretization scheme. The results obtained from the stability analysis match up well with what is obtained experimentally. This implies that the z -transform is a valuable tool in determining the stability of a numerical method and predicting when instabilities wi l l occur. When kAt * the Forward and Backward Euler methods, Midpoint rule, and the e -scheme are stable while the Trapezoidal rule, Simpson's rule, and half-step scheme are unstable. As stated before, the Trapezoidal rule and e -scheme are closely related. For kAt * the role of the parameter e is to pull the one pole of the transfer function Tj{z) that is on the unit disk inside the unit disk, stabilizing the Trapezoidal rule and making the e -scheme a useful method for values of e such that 0 < e < 1. The e -scheme improves the performance of the Trapezoidal scheme, which agrees with numerical evidence [3]. The half-step scheme is not a viable approximation method for this problem because it is exponentially unstable and the solution explodes very quickly. For 51 the indirect boundary element method, the half-step scheme had superior stability properties but here it was found to be unstable for the direct boundary element method. This is consistent with numerical results obtained in [3]. Knowing that the occurrence of an instability can be predicted accurately for this one-dimensional model problem, it is likely that the z -transform can be applied to more complicated problems and give satisfying results. Understanding what causes the instabilities and developing criteria to ensure that a numerical solution is stable when applied to a particular problem are the first steps in allowing boundary element methods to become a powerful tool in approximating the solutions of the elastodynamic equations. 52 Bibliography [1] D.J. Andrews, "Dynamic growth of mixed-mode shear cracks", Bull. Seism. Soc. Am., vol. 84, 1184-1198,(1994). [2] P .K. Banerjee and R. Butterfield, Boundary element methods in engineering science, McGraw-Hil l , Maidenhead, (1981). [3] Birgisson, A .P . Peirce and E . Siebrits, "Elastodynamic direct boundary element methods with enhanced numerical stability properties", Int. J. Num. Meth. Engng., vol. 46, 871-888, (1999). [4] J. Dominguez, Boundary elements in dynamics, Computational Mechanics Publications, Southampton, (1993). [5] C H . Edwards and D.E. Penney, Calculus with analytic geometry : early transcendentals, Prentice-Hall, New Jersey, (1998). [6] A . Frangi and G. Novati, "On the numerical stability of time-domain elastodynamic analyses by B E M " , Comp. Meth. Appl. Mech. Engng., vol. 173(3-4), 403-417, (1999). [7] A . J . Jerri, Integral and discrete transforms with applications and error analysis, Marcel Dekker, New York, (1992). [8] P. Lancaster and K . Salkauskas, Transform methods in applied mathematics : an introduction, John Wiley & Sons, New York, (1996). [9] G.D. Manolis and D.E. Beskos, Boundary element methods in elastodynamics, Unwin Hyman, London, (1988). [10] G.D. Manolis and T.G. Davies, Boundary element techniques in geomechanics, Computational Mechanics Publications, Southampton, (1993). [11] A .P . Peirce and E. Siebrits, The stability properties of time domain elastodynamic boundary element methods, unpublished monograph, (1995). 53