A MULTIGRID M E T H O D OF F O R D E T E R M I N I N G LITHOSPHERIC T H E D E F L E C T I O N PLATES By Paul M . Carter B . Sc., T h e U n i v e r s i t y of Sheffield, 1985 A T H E S I S 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 O F T H E REQUIREMENTS FOR T H ED 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 O F G R A D U A T E STUDIES INSTITUTE O F APPLIED MATHEMATICS W e - accept this thesis as c o n f o r m i n g t o the required s t a n d a r d T H E UNIVERSITY O F BRITISH C O L U M B I A J u n e 1988 © P a u l M . C a r t e r , 1988 In presenting degree this thesis in partial fulfilment at the University of British Columbia, freely available for reference and study. copying of this department or thesis by for scholarly his or publication of this thesis MATHEMATICS DE-6 (2788) 24 th JUNE 1988 I further agree that permission for extensive purposes may be granted It is by the head understood that for financial gain shall not be allowed without The University of British Columbia Vancouver, Canada Date for an advanced I agree that the Library shall make it her representatives. permission. Department of of the requirements of my copying or my written Abstract Various models are currently i n existence for determining the deflection of lithospheric plates under an applied transverse load. T h e most popular models treat lithospheric plates as t h i n elastic or t h i n viscoelastic plates. The equations governing the deflection of such plates have been solved successfully i n two dimensions using integral transform techniques. Three dimensional models have been solved using Fourier Series expansions assuming a sinusoidal variation for the load and deflection. In the engineering context, the finite element technique has also been employed. T h e current a i m , however, is to develop an efficient solver for the three dimensional elastic and viscoelastic problems using finite difference techniques. A variety of loading functions may therefore be considered w i t h m i n i m u m work involved i n obtaining a solution for different forcing functions once the m a i n program has been developed. T h e proposed method would therefore provide a valuable technique for assessing new models for the loading of lithospheric plates as well as a useful educational tool for use i n geophysics laboratories. The m u l t i g r i d method, which has proved to be a fast, efficient solver for elliptic partial differential equations, is examined as the basis for a solver of b o t h the elastic and viscoelastic problems. T h e viscoelastic problem, being explicitly time-dependent, is the more challenging of the two and will receive particular attention. M u l t i g r i d proves to be a very effective method applicable to the solution of b o t h the elastic and viscoelastic problems. n Table of Contents Abstract ii List of Tables v List of Figures vii Acknowledgement viii 1 Introduction 1 2 Derivation of the Governing Equation 9 2.1 Lithospheric Plates 9 2.2 Classical E l a s t i c i t y and Viscoelastic Theory 10 2.3 T h i n E l a s t i c Plates 13 2.4 Transverse L o a d i n g 19 2.5 Nondimensionalization 20 2.6 B o u n d a r y and Initial Conditions 21 2.7 Reformulation of the governing equations 22 3 The Multigrid Method 3.1 25 T h e basic method 25 3.1.1 Relaxation 28 3.1.2 Restriction 31 3.1.3 Prolongation 31 iii 3.1.4 4 5 6 The Coarse G r i d Operator 32 3.2 Basic M u l t i g r i d A l g o r i t h m s 32 3.3 M u l t i g r i d for coupled equations 35 3.4 M u l t i g r i d for time-dependent problems 37 The Elastic Problem 41 4.1 A n a l y t i c Solution 41 4.2 N u m e r i c a l Solution 46 4.3 Numerical Results 48 4.3.1 Periodic L o a d i n g 48 4.3.2 M u l t i g r i d vs. Gauss-Seidel 52 4.3.3 Gaussian L o a d i n g 53 The Viscoelastic P r o b l e m 56 5.1 A n a l y t i c Solution 56 5.2 N u m e r i c a l Solution 58 5.3 N u m e r i c a l Results 61 5.3.1 Periodic L o a d i n g 61 5.3.2 Integration w i t h steady forcing 69 5.3.3 Geophysical A p p l i c a t i o n 71 Conclusion 74 6.1 Results 74 6.2 Future W o r k 76 Bibliography 79 iv List of Tables I Degree of support (Elastic Model) 45 II l'/D' as a function of plate length 47 III L = 2.5 x 1 0 m , Periodic Loading 50 IV L = 5.0 x 1 0 m , Periodic Loading 50 V L = 1.0 x 1 0 m , Periodic Loading 6 50 VI M u l t i g r i d vs. Gauss-Seidel 52 VII L = 2.5 x 10 m , Gaussian Loading 54 VIII L — 5.0 x 10 m , Gaussian L o a d i n g 54 IX L = 1.0 x 10 m , Gaussian L o a d i n g 54 X Degree of support (Viscoelastic M o d e l ) 58 XI B a c k w a r d Euler, L M G vs. M F M G 63 XII C . D . F . method, L M G vs. M F M G 63 XIII Two-step B . D . F . , L M G vs. M F M G 63 XIV B a c k w a r d Euler, XV C . D . F . method, XVI Two-step B . D . F . , XVII B a c k w a r d Euler, At = 5 M a 67 XVIII C . D . F . method, A< = 5 M a 67 XIX Two-step B . D . F . , A t = 5 M a 67 XX Steady state integration : case 1 70 XXI Steady state integration : case 2 70 5 5 5 5 6 At = lOMa At = lOMa A t = lOMa v -. . 65 ". . 65 65 XXII T i m e dependence i n terms of m XXIII M a x i m u m deflection (metres) for case (i) XXIV M a x i m u m deflection (metres) for case (ii) vi List of Figures 2.1 Elements of Stress 11 2.2 T h i n plate deformation 14 2.3 H y d r o s t a t i c Restoring Force 19 3.4 T h e m u l t i g r i d two-level cycle 27 4.5 Degree of support as a function of nondimensional topographic wavelength 44 vii Acknowledgement I would like to acknowledge the enthusiasm, encouragement and advice offered by my supervisors D r . U r i Ascher and D r . M a t t Y e d l i n during the research and w r i t i n g of this thesis. vm Chapter 1 Introduction T h e a i m of this dissertation is to provide a fast, efficient, numerical method for determ i n i n g the deflection of lithospheric plates under an applied load. In particular, the m u l t i g r i d method is employed i n this capacity. We w i l l examine plates having horizontal dimensions ranging from those of sedimentary basins to continental plates. We w i l l also consider a wide range of other physical parameters. Therefore, our atten- tion is focussed on developing a method capable of providing a solution to a variety of geophysical problems rather than on developing a good model to one particular geophysical situation. In other words, we shall concentrate m a i n l y on the numerical aspects of the topic while still maintaining a geophysical framework w i t h i n which our forcing functions and model parameters w i l l lie. Our discussion will begin w i t h a justification of the treatment of lithospheric plates as t h i n elastic or viscoelastic plates. In elastic theory one assumes that at any given time the plate responds only to the applied load at that time so that once the load is removed, the plate instantaneously reverts to its undisturbed state. In viscoelastic theory the assumption is that the plate responds not only to the present load but also to all previous loads so that it has the capacity to remember previous states of deformation. O n removal of the load, a viscoelastic plate slowly relaxes to its undisturbed state. We thus find ourselves in the regime of t h i n plate deflection to which much attention has been given in the engineering literature. T h e equations governing t h i n elastic or viscoelastic plate deformation have already been solved numerically, quite successfully, 1 Chapter 1. Introduction 2 using finite element packages, but such methods tend to be very expensive. In a few special cases, an analytic solution may also be obtained allowing for a more objective evaluation of the numerical results. T h e application of the current work to the evolution of sedimentary basins is of particular interest since such areas tend to be rich i n hydrocarbons. Sedimentary basins are generally considered to be-formed by the combination of a load due to thermal contraction following a sub-surface heating event and a subsequent load due to infilling of the deflection by sedimentary layers. It is the load due to thermal contraction which is generally considered to be responsible for the i n i t i a l deflection in the plate and incidentally, the thermal anomaly causing this contraction also plays an important role in the m a t u r a t i o n of petroleum, see N u n n , Sleep and M o o r e [8]. Various models for the evolution of sedimentary basins have been developed for b o t h elastic and viscoelastic plates. Beaumont [2] examined the formation of sedimentary basins due to the loading of elastic and viscoelastic lithospheres. T h e problem was considered i n cylindrical co-ordinates w i t h the deflection depending only on the radial component. A simple ^-function was used to represent the spatial variation (point loading) in the forcing term for the elastic model. For the viscoelastic problem, the forcing function was augmented w i t h the heavy-side-step function to represent the time dependence. M o d e l l i n g of the plate loading was therefore very simple, however, mathematical tractability was maintained i n the process. B o t h the elastic and viscoelastic problems were solved using integral transform techniques. Beaumont also addressed the problem of basin i n i t i a t i o n , i.e. that mechanism by which an i n i t i a l deflection is created in the lithosphere i n which sediment subsequently accumulates. Six possible mechanisms were proposed i n c l u d i n g thermal contraction following a sub-surface heating event, a subsiding graben and necking due to stretching of the lithosphere. Deflection of the lithosphere due to the sedimentary infill of an Chapter 1. 3 Introduction instantaneously formed graben was examined. T h e case of an exponentially subsiding graben was also considered where subsidence was assumed to occur over a period of lOOMa. Results from the viscoelastic, rather than elastic, model were found to be i n best agreement w i t h seismic data obtained from the N o r t h Sea basin. W a t t s , K a r n e r & Steckler [12] examined a two-dimensional model where thermal contraction was assumed to control the tectonic subsidence of the plate, i.e. that subsidence which is independent of subsequent deflection due to sedimentary infill. Fourier transform techniques were used to solve the governing equations which are simplified by the assumption of two-dimensionality. Sedimentary loading was assumed to be constant over discrete time intervals so that at the beginning of each new period a layer of sediment was assumed to instantaneously infill the cavity caused by tectonic subsidence. T w o models were considered for tectonic subsidence the latter being attributed to thermal contraction i n b o t h cases. T h e first assumed that the depth of the basin increases as the square root of time since basin i n i t i a t i o n . A g a i n , formation of the basin was assumed to occur over a p e r i o d of lOOMa. T h e second model was based on the assumption that the lithosphere is stretched laterally at the time of basin i n i t i a t i o n resulting i n plate t h i n n i n g and heating which in t u r n r e s u l t / i n plate subsidence on cooling. A variation i n elastic thickness, and therefore flexural rigidity, w i t h time was allowed for i n b o t h the elastic and viscoelastic. models between time periods. Spatial variation of the flexural rigidity was also considered. T h e latter problem was solved using a numerical finite difference scheme rather than Fourier Transform techniques. It should be noted that this problem is simplified considerably i n two-dimensions, the extension to three dimensions being non-trivial. It was concluded that the best model fit to overall basin geometry arises from an elastic model i n which the flexural rigidity increases with Chapter 1. Introduction 4 time after basin i n i t i a t i o n . T h i s increase is due to an increase i n elastic thickness w i t h age as the plate cools. It should be noted that the conclusion here is at variance w i t h Beaumont who favoured the viscoelastic model although the latter model assumed the flexural rigidity to be constant. It was again deduced that the dominant mechanisms i n the formation of sedimentary basins are thermal contraction, which controls the overall shape of the basin, and sedimentary loading. N u n n & Sleep [7] considered the load due to b o t h thermal contraction and sedimentary loading for an elastic and viscoelastic lithosphere. A sinusoidal spatial variation was assumed for b o t h and solutions were obtained using Fourier Series expansions. A linear deposition of sedimentary layers over discrete time intervals was considered where a variation in sediment density from one layer to the next was incorporated i n the model. A n estimate of the thickness of sedimentary layers obtained from lithostratigraphic data for the M i c h i g a n Basin was used to compute the Fourier coefficients for the expansion representing the load due to sedimentary infill. Fourier coefficients for the load due to thermal contraction were computed using the fact that the deflections due to thermal and sedimentary loading must add to yield the observed deflection at some time, t , after basin i n i t i a t i o n . These coefficients were therefore re-computed for max each set of rheological parameters. It was estimated that the load due to sedimentary infill constitutes approximately 75% of the total driving force. T h e best model fit to observations was obtained for a viscoelastic lithosphere w i t h low flexural rigidity. However, it was pointed out that w i t h the d a t a available at present, either model yields acceptable results, and a better understanding of basin i n i t i a t i o n , sediment budget and sediment compaction is necessary before any definitive statements can be made. We w i l l examine the deflection due to the loading of b o t h elastic and viscoelastic Chapter 1. 5 Introduction plates, all models being three dimensional. Simple models which account only for sedimentary loading w i l l be considered. We therefore assume that some process such as thermal contraction is responsible for the formation of a cavity in which the sediment accumulates. In the viscoelastic case, we allow the basin to form gradually so that the sediment accumulates i n a continuous fashion. The infilling sediment is assumed to have constant density. We now present an overview of material to be found i n the following chapters. In C h a p t e r 2 we w i l l call u p o n material found i n the engineering literature to provide a derivation of the governing equations for the deflection of a viscoelastic plate. In doing so we w i l l closely follow a derivation of the equation governing deflection of an elastic plate replacing the elastic constitutive equations w i t h those for a viscoelastic solid. W e w i l l also present appropriate boundary and initial conditions along w i t h a nondimensionalized form of the governing equations. To conclude the chapter we w i l l consider an appropriate reformulation of the governing equations which w i l l make their solution easier numerically. For b o t h the elastic and viscoelastic problems, we are faced w i t h solving a pair of coupled, linear, elliptic equations. A n efficient solver for such a system is therefore required. In Chapter 3 we present an introduction to the multigrid method and consider its application to the two problems i n hand. The presentation includes a discussion of the m a i n components of m u l t i g r i d as well as some basic algorithms. T h e fundamental difference between the elastic, and viscoelastic problems, from a numerical standpoint, is that the viscoelastic model is time-dependent whereas the elastic problem is not. C h a p t e r 3 is therefore concluded w i t h a discussion of m u l t i g r i d for time-dependent problems. T w o algorithms are considered, one of which follows immediately from the time-independent, elastic problem. T h e applicability of the two algorithms is discussed and some attempt is made to ascertain the conditions under w h i c h one is expected Chapter 1. 6 Introduction to perform more efficiently than the other. T h e performance of the two algorithms is tested numerically in Chapter 5. In C h a p t e r 4 we examine the elastic problem on a square domain. T h e assumption of a square domain places some restriction on the problem geophysically, but it still maintains a definite advantage over the assumption of two-dimensionality which has been used by other authors. However, the regular d o m a i n leads to an easier treatment of the problem from a numerical standpoint since irregular domains involve added computational difficulties near the boundaries and hinder the development of fast, efficient solvers. We begin w i t h the presentation of an analytic solution to which later numerical results w i l l be compared. T h i s analytic solution is used to examine the effect of altering the elastic properties of the plate on subsequent, plate deflection. A relation expressing the degree to which a load is supported by the flexural rigidity of the plate will be derived. It w i l l be shown that the problem is potentially tougher to solve in cases where the degree of support is either very high or very low. The rate of convergence and accuracy of the m u l t i g r i d solution w i l l be examined for various plate dimensions and elastic constants representing a wide spectrum of degrees of support. A comparison between m u l t i g r i d and a simple Gauss-Seidel iteration will be carried out showing that m u l t i g r i d performs relatively more efficiently except in cases where the degree of support is very low i n which case the numerical solution becomes t r i v i a l . We conclude the chapter by considering a Gaussian distribution for the loading for which an analytic solution is not available. In Chapter 5 we turn to the viscoelastic (time dependent) problem which is certainly the more interesting of the two problems in hand. We begin w i t h the presentation of an analytic solution representing an exponentially increasing forcing function w i t h periodic spatial dependence. T h i s representation has some major shortcomings i n the Chapter 1. Introduction 7 geophysical context but it is nonetheless considered acceptable i n so far as it provides a means for assessing the numerical results more objectively. T w o m u l t i g r i d algorithms for time-dependent problems w i l l be contrasted and compared. After an examination of the results and a consideration of the work involved i n each of the algorithms, one w i l l be eliminated on the grounds of b o t h efficiency and accuracy. Some analysis of the reason for the efficiency and accuracy of one algorithm over the other w i l l also be presented. H a v i n g chosen a suitable algorithm we then examine the performance of three different approaches to dealing w i t h the time discretization, namely, B a c k w a r d E u l e r , a centred difference scheme and a two-step backward difference formula. T h e first of these provides a first order a p p r o x i m a t i o n i n time to the solution while the other two provide second order approximations. It w i l l be seen that the centred difference scheme, while providing a second order approximation, involves no more work or storage than Backward Euler, the first order method, m a k i n g it an extremely attractive method from the point of view of b o t h accuracy and efficiency. We examine the effect of reducing the time step and, in consideration of the results, the centred difference scheme is selected for use i n future numerical experiments. Next we consider the effect of allowing the forcing to reach a steady state. A n a l y t ically we expect all plates, regardless of size or elastic properties, to relax u n t i l they are i n hydrostatic e q u i l i b r i u m w i t h the applied load. We conduct two experiments for plates of different sizes and different elastic properties and show that if integration is continued for long enough after the forcing reaches a steady state then the numerical solution indeed appears to converge to the solution representing an adjustment of the plate to hydrostatic e q u i l i b r i u m . Finally, we consider a forcing function which is more reasonable geophysically. T h e spatial variation is represented by a Gaussian distribution and the time dependence is Chapter 1. Introduction 8 such that forcing is initially zero and reaches a steady state asymptotically. T w o sizes of plate w i t h different elastic properties w i l l be considered and again the numerical results indicate that if integration is continued for a sufficiently long period, the solution converges to one representing hydrostatic e q u i l i b r i u m . Conclusions arising from the numerical experiments are presented i n Chapter 6. Suggestions as to the direction of future research are also proposed. Chapter 2 Derivation of the Governing 2.1 Lithospheric Equation Plates In considering the deflection of lithospheric plates, we will be encountering horizontal dimensions ranging from 250km to 1000km. T h i s covers anything from the formation of sedimentary basins to the deflection of continental plates due to topographic loading. In order to derive an equation governing plate deflection we make the following observations. F i r s t , that the horizontal extent of the plate is not so large as to make the effects of the earth's curvature of importance. W e therefore assume that i n their undisturbed state, the plates are completely flat. Second, the thickness of the plate is very small i n comparison to the horizontal dimension. In general, lithospheric plates have an elastic or mechanical thickness of between 25 and 50km depending on whether the plate is oceanic or continental. T h i s is considerably less than the seismic thickness which is of the order of 100km, see N u n n & Sleep [7, p.590]. It is the mechanical thickness w h i c h is of importance i n an elastic or viscoelastic model so we therefore assume that lithospheric plates may be treated as thin elastic or viscoelastic plates. In elastic theory, we assume that at any given time the plate responds only to the applied load at that time so that on removal of the load, the plate instantaneously returns to its undisturbed state, having no memory of the previously applied load. Should we wish to examine the deflection of the lithosphere under an applied load at a given point in time, the present day for example, then elastic theory provides us 9 Chapter 2. Derivation of the Governing 10 Equation w i t h an adequate model for achieving this. However, if we wish to study the of the plate's deflection then a viscoelastic model is more appropriate. evolution In viscoelastic theory the elastic properties are assumed to be a function of time so that materials have the capacity to remember, and therefore respond to, previously applied loads as well as the current one. If the load is removed at a given time, a viscoelastic plate does not instantaneously revert to its undisturbed state but rather relaxes, possibly to a steady state, as time goes on. Obviously, from a geophysical standpoint, we w i l l not be considering the complete removal of the load but, i n the case where the load reaches a steady state, the subsequent evolution of the plate's deflection is of interest. 2.2 Classical Elasticity and Viscoelastic Theory H a v i n g justified the treatment of lithospheric plates as t h i n elastic or viscoelastic plates we now derive an expression governing the deflection of t h i n viscoelastic plates under the action of transverse loading. We begin with an overview of classical elasticity introducing the ideas of stress, strain and subsequent particle displacement. Consider a volume element of an elastic material as shown i n Figure 2.1 w i t h Cartesian axes-in the given orientation. This volume element is subject to normal and shear stresses due to its interaction w i t h neighbouring volume elements. Figure 2.1 indicates the direction i n w h i c h these stresses are taken as positive. Note that stresses acting on surfaces normal to the x - d i r e c t i o n 2 have been omitted for clarity. Subscripts on the stresses are in the following convention: (i) the first represents the direction normal to the surface on which the stress is acting, (ii) the second indicates the direction parallel to which the stress is acting. W e adopt a suffix notation w i t h respect to Cartesian axes (x , x , x ) and assume the E i n s t e i n 1 2 3 S u m m a t i o n Convention to be i n operation over repeated suffices. T h e notation (. )_j Chapter 2. Derivation of the Governing 11 Equation ^33 Figure 2.1: Elements of Stress represents differentiation w i t h respect to x r Following V i n s o n [11], we now present the governing equations for the deformation of an elastic, body. A more in-depth discussion than that given i n V i n s o n may be found in Hunter [6]. Consideration of conservation of m o m e n t u m for a static, elastic body i n the absence of b o d y forces (such as gravity) yields the following e q u i l i b r i u m equations °~H,j = (2-2.1) 0 Conservation of angular m o m e n t u m requires that a, = on (2.2.2) ? in other words that the stress tensor be symmetric. A deformed elastic body has associated w i t h it not only stresses but also strains and a displacement field. A s s u m i n g small displacements, it can be shown that the strains e, (suffix convention being the same as for the stresses a ) are related to the 7 t] Chapter 2. Derivation displacement field of the Governing 12 Equation i n the following way H = \( *J + 3,i) £ u (2.2.3) U F i n a l l y we have the constitutive or strain-stress relations. A s s u m i n g the material to be isothermal, isotropic and homogeneous we have: £ij 1 (1 + v)uij - v u E kk 8 { (2.2.4) where the superscript e refers to an elastic solid, E = Young's M o d u l u s , v = Poisson's R a t i o and 1 6{j = Kronecker D e l t a •0 if % = j, otherwise. Drucker [3, p.267] shows that the strain-stress relations for a linear elastic solid lead directly to those for a linear viscous solid w i t h the replacement of the strains e\- by the rate of strain z - where the dot denotes differentiation w i t h respect to time and the v superscript v refers to a viscous solid. For an isotropic material we replace Young's M o d u l u s E, by the coefficient of viscosity, p. Hence, for a linear viscous solid we have: £ii = ~~ K + ) A 1 4 v ~ v t>ij} (2.2.5) where the assumption of isotropy and homogeneity has been maintained. The linear viscoelastic model consists of a combination of the above linear and viscous models. T w o possible formulations exist, the M a x w e l l model and the K e l v i n model. T h e M a x w e l l model assumes that the strain for a viscoelastic body is given by the sum of the strains arising from the elastic and viscous models. T h e K e l v i n model, however, assumes that the stress for a viscoelastic body is given by the sum of the stresses arising from these models. A s pointed out by Turcotte [9, p.70], the usual approach, at least in the geophysical context, is to deal w i t h the problem using the M a x w e l l model. Chapter 2. Derivation of the Governing 13 Equation We therefore assume that the rate of strain, e^, is given by hence, using equations (2.2.4) and (2.2.5) we have: £ij = (1 + z/) (7,;, U. E &kk Pkk ~E /T (2.2.6) As w i l l be seen shortly, for our purposes it will be of use to have a form of the constitutive equations where summation from only 1 to 2 and not to 3 is implied over repeated suffices. Inverting equation (2.2.6) w i t h this in m i n d we have: <7„ + o~ij T _ E (1 + 1/ ^ + (1 - ") [2.2.7) - kk^ij where T _ (2.2.8) E is k n o w n as the relaxation constant. 2.3 T h i n Elastic Plates V i n s o n [ l l ] derives an equation for the deflection of a t h i n elastic plate. We shall follow along the same lines using the constitutive equations for a viscoelastic solid i n place of those for an elastic solid, hence deriving an equation governing the deflection of a t h i n viscoelastic. plate. Consider a t h i n , rectangular, isotropic plate w i t h horizontal dimensions x\, x\ and thickness h. We assume that the conditions h -C x\ and h <C x% define a t h i n plate. We now position a set of Cartesian axes at the centroid of the plate as shown i n Figure 2.2 so that the mid-plane of the plate is the plane x = 0. 3 Consider a cuboid element of the plate extending through the plate thickness as shown. W e assume that on the application of a transverse load this element remains Chapter 2. Derivation of the Governing Equation 14 CUBOID E L E M E N T MID-PLANE Figure 2.2: T h i n plate deformation normal to the deformed mid-plane, i n other words, it undergoes at most a translation and rotation w i t h respect to the co-ordinate system. T h e element therefore remains straight upon load application. T h i s means that sides of the cuboid initially l y i n g i n the X ] a n d x -directions remain perpendicular to sides l y i n g i n the x - d i r e c t i o n after 2 deformation. 3 H u n t e r [6, p . I l l ] interprets the strains e tJ (i ^ j) as representing a measure of the deviation from perpendicularity, after deformation, of sides initially pointing i n the i a n d j - d i r e c t i o n s . Since i n the current problem sides i n the X j a n d x -directions remain perpendicular to sides i n the x - d i r e c t i o n , we m a y deduce that 2 e 13 3 — £ 23 = 0- W e also assume that an applied transverse load results i n plate bending but does not lead to compression or tension i n the plate thickness direction. T h i s is i n contrast to the behaviour of a sponge where a transverse load is absorbed by compression i n the thickness direction and not by bending. A g a i n referring to Hunter [6, p . I l l ] , we see that the strains £,-j (i — j) may be interpreted as representing some measure of Chapter 2. Derivation of the Governing 15 Equation the elongation or compression of an element in the z-direction. Since i n the current application there is no elongation or compression i n the ^ - d i r e c t i o n we therefore deduce that e 33 = 0. E q u a t i o n (2.2.3) then implies that u (the displacement i n the 3 x3 direction) is independent of x . 3 Hence, the admissible form of the subsequent displacement field for a t h i n elastic plate is: u = u (x ,x ,x ,t) = u t(x ,x ,t) + u = u (x ,x ,x ,t) = u° (x ,x ,t) + x a (x x ,t) x 2 1 2 1 2 1 2 3 3 u 2 — 1 2 2 1 = 3 where u°(xi,x ,t) < 2 x a (x ,x ,t) 3 3 2 u 2 1 1 2 > (2.3.1) u (x ,x ,f) 3 1 2 u , ( x i , x , 0, t) for i = l , 2 and a (xi,x ,t) 2 1 2 represent rotations to be defined later. F r o m here we shall adopt modified forms of the suffix notation and s u m m a t i o n convention which have been used to this point. W e shall assume that a l l suffices run only from 1 to 2 and not to 3. A n y variable w i t h a suffix 3 will be w r i t t e n explicitly. We now derive expressions for the net moments and forces acting on an element across the plate thickness. T h e moments, per unit area about the mid-plane are given by: = cr jX % 3 Integrating these across the plate thickness we have: h — J\ Mij a zj 3 x dx 3 (2.3.2) The transverse shears, q , i n the x - d i r e c t i o n are given by: l 3 and so integrating over the plate thickness we have: h Q % = f\ a l3 dx 3 (2.3.3) Chapter 2. Derivation of the Governing 16 Equation We now derive an integrated form of the equilibrium equations (2.2.1). M u l t i p l y i n g (2.2.1) by x and integrating over the plate thickness we have 3 0 - ^ 2 : 3 dx 3 + M ijtj + cr ,3X / l3 3 3 o OX l3 =*> = 0 dx %3 3 dx = 0 3 Mijj - Q , . = 0 assuming there are no shear stresses acting on the surface of the plate. (2.3.4) Integrating (2.2.1) for t = 3 we have 12 2 h o 3h] dx 3 + 0-33,3 dx I Qj,i 3 + = 0 o- 33 Qj,j +p = 0 =4- (2.3.5) where p = p(x ,x ,t) 1 2 = o- 33 represents the net force per unit area acting n o r m a l to, the plate surface. Combining equations (2.3.4) and (2.3.5) we have Mijj + p = 0 ti We now relate the moments M l3 (2.3.6) to the vertical displacement of the plate. C o m b i n i n g equations (2.2.3) and (2.2.7) we have Chapter 2. Derivation of the Governing 17 Equation where the zero subscript has been dropped from E and r . M u l t i p l y i n g through by x , 3 integrating and m a k i n g use of (2.3.1) we have: h 2=3 2 E / • 0 , - 0 \ , r dx-. 2 c dx-- (2.3.7) u (x ,x ,t) which •0 V [l + u) + 1 E .2 (1 + We now relate the rotations from here w i l l be denoted by that e 13 = e 23 a -[Ct h l 2 w = (1-^) to the vertical deflection a (xi,x ,t) % ) X3 + w(x , 1 2 Equations (2.2.3) and the assumption x ,t). 1 3 2 = 0 imply: 0 dxi £23 = dw dx J - ' dw =±> 9xi 0 dw dx a, 2 2 Notice we have dropped the c o m m a notation for derivatives and w i l l continue to do so for the remainder of this discussion. E q u a t i o n (2.3.7) may now be written dw dw + v dxjdx., dxf.dx 2 -D Mi. 1- V 2 (2.3.8) ZJ k where D , the flexural rigidity of the plate, is given by D = E h 3 12(1 - v ) 2 Equations (2.3.6) and (2.3.8) represent a closed system of equations for the vertical deflection w, i n terms of the applied transverse load p. E l i m i n a t i n g M l3 equations we have dw 2 dx\ I D dw D(l + 2 dxidx dx \ 2 | 2 2 dw 2 dx-t dx* from these Chapter 2. Derivation of the Governing dw dw .dx\ dx\ 2 dx\ 18 Equation 2 p , = - • +P T which may be written - (1 - u)0\D,w) V {D(V w)) 2 2 = V -+p (2.3.9) T where dx\ dx\ and d Dd 2 0 (D,w) 4 = dD : dw 2 w dx\ dx\ dx dx 1 2 dxidx 2 dD 2 + 2 dw 2 dx\ dx\ R e c a l l that we have assumed the plate to be homogeneous so that E and v are constant. If we further make the assumption that the plate thickness is uniform we may deduce that D is independent of position (although it may still m a i n t a i n some time dependence if h changes w i t h time) and so we have: DV w — P hp 4 (2.3.10) r wnere V\i; = V (V w) 2 2 Inverting equations (2.2.4) and using the resulting equations in place of (2.2.7) i n the above derivation leads to an equation relating the vertical deflection of an elastic, rather than viscoelastic, plate to transverse loading. V i n s o n [11, p.16] shows that the equation thus obtained is DV w=p 4 instead of (2.3.10). (2.3.11) Chapter 2. Derivation u of the Governing 19 Equation Surface Rock Pc Pa w h. •m Mantle B Fluid Mantle A Pm Figure 2.3: Hydrostatic Restoring Force 2.4 Transverse Loading W h e n considering the deflection of a lithospheric plate due to a transverse load, p, we must take into account the hydrostatic forces which result from such deflections. If a lithospheric plate having an overlying layer of sediment is deflected downward, the sediment which is usually less dense than the mantle, will infill the cavity thus formed, resulting i n a net hydrostatic restoring force. Following Turcotte & Schubert [10, p.121] we consider the case shown i n figure 2.3 where the continental lithosphere of density p m and thickness h m rides on the fluid mantle, also of density p . Overlying the lithosphere m is a layer of sediment of thickness h and density p . A n applied transverse load p c c a causes a deflection w i n the lithospheric plate. We shall now calculate the hydrostatic restoring force and hence the net transverse load. T h e pressure acting at the point A in F i g u r e 2.3 below the deflected plate is given by Chapter 2. Derivation of the Governing 20 Equation whereas the pressure at the point B under an u n d i s t u r b e d p o r t i o n of the plate is given by: p c g h + c p m g ( h m (2.4.2) + w) S u b t r a c t i n g (2.4.2) from (2.4.1) we have the hydrostatic force acting on the deflected plate: (Pc - Pm)gw (2.4.3) Hence, the net transverse load is given by P = where 7 = ( p — p ) g . m c p a + {p - c p m ) g W = p a - (2.4.4) JW Hence, the governing equation for deflection of a viscoelastic plate under an applied load p may now be w r i t t e n : a DV w A + 7 ^ - + iu) = — +p (2.4.5) a and for an elastic plate we have: DV w 4 + j w = p (2.4.6) a F r o m hereon, the subscript 'a' on the loading, p , will be dropped so that p w i l l be understood to mean the applied load. 2.5 Nondimensionalization For the purposes of numerical computation w e . w i l l now express equations (2.4.5) and (2.4.6) in nondimensional form. A s our solution domains will be square, we scale all lengths by a factor L . A l l times will be scaled by a factor T . U s i n g primes to denote nondimensional parameters, we have: w = w'L, p = p'/L, D = D'L 2 Chapter 2. Derivation of the Governing 7 = 7'/L , T — 2 d dx! Id d Ldx\\ 21 Equation T'T Id dx d i d Ldx' ' 2 dt 2 T dt' Hence, equation (2.2.5) becomes Cancelling the common factor 1/LT from each of the terms and dropping the primes we have: 4q - . ( , w , P -'\ DX7 w + 7 ^ - +wj= 7 (2.5.1) - + p where now, all variables and parameters are nondimensionalized w i t h respect to a length scale L and a time scale T . Similarly, equation (2.4.6) remains essentially unaltered i n nondimensional form. 2.6 B o u n d a r y and I n i t i a l C o n d i t i o n s Consider a solution domain Q, w i t h boundary d£l. Following N u n n & Sleep [7, p.597] we assume the plate to be simply supported at the edges so that there is zero deflection and zero curvature at the boundary, .w = 0, <9fi. Hence, V u> = 0 on dn 2 ' (2.6.1) For the time-dependent (viscoelastic) problem, we also need some i n i t i a l condition. A n y condition consistent w i t h (2.6.1) at t = 0 w i l l suffice. In cases where the exact solution, w , is known we w i l l take: e w = w {t = 0) on Q, e as our i n i t i a l condition where it is assumed that w e (2.6.2) satisfies (2.6.1). In cases where the exact solution is not known, we follow N u n n & Sleep [7, p.630] i n assuming that the Chapter 2. Derivation of the Governing 22 Equation deflection is initially zero so we have: w(t = 0) = 0 on n (2.6.3) A t this point we shall draw attention to an alternative formulation of (2.6.1) which w i l l be of use to us later. A s the boundary conditions (2.6.1) are assumed to hold for all time they may be w r i t t e n equivalently as: w = 0, Vw 2 = 0, on dft (2.6.4) Note that while the b.c.'s do not depend on time, the forcing function does. In the following chapter we shall consider possible solution methods for equations (2.4.5) and (2.4.6) w i t h the above boundary and i n i t i a l conditions. 2.7 Reformulation of the governing equations P a r t i c u l a r analytic solutions are readily available for certain choices of the forcing function. Examples of these will be encountered later where an analytic solution w i l l be used to assess the accuracy of a numerical solution. N u n n & Sleep [7] provide an in-depth discussion of the solution for periodic loading using the elastic and viscoelastic models. T h e current aim however, is to provide a means for solving the elastic or viscoelastic problem numerically i n as efficient a manner as possible. We propose to solve the equations using finite difference discretizations, however, as b o t h problems have 4 th order derivatives i n space, we are presented w i t h several computational difficulties if they are left in their current formulation. W e w i l l alleviate some of the difficulties by re-writing each equation as a pair of coupled, second order, p.d.e.'s as follows. E q u a t i o n Chapter 2. Derivation of the Governing 23 Equation (2.4.5) for the viscoelastic problem becomes: + 7 ^ — + uj = - + p DV v 2 Vu = v u = w 2 (2.7.1) w i t h boundary conditions: u = 0, v = 0 on dft (2.7.2) and initial condition: w(t = 0) = 0 on 9, (2.7.3) assuming the exact solution is unknown. E q u a t i o n (2.4.6) for the elastic problem may be written: D V u + jw = p 2 (2.7.4) V u> = v 2 w i t h boundary conditions: w = 0, v = 0 on dn (2.7.5) Notice that the boundary and initial conditions arise from a reformulation of those presented i n section 2.6. It is clear that the elastic, rather t h a n viscoelastic, problem is the easier of the two to solve. T h e solution of the elastic problem requires an efficient solver for coupled elliptic p.d.e.'s of the form: £ (u ,« ) = / 1 1 2 C (u\u ) 2 1 on = 0 2 n (2.7.6) where u = u = 0 on dn 1 2 (2.7.7) Later it will be seen that having chosen a discretization for the time derivative, the viscoelastic problem may, i n part, be written in the same form. So, once an efficient Chapter 2. Derivation of the Governing Equation 24 solver has been developed for the elastic problem, it may also be used to solve part of the viscoelastic problem which, on the grounds of efficiency and modularity, is desirable. Chapter 3 The Multigrid Method 3.1 T h e basic m e t h o d Consider the solution of Cu = f on ft (3.1.1) where C is an elliptic differential operator, and u — 0 on dtt (3.1.2) Define a grid £7^ on which a discrete approximation to u is to be computed. T h e n , u^, the exact solution of the discrete ^problem, satisfies: • Cu h h = f h on Q (3.1.3) and Uh = 0 on dn where C] is a. discretization of the differential operator C. and fh is a discrete approxx imation to / . T h e discrete problem (3.1.3) may be regarded as a large sparse system of equations for the solution values at grid points. Direct methods tend not to be efficient solvers for such problems owing to fill-in during the decomposition stage. T h i s fill-in is p a r t i c u l a r l y t a x i n g as regards storage ^requirements. More efficient methods, i n terms of storage, may be found amongst simple iterative schemes such as GaussSeidel. However, the asymptotic convergence of such methods tends to be slow and 25 Chapter 3. The Multigrid 26 Method for very large systems such methods are extremely slow and therefore expensive. T h e m u l t i g r i d method, about to be described, proves to be a fast, efficient solver for such a system. T h e m u l t i g r i d method can be viewed as a defect-correction method. Given an a p p r o x i m a t i o n Uh to Uh, the solution of (3.1.3), we define a defect, d^, by d = hh (3.1.5) £hU h and a correction Vh satisfying: C {U h h + v) = d + CU h h h h (3.1.6) which, i n the case of a linear operator Ch, reduces to = d Cv h h (3.1.7) h where v = 0 on dfl h (3.1.8) Hence, our original problem (3.1.3) can be w r i t t e n i n terms of a defect, dh. and a correction, Vh- H a v i n g computed dh and solved exactly for Vh, our corrected solution becomes u :=Ui h so that u h now satisfies (3.1.3) exactly. + v h (3.1.9) A s things are right now, we stand to gain nothing from our reformulation of the problem i n terms of :a defect-correction scheme since solving (3.1.6) or (3.1.7) does not appear easier than solving (3.1.3). However, as yet we have not introduced the idea of m u l t i g r i d d i n g . In the m u l t i g r i d method, we propose to construct a hierarchy of successively coarser grids on w h i c h the defect for the problem on the previous finer grid may be represented. For simplicity, i n this initial discussion we shall confine ourselves to only two grids, a fine grid 0 ^ and a coarse one £7 , where, for example, H = 2h and H C fit,,. Chapter 3. The Multigrid R E L A X Method 27 RELAX^ 1 C-h.Uk = fh 2 £h.Uh = fh # 4 p R E S T R I C T R O L O N G A T •E SOLVE C-HVH F i g u r e 3.4: T h e m u l t i g r i d two-level cycle W i t h a coarser grid defined, it is proposed to solve (3.1.6) on this coarser grid rather t h a n the finer one as this w i l l , of course, be c o m p u t a t i o n a l l y cheaper. However, the solution on QH cannot compensate for the high frequencyicomponents of the error which can be seen only on the fine grid. Hence, in order to achieve an approximate solution on the fine grid, the coarse grid correction or G G C is complemented by relaxation iterations on the fine grid. T h e elements of a m u l t i g r i d two-level cycle are presented i n Figure 3.4 where it has been assumed that is a linear operator. A s can be seen, the m e t h o d consists of four m a i n elements: 1. - Relax": generally consists of v sweeps of a simple iterative method such as GaussSeidel which serves to update the solution and to smooth out high frequency components of the error in the solution. It is applied before and after the coarse grid correction. 2. Restrict: computes the defect on the finer grid and transfers it to the coarse grid using an appropriate restriction operator 1^. 3. Solve: an exact, solver for the correction on the coarsest grid. Chapter 3. The Multigrid 28 Method 4. Prolongate: transfers the correction from the coarse grid to the fine grid using an appropriate prolongation operator The essential philosophy behind this two-grid method is that a simple relaxation scheme such as Gauss-Seidel can be a very efficient smoother even though it is slow to converge as a stand-alone iterative method. Before considering these components i n more detail we must decide on a grid coarsening process. We will restrict ourselves to square domains on which we have a uniform grid. A s s u m i n g that the finest grid contains an odd number of grid points i n each direction then a coarser grid can be defined by doubling the grid spacing i n each direction so that a l l points on the coarse grid are common to the fine grid and the grid remains uniform. M a n y other grid coarsening processes exist but the one described above is the easiest, to implement and will suit our purposes adequately. W i t h this assumed, we w i l l now consider the m a i n components of the m u l t i g r i d cycle in more detail. 3.1.1 Relaxation The purpose of the relaxation routine is to smooth out high frequency components of the error: (3.1.10) e = U -u h The h h generation of other error components is also undesirable and so the relaxation m e t h o d , i n itself, should be a convergent iterative scheme. We first describe some popular iterative schemes and their convergence properties. For ease of notation we define: A = C , x = u , b= f k h h so we are interested i n solving the system Ax — b. Let A be an approximation to A such that B = A~ x is easily computed. W e may now write an iterative scheme for Chapter 3. The Multigrid 29 Method solving the given system as follows: nO+i) The matrix M : = (I - BA)xW + Bb (3.1.11; := I — BA is called the amplification m a t r i x . The methods we are about to describe are determined by the choice of A. W e begin by considering a splitting of the m a t r i x A into three components: (3.1.12) A = D — E where D is diagonal, —E is strictly lower triangular and —F is strictly upper triangular. T w o of the most common iterative schemes of this type are: 1. The Jacobi Method where we choose A = D. We then have M = D' (E 1 + F) so that pointwise the iteration may b.e written: .(j'+i) 1 X] kl- Si) (3.1.13) a Q>kk 1=1 i^k 2. The Gauss-Seidel Method where we choose A D - E so that M = (D- E)~ F l and the pointwise iteration may be written: 1 o-kk Yl i<k a k i x i 3) (j+i) (3.1.141 l>k B o t h the Jacobi and Gauss-Seidel iterations can be shown to converge provided the m a t r i x A is strictly diagonally dominant. Gauss-Seidel, however, generally provides faster convergence than Jacobi. We now consider, the smoothing properties of these iterations but first let us consider why this should be of concern to us at a l l . The restriction operator, iff, may be viewed as an rnxn m a t r i x w i t h m < n. T h e kernel of this m a t r i x , k e r ( l j f ) , is therefore n o n - t r i v i a l . Hence, if dn •= iff dh represents Chapter 3. The Multigrid 30 Method the defect on the coarse grid, it is possible that a non-zero defect on the fine grid could give rise to a zero defect, djj, on the coarse grid. T h i s has disastrous consequences when solving the equation Cjjvjj = dn since if djj — 0 then VJJ = 0 and the iteration will converge to the wrong solution. It is therefore extremely important to determine those dh E ker(J^f) and find a suitable means of dealing w i t h them. A s observed by Wesseling [13, p.42] and as w i l l be seen here shortly, the restriction operator is a weighted average w i t h positive weights, and so any dh (E k e r ( i ^ ) must have frequent sign changes and therefore be highly oscillatory. Notice that the error, := Uh — u^., satisfies the relation: Ce h h = d h (3.1.15) Since C- arises from the discretization of a differential operator, it does not, as noted by h Wesseling [13, p.46], have a smoothing effect on e^,. Hence, high frequency components of the error, e^, w i l l give rise to h i g h frequency components of dh which as we have seen are undesirable. We therefore require some means of smoothing out high frequency components of ehIt can be shown that Gauss-Seidel iteration acts as an effective smoother for high frequency components of the error whereas the Jacobi method is not as effective i n this capacity. However, it should be rioted that variations of J a c o b i , i.e. damped J a c o b i , can be used as effective smoothers w i t h an appropriate choice of the d a m p i n g coefficient. The interested reader is referred to Hackbusch [5] for further discussion. Hackbusch also describes other possibilities which exist for performing the smoothing iteration such as Incomplete L U decomposition, Successive-Over-Relaxation, A l t e r n a t i n g Direction Implicit and Conjugate Gradient but we w i l l not discuss these methods here as GaussSeidel will prove to be adequate for our needs. Chapter 3. The Multigrid 3.1.2 31 Method Restriction T h e restriction procedure serves to transfer the defect from the fine to the coarser grid. T h e two most popular methods are: 1. Straight where the defect is injected through to the coarser grid at points Injection common to b o t h grids. 2. Full Weighting where a weighted average of defects at points neighbouring those common to both grids is transferred to the coarse grid. In stencil form this m e t h o d may be written: I H — 16 1 2 1 2 4 2 1 2 1 (3.1.16) In most circumstances, full weighting proves to be the better choice of the two since it uses information from all points on the fine grid rather than just those points common to b o t h grids. 3.1.3 Prolongation T h e purpose of the prolongation routine is to transfer the correction from the coarse grid to the finer one. T h i s may be achieved by injecting the corrections up to the finer grid at points c o m m o n to b o t h grids and then performing a bilinear interpolation. So, in stencil form we have: 1 4 1 2 1 2 4 2 1 -2 1 (3.1.17) Higher order interpolation may of course be carried out but, for our purposes, bilinear interpolation will be sufficient. Chapter 3. The Multigrid 3.1.4 32 Method T h e Coarse G r i d Operator F i n a l l y we must define the coarse grid operator CH- There are two methods available to us for achieving this. 1. W e may use the same discretization of the differential operator on the coarse grid as was used on the fine grid. So for example, if a 5-point scheme was used to discretize a Laplacian operator on the fine .grid, then the same 5-point scheme could be used on the coarse grid to define CH2. T h e G a l e r k i n method defines the coarse grid operator i n terms of the fine grid operator and the prolongation and restriction operators as follows: C H • •= I K C J H (3.1.18) T h e first method is the most popular for use w i t h finite difference discretizations owing to its shear simplicity. However, for finite element methods the second m e t h o d is the natural choice. It should also be noted that i n 2 dimensions, (3.1.18) gives rise to a 9-point stencil even if Ch is only a 5-point stencil. We shall be dealing only w i t h finite differences and so the first of these options w i l l be employed. 3.2 Basic M u l t i g r i d Algorithms E x t e n d i n g the concept of the two-grid method outlined i n section 3.1, we may now define a hierarchy of grids such that grid spacing on each coarser grid is double that of the preceding one. If our finest grid is square and contains 2 — 1 interior grid points n i n each direction then we could continue coarsening the grids until the coarsest grid contains only 1 interior grid point provided that the correct physical problem is still represented on the coarsest grid. We assume this to be the case so that the exact Chapter 3. The Multigrid 33 Method solution of the discrete problem on the coarsest grid is easily obtained by one sweep of Gauss-Seidel. E v e n though this is generally considered to be an iterative m e t h o d , when used i n this capacity w i t h only one grid point where the solution is u n k n o w n , it i n fact acts as an exact solver for the discrete problem. Should the coarsest grid contain more than one interior grid point, we simply postulate the existence of an efficient exact solver for the p r o b l e m on that grid. A s s u m i n g the'hierarchy of grids consists of p levels, I = p being the finest and / = 1 the coarsest, we have defects di, corrections vi and operators for each level. W i t h this n o t a t i o n we may now define a basic m u l t i g r i d algorithm, L M G , for the problem C u p = f p p as-follows: Algorithm 3.1 (LMG) Procedure L M G ( / , 7 , ui, u )\ 2 Begin if / = 1 then S O L V E ^ U i = e^) else begin for i:= 1 to u R E L A X ( £ ^ = d\); x RESTRICT; {d^ <- j / ( d j - 1 x £/Uj)} for j : = 1 to 7 L M G ( / - 1, 7, u , z PROLONGATE; {v t u )\ 2 <- v + t I^vi^} for i : = l to u R E L A X ( £ u , = di); 2 ; end; End; Notice that the recursive call to L M G is made 7 times at each of the coarser levels. Usually, on the grounds of efficiency, 7 takes on only the values 1 or 2. T h e case 7 = 1 Chapter 3. The Multigrid is k n o w n as the V(u u ) 1} 34 Method cycle while 7 = 2 is known as the W(ui,u ) 2 2 cycle where v x and v are the number of pre and post C G C relaxations. 2 A n alternative approach is to begin not on the finest level as i n L M G but on the coarsest level. T h i s method is commonly known as full m u l t i g r i d or F M G . H a v i n g obtained an approximation to the solution on the coarsest level this solution is interpolated, using a high order interpolation / / _ ! , to the next finer level. L M G is then applied at this level after which the solution is interpolated to the next finest level and so on. Hence, our F M G algorithm may be written: A l g o r i t h m 3.2 ( F M G ) Procedure F M G ( / , 7, v , u ); x 2 Begin SOLVE(£ u 1 1 = / ); a for I := 2 to p do begin {u <— f j ^ u ^ j } I N T E R P O L A T E ; t L M G ( / , 7, u u ); u 2 end; End; Finally, it should be noted that for nonlinear operators, (3.1.6) may have to be used i n place of (3.1.7) thus requiring that the solution-as. .well as the defect be transfered to each coarser level. T h i s method is known -as the full approximation scheme or F A S . However, we shall be dealing exclusively w i t h linear operators and therefore the. nonlinear case w i l l not be considered further. Chapter 3. The Multigrid 3.3 35 Method M u l t i g r i d for c o u p l e d equations We now consider the application of the multigrid method outlined above to the following problem. £V.« ) = / 2 £ (u\u ) 2 u 1 1 (3- - ) 3 (3.3.2) = 0 2 1 = u = 0 on dQ 2 (3.3.3) Recall that both the elastic and, i n part, the viscoelastic problem may be w r i t t e n i n this form. Discretizing these equations we may define a defect for each in the following way. O n the finest level we set d\ = f\ and d = 0 so that 2 h d\-.= dl-C\(ulu ) h (3.3.4) d :=d -C (ulu ) (3.3.5) 2 2 h 2 2 h 2 h h and a coupled system of equations for the corrections on n^ is also obtained as follows: 4(«) (3-3.6) = 4 (3.3.7) C l ( v l v ) = dl 2 h where vl = v = 0 on dn 2 h h (3.3.8) The transfer of defects and corrections between coarse and fine grids remains standard since i n our case the system (3.3.1) to (3.3.2) is only weakly coupled. However, the relaxation component requires some thought. A s we have a coupled system of equations there is some r o o m for manoeuvre i n this area. Basically we have to choose from the following: • Sequential Smoothing Chapter 3. The Multigrid 36 Method • Collective Smoothing We now consider each of these i n more detail. Sequential is performed by successively smoothing one component of the smoothing solution while keeping the other fixed and then reversing the roles. So, i n consideration of (3.3.6) and (3.3.7) for example, at each grid point we might first update v (3.3.6) while keeping v 2 fixed and then update v 2 x using using (3.3.7) keeping Vj fixed. A n alternative strategy is to update one component over the entire domain using (3.3.6) while keeping the other fixed and then reverse the roles using (3.3.7) i n place of (3.3.6). We should expect these methods to work well when the variables v 1 and v 2 are not strongly coupled. A t any grid point equations (3.3.6) and (3.3.7) may be written: Vi \ Cj (3.3.9) A \ v 2 \c 2 ) where the 2x2 m a t r i x A and the right hand side ( c , c ) 1 2 T are particular to the chosen discretization. If the m a t r i x A is strongly diagonally dominant then we may deduce that the variables v\ and v 2 are not strongly coupled and that sequential smoothing may be considered for the relaxation procedure. Collective smoothing is performed by solving the above 2x2 system directly at each grid point so that the two variables are updated simultaneously. Computationally this is more expensive than sequential smoothing but in cases where the variables are strongly coupled i n one or other of the equations, this method is a necessary alternative to sequential smoothing. In our case, the added expense i n using this method is almost negligible as we have only a 2x2 system and being a direct solver, it is more accurate than sequential smoothing which provides only an approximate solution. Chapter 3.4 3. The Multigrid 37 Method Multigrid for time-dependent problems Consider the time-dependent problem Cu(t) = f(t) on ft (3.4.1) w i t h the time-independent boundary condition: u = 0 on dn (3.4.2) A s s u m i n g a uniform discretization i n time w i t h time-step At, the discrete problem at the (n + l ) s t time step may be written: ^W + 1 < + = A n+1 on n (3.4.3) h where 1 = 0 on dtt (3.4.4) h A t each time step we may consider using L M G to solve the problem for u^. A t the first time step we take u° = 0 on the entire domain as an i n i t i a l guess to the solution. h C l e a r l y at each successive time step it would be pointless to use zero as an initial guess since u n can be considered to be a better approximation to u n + 1 than zero. If the time step is small enough, these approximations may be good enough to allow us to use only one L M G iteration per time step while still m a i n t a i n i n g reasonable accuracy. A n alternative to L M G is the F M G algorithm. Gendler [4] (following B r a n d t ) and H a c k b u s c h [5, p.272] suggest modified F M G F A S algorithms for dealing w i t h time dependent problems. Hackbusch [5, p.273] also presents criteria under which the proposed algorithm will converge. In our case, the F A S component is not required since we are dealing w i t h a linear problem, however, the procedure outlined below follows the same philosophy. Chapter 3. The Multigrid 38 Method A g a i n the question uppermost i n our minds is how information from the previous time step can best be used as an i n i t i a l guess for the current time step. Consider the following: 1. Inject the solution u n down to every level. 2. Solve exactly on the coarsest level to obtain u 3. Perform a high order interpolation of the n + 1 correction, u n + 1 — u, n to this solution and add it to the solution on the next finest level. 4. Perform one L M G cycle to obtain u . n+1 5. Repeat 3. and 4. until arriving at the finest level. Hence, at each time step, the problem may be viewed in terms of obtaining a correction at each level to the solution from the previous time step. T h i s prompts us to define the following algorithm: A l g o r i t h m 3.3 ( M F M G ) U ); Procedure M F M G ( / , 7 , 2 Begin for / :=p - 1 to 1 do un R E L A X ( £ 1 INJ« ); + 1 U 1 = /i); uf for / := 2 to p do begin u n+l n+l ;= u? + INT(u7-1 + 1 ] Chapter 3. The Multigrid L M G ( / , 7, u u 39 Method u ); {-> 2 end; End; where INJ represents the straight injection operator and I N T is a high order interpolation. T h i s algorithm involves more work than the L M G a l g o r i t h m but only on the coarser grids so the cost is not significantly higher. It is obviously hoped that the additional expense results i n a more accurate solution - this is yet to be established. In the M F M G algorithm we are m a k i n g the assumption that the 0(2h) tion error arising from our a p p r o x i m a t i o n to u initial guess to u step. n+1 2 discretiza- on the next to finest level is a better n + 1 than the 0(h) a p p r o x i m a t i o n to u obtained at the previous time 2 n It is not clear why this should be the case. In the limit as At —> 0 we can be certain that L M G should provide the more accurate solution since i n this case the O(Ai) error becomes insignificant i n comparison to the 0(2h) 2 error. A n estimate of the 0(2/i) discretization error would therefore be of use. T h e discretization errors, 8^2 and Sfj , at each grid point on the finest and next-to-finest levels respectively, are 1 given by % := - u 6}} ':= u% - ^ h (3.4.5) 2 i: ClJ ~ (2h) 2 U i j Together these equations lead to the following: ClJ (3.4.6) Chapter 3. The Multigrid 40 Method We w i l l use an r.m.s. representation of so that the discretization error on the next- to-finest level is approximated by: ^ = = ^ ( E ( < - < ) 2 ] (--) 3 4 9 Notice that such an a p p r o x i m a t i o n is easily obtained using the M F M G a l g o r i t h m where the solution (not just a correction to the solution) is computed at every level. In test cases where the exact solution is known, we can compare the exact r.m.s. change i n the solution from one time step to the next to the discretization error, 8 , on the 2h next-to-finest level. In cases where the discretization error is larger t h a n the change i n the solution we expect L M G to provide a better solution t h a n M F M G , i n Chapter 5 we w i l l show this to be the case., Chapter 4 The Elastic P r o b l e m Analytic Solution 4.1 First let us recall the governing equation w i t h its boundary conditions: D V w + jw = p (4.1.1) w ~ V"tt> = 0 on dfl (4.1.2) where In this section we shall seek an analytic solution to which we may later compare our numerical results. First we must choose a forcing function, p. Consider the case of periodic loading on a square d o m a i n where two adjacent sides of the domain lie along a set of Cartesian axes x and y. Assume we have a periodic topography with wavelength \ x i n the re- direction and X i n the y-direction, so that the topographic displacement, hr, at any y point (x,y) is given by: (4.1.3) where ho is the m a x i m u m topographic amplitude. Hence, if p represents the density c of the rock forming the surface topography we have: We now assume that the plate deflection responds in such a way that: (4.1.5) 41 Chapter 4. The Elastic where w 0 42 Problem is a constant to be determined. In other words, it is assumed that the plate deflection is also periodic w i t h a period equal to that of the loading function. Notice that i n order to satisfy the boundary conditions (4.1.2) we must insist that 2L/X X a n d 2L/X y be positive integers where L represents the length of the domain. Now, substituting (4.1.4) a n d (4.1.5) into (4.1.1) we obtain: -i-i 2 (4.1.6) which represents the amplitude of the deflection. W h i l e the application here is obviously to topographic loading such as mountain ranges a n d valleys, equation (4.1.4) may also be seen to represent a simple model for the formation of sedimentary basins. T h i s requires the assumption that A = X = 2L x y so that the forcing increases towards the centre of the plate. We therefore assume that a basin of shape hj is formed by some process such as thermal contraction so that w represents the subsequent, deflection of the plate due to sedimentation. However, we shall proceed w i t h the assumption that we are attempting to model deflection due to topographic loading. For most practical applications we will be considering the case of m o u n t a i n ranges which extend the length of the domain i n the y-direction a n d vary periodically i n the x-direction so that X = 2L a n d nX = 2L hence, X = nX where n is a positive y x y x integer. W e may therefore write (4.1.6) i n the form: w 0 ho 2TT\4 K.) D \ ( \ Pc9J V 1 y 1 n / 2 / 7 \p g (4.1.7) c where now the two terms i n the denominator are nondimensional. Turcotte & Schubert [10, p.123] point out that the quantity (D/pcg) / has dimensions of length a n d is 1 4 proportional to the natural wavelength for flexure of the lithosphere. Let us first consider the case where the wavelength of the forcing function is much greater t h a n Chapter 4. The Elastic 43 Problem the n a t u r a l wavelength so that: (1 + (4.1.8) then we have: Peg h•0 { (4.1.9) 7 where w represents the amplitude of deflection resulting from an isostatic adjustment l of the lithosphere to the loading. In other words, the lithosphere may be considered to have no rigidity and to be i n hydrostatic equilibrium. Alternatively, if the wavelength of the forcing is much less than the n a t u r a l wavelength, we have: so that: < and h 0 the deflection of the lithosphere is negligible. In other words, the rigidity of the plate totally supports the load. The above discussion is crucial to the understanding of future numerical results. In cases where the plate's rigidity is of importance, the biharmonic operator plays a heavy role in the governing equation and we may therefore consider the problem to be tougher to solve. A t the opposite end of the spectrum, i n cases where the plate's rigidity is not significant, we have the equivalent of a small parameter m u l t i p l y i n g the highest derivative in the equation. T h i s yields a singularly perturbed problem which, although a very simple one, could again give rise to difficulties i n obtaining a numerical solution. A s we proceed, we shall therefore bear in m i n d that some difficulty may be encountered in these extreme cases. Turcotte &i Schubert [10, p.123] examine the degree to which a load is compensated by the plate's flexural rigidity. In a similar vein, Cha.pter 4. The Elastic 44 Problem o o o o 0.0 0.5 NONDiMENSIONflL WRVELENGTH 2.5 3.0 Figure 4.5: Degree of support as a function of nondimensional topographic wavelength we will compute the degree to which the plate's rigidity supports the load as this will provide us w i t h some measure of the importance of the biharmonic operator i n the solution of the equation. Let R = (u>, — w )/w. so that for iy — 0, R. ~ 1 and we have 0 i 0 the case where the plate's rigidity totally supports the load, the effects of buoyancy are low and the biharmonic operator plays an important role in the solution of the governing equation. However, for w 0 — to,-, R — 0 and we have the case where the plate's rigidity is ineffective against the loading, the effects of buoyancy are significant and the biharmonic operator does not play an important role i n the solution of the governing equation. Substituting for WQ and w, we have: (4.1.11) Figure 4.5 represents the degree of support, R, as a function of nondimensional Chapter 4. The Elastic 45 Problem L(m) 1 x 10 I 2.5 x 10 5.0 x 10 1.0 x I O 5 6 = = D (Nm) 1 x 10 23 1 x 10 92.7 11.3 0.8 0.05 5 L D 21 25 99.9 98.8 83.2 44.3 4.7 L e n g t h of Plate Flexural Rigidity Table I: Degree of support (Elastic M o d e l ) topographic wavelength. Hence, for a 50% support of the loading by the plate's rigidity we have: 1 + = 2*1-1 rt j (4.1.121 \ 7 2 We may also express R as a function of the plate length, L, as follows -l 1 R = r VD/ Vl (4.1.13) n' Table I gives appropriate values of R for various values of D and L where a value of n — 2 has been assumed. T h e range of values of D covers most of those found i n the literature for geophysical applications. T h e range of values of L covers everything from sedimentary basins to continental plates. T h e degree of support depends directly on n , hence, for a given length, L, and rigidity, D, larger values of n correspond to greater support. Following Turcotte & Schubert [10, p. 123], we take p , the density of m the mantle, to be 3300 kg m ~ m~ . 3 3 and p , the density of the surface rock to be 2800 kg c Hence, 7 has the value 4900 k g m ~ . 2 2 Chapter 4. The Elastic 4.2 Numerical 46 Problem Solution Recall that we propose to solve the problem numerically by reducing (4.1.1) to two coupled equations so that at the grid point D^Vij+'ywij V Wij 2 we have: = (4.2.1) Pij = Vij (4.2.2) w i t h boundary conditions to = w = 0 o n 80, (4.2.3) We discretize the L a p l a c i a n operator using the standard five-point formula: V^u' t ; = —(Wij-x + w lJ+1 + w_ z ltj + w l+li3 - 4w ) li0 (4.2.4) where h represents the grid-spacing which we have assumed to be uniform over the entire domain. W e use full weighting for restriction a n d bilinear interpolation for prolongation. R e l a x a t i o n is carried out using Gauss-Seidel w i t h red-black ordering, see Hackbusch [5, p.51]. We assume that the coarsest grid contains only one interior grid point so that one sweep of Gauss-Seidel solves the discretized problem exactly on the coarsest level. Finally, we must decide on sequential or collective smoothing as outlined in section 3.3. Recall that our choice depends on the degree to which the variables are coupled i n each of the equations. H a v i n g chosen an appropriate discretization of the L a p l a c i a n , we may now write our equations i n the form (4.2.4) where the 2x2 m a t r i x , A, becomes: ' 4/h 2 - i j U \ (4.2.5) 1 4/A / 2 Notice that the primes have been brought back into use t o stress that numerically we deal w i t h nondimensional quantities, see section 2.5. Sequential smoothing should be Chapter 4. The Elastic L(m) 2.5 x 10 1.9 ijU 47 Problem L 7'/D' 5.0 x 10 31.0 5 = = 5 1.0 x I O 4.9 x 10 6 2 5.0 x I O 3.1 x I O 6 6 L e n g t h of P l a t e Nondimensionalized j/D. Table II: j' / D' as a function of plate length effective when A is diagonally dominant, this requires that: h D' 2 anc J ? > 1 T h e second inequality is clearly satisfied for all sensible choices of h. In the worst case (i.e. on the coarsest level) we have 4 / / i = 16. For sequential smoothing to be effective 2 we therefore require j'/D' < 16. Table II indicates values of "y'/D' for various choices of L where D is taken to be 1 x 1 0 25 (i.e. the best choice for which we might expect the inequality to hold). It can be seen that sequential smoothing should be considered for only the smaller domains since the variables become strongly coupled for L > 500 k m . Notice that smaller values of D would result i n an even stronger coupling for a given L. It is therefore clear that sequential smoothing should be avoided and hence we choose collective Gauss-Seidel for the relaxation procedure throughout. A s was pointed out earlier, this does not lead to any significant additional expense since we are dealing only w i t h a 2x2 system. There is also the added benefit of obtaining an exact solution of the 2x2 system at each grid point. Chapter 4.3 4. The Elastic 48 Problem N u m e r i c a l Results We now examine the performance of the suggested m u l t i g r i d procedure for various lengths of domain and choices of flexural rigidity. We first consider periodic loading for which an analytic solution is available making it possible to determine the accuracy of the numerical results. We shall then proceed to examine Gaussian loading for which an analytic result is not available. 4.3.1 Periodic L o a d i n g In examining the performance of a numerical procedure there are two things to consider. F i r s t , the rate of convergence of the numerical solution to the exact solution of the discrete problem and second, the accuracy of the numerical solution which, given convergence of the solution, depends on the discretization error or in other words the difference between the exact solutions to the discrete and continuous problems. Ide- ally, we would like to have a method which provides us w i t h a good combination of b o t h - clearly one is not much use without the other. A s we have an analytic solution available, we can easily find a good measure of the accuracy of the numerical solution. A root mean square error, E , expressed as a percentage of the m a x i m u m displacement w i l l be used in this capacity so that: (4.3.6) where w l] is the numerical solution, u>, is the exact solution and w 7 0 is the amplitude of the deflection derived from the exact solution. W h e n considering periodic loading where the solution may become zero at interior points of the domain we must be careful how we measure convergence. T h e easiest Chapter 4. The Elastic 49 Problem method is to examine convergence at one particular point and take this to be representative of the overall convergence. Hence, it is imperative that we choose a point where m a x i m u m deflection, rather t h a n zero deflection, is expected - knowing the exact solution, this w i l l be possible. We w i l l terminate the iterative procedure at the k th iterate when (k) (k- < (4.3.7) TOL (k) where the grid point all cases TOL is a location where m a x i m u m deflection is expected. is chosen to be 1 x I O - 4 In so that the iterative procedure is terminated when the solution at a point of m a x i m u m deflection has converged to w i t h i n one ten thousandth of its computed value. In Tables III, I V and V , we present results for various sizes of plates, rigidities and topographic wavelengths. flexural T h e number of iterations, I T N , required for convergence, the r.m.s. error, E , and the predicted, R , and exact, R , degrees of p e support, are given in each case. T h e m u l t i g r i d hierarchy of grids consists of 4 levels, the finest, having 15 interior grid points in each direction and the coarsest having only 1 interior grid point. T h e L M G algorithm w i t h V ( 2 , 2 ) cycles was employed throughout. Our results clearly indicate the relevance of the degree of support to the accuracy and convergence rate of the numerical method. We see that for small degree of support, the m u l t i g r i d iteration generally converges quickly to a solution which is a very good approximation to the exact solution. T h i s is worthy of note since for a small degree of support, the nondimensionalized equations have a small parameter m u l t i p l y i n g the highest derivative i.e. the biharmonic operator. We therefore have a singularly per- turbed problem for which our m u l t i g r i d method provides a good solution. However, in consideration of the discretized equations we see that i n the limit as the degree of Chapter 4. The Elastic D n 1 1 1 1 1 1 2 4 50 Problem x x x x x x 10 10 10 10 10 10 ITN E(%) Rp(%) 4 5 5 4 4 4 0.13 1.09 1.18 3.12 5.42 5.46 11.1 92.6 99.9 57.2 99.2 99.9 21 23 25 21 23 25 R (%) e 11.3 92.7 99.9 59.5 99.3 99.9 Table III: L = 2.5 x 1 0 m , Periodic L o a d i n g 5 n D 1 1 1 1 1 1 2 4 x x x x x x 10 10 10 10 10 10 ITN E{%) 5 3 5 4 4 4 0.01 0.52 21 23 25 21 23 25 1.16 0.42 4.87 5.46 R (%) P 0.8 43.7 98.7 7.7 89.3 99.9 R (%) e 0.8 44.3 98.8 8.4 90.2 99.9 Table I V : L = 5.0 x 1 0 m , Periodic L o a d i n g 5 n ITN D 1 1 1 1 1 1 2 4 x x x x x x 10 10 10 10 10 10 21 23 25 21 23 25 3 3 ' 5 3 4 4 . E(%) 0.001 0.04 0.98 0.03 1.87 5.36 Rp(%) 0.05 4.7 82.9 0.52 34.3 98.1 R (%) e 0.05 4.7 83.2 0.57 36.5 98.3 Table V : L = 1.0 x 1 0 m , Periodic L o a d i n g 6 n D ITN E R R p e = — = = = = 2L/X Flexural Rigidity N u m b e r of iterations for convergence Relative error Predicted degree of support E x a c t degree of support X Chapter 4. The Elastic 51 Problem support tends to zero, equation (4.2.2) plays no role in the c o m p u t a t i o n and equation (4.2.1) reduces to: w tj = ^ (4.3.8) 7 which is the solution of the reduced continuous equation representing an isostatic adjustment of the lithosphere to the applied load. T h e iterative Gauss-Seidel procedure therefore becomes an exact solver for the problem which is further simplified by the fact that p = 0 on the boundaries so that there are no boundary layers. It is therefore no surprise that the numerical scheme deals adequately w i t h the singularly perturbed problem. It should be noted that if we had not had p = 0 on the boundaries, the procedure would still have been good owing to the local nature of the Gauss-Seidel iteration. As the degree of support increases for a given L the procedure requires a larger number of iterations to converge and the accuracy of the solution is seen to decline. It is also observed that for a given L and D , increasing n also increases the error i n the solution. T h i s is likely caused by a combination of two factors. F i r s t , increasing n increases the degree of support and therefore the importance of the biharmonic operator, and second, it increases the frequency of the forcing function so that it is not as well represented on the finest grid. O v e r a l l , we can be pleased w i t h the accuracy and convergence rate. Note that 5 V ( 2 , 2 ) cycles require only 20 Gauss-Seidel iterations on the finest grid so we have a very efficient solver in hand. A l s o , the r.m.s. error i n the solution is at worst 5% of the m a x i m u m deflection which again can be considered as reasonable. A reduction in this error may be achieved by decreasing the stepsize on the finest, grid. Chapter 4. The Elastic L 1.0 5.0 2.5 X X X 52 Problem D 10 10 10 L D GS-itn MG-itn R e 1 x 10 1 x 10 1 x 10 6 5 5 — = = = — GS-itn MG-itn 5 58 94 3 3 5 21 23 25 R {%) e 0.05 44.3 99.9 Length of P l a t e Flexural Rigidity N u m b e r of Gauss-Seidel iterations N u m b e r of M u l t i g r i d iterations E x a c t degree of support Table V I : M u l t i g r i d vs. Gauss-Seidel 4.3.2 M u l t i g r i d vs. Gauss-Seidel In order to appreciate the efficiency of the multigrid solver, we compare it to straight Gauss-Seidel where we repeatedly call the R E L A X routine on the finest level and do not perform coarse grid corrections. Three cases have been tested covering an extensive range of values of R. In all cases we take n = 2 and m a i n t a i n TOL at 1 x 10~ . T h e 4 results are presented i n Table V I where G S - i t n is the number of straight Gauss-Seidel iterations and M G - i t n is the number of m u l t i g r i d iterations required for convergence. The work involved i n relaxing the 4 levels of one m u l t i g r i d V ( 2 , 2) cycle is equivalent to approximately 5 Gauss-Seidel iterations on the finest, level - this gives us some means of comparing the two methods. We see that for a small degree of support where the biharmonic operator is not significant, straight Gauss-Seidel proves to be more efficient than multigrid. However, we suspect that the degree of support does not have to become too large before multigrid performs"relatively more efficiently. Multigrid certainly proves to be more efficient when the degree of support is at 44% or higher. We can therefore claim w i t h some confidence that i n cases where the flexural rigidity is significant, the m u l t i g r i d method described proves to be a very efficient solver for the problem. Chapter 4. 4.3.3 The Elastic 53 Problem Gaussian Loading We now t u r n our attention to a loading function for which an analytic solution may not be obtained. Consider a Gaussian distribution for the forcing so that (4.3.9) where (x - 0.5) + (y- 0.5) 2 2 T h e application here is to the formation of sedimentary basins where a basin of depth hr is formed by thermal contraction. W e therefore attempt to predict the subsequent deflection, w, due to the sedimentary infill of this basin. T h e parameter o is chosen so that the loading function at the edge of the domain is close to zero and the function hr{r) is close to continuous at r — 0.5. T h e parameter a in some sense represents the wavelength of the loading. Notice that m a x i m u m deflection w i l l occur at the midpoint and so again we may compute the degree of support although we no longer have an exact value to which we may compare it. W e again expect that smaller wavelengths i n the loading result in greater support from the plate's rigidity. W e therefore consider two choices for a, namely 0.04 and 0.08. B o t h choices are reasonable from a numerical standpoint in that the function / i j ( r ) is not too eratic when seen at the edges of the finest grid. T h e results are presented in Tables V I I , V I I I and I X . We again observe that in some cases faster convergence is achieved when the degree of support is lower. It can be seen that the larger value of a yields lower degrees of support. T h i s is as expected since a larger value of a results i n a larger wavelength in the forcing which in t u r n results in a lower degree of support as was seen for periodic loading. A g a i n we can be satisfied w i t h the convergence rate but of course we no longer have any means of assessing the accuracy of the solution, however, the fact Chapter 4. The Elastic 54 Problem a 0.04 0.08 D ITN 4 5 5 4 5 5 1 1 1 1 1 1 x x x x x x Rp{%) 10 10 10 10 IO 10 21 23 25 21 25 25 28.3 85.5 99.8 13.3 76.9 99.7 Table V I I : L = 2.5 x 10 m , Gaussian L o a d i n g 5 a 0.04 0.08 ITN D 5 5 5 5 5 5 1 1 1 1 1 1 x x x x x x R (%) v 10 10 10 10 10 10 21 23 25 21 25 25 4.6 52.7 96.8 1.4 32.4 94.9 T I L L = 5.0 x 10 m , G aussian . s a 0.04 0.08 5 ITN 3 4 4 3 5 4 D 1 1 1 1 1 1 x x x x x x 10 10 10 10 10 10 21 23 25 21 25 25 Rp(%) 0.3 17.2 74.8 0.1 6.9 60.6 Table I X : L = 1.0 x IO m , Gaussian L o a d i n g 6 a ITN D R p = — = = Forcing function parameter, see (4.3.9) N u m b e r of iterations for convergence Flexural Rigidity Predicted degree of support Chapter 4. The Elastic Problem 55 that it behaves qualitatively as expected seems to indicate that there are no serious difficulties. Chapter 5 The Viscoelastic Problem 5.1 Analytic Solution W e begin by restating the governing equation w i t h its initial and boundary conditions: , _-.\ ( w 74 • , DV*w + 7 ^ - +wj P (5.1.1) = ^ + p where, assuming the exact solution is known, w = (t = 0) We o n (1 (5.1.2) and w = V u ; = 0 on dn (5.1.3) 2 A s for the elastic problem, we seek an analytic solution to which we may compare future numerical results. W e again consider a forcing function which varies periodically over the domain but we must now include some time dependence. In geophysical applications such as the formation of sedimentary basins, the load increases from zero w i t h time. A s an approximation to this situation, bearing in m i n d that we need an analytic solution, we assume a loading function of the form: p = gh e^-^ Pc 0 s i n s m ^ ( 5 . L 4 ) for which we have a corresponding deflection: w = sin ( ~ ) 56 sin (5.1.5) Chapter 5. The Viscoelastic 57 Problem where w = (m + 0 and t e l)p gh c Q 16ZW 4 . o 1 + — 1 ) +(m + l ) (— -I -1 (5.1.6) 7 is the time at which the integration is terminated. Notice that at t = t , the e forcing reaches its m a x i m u m value of . Pmai = Pcgho sin (2irx \ . J s m \ (2-Ky\ (5.1.7) j which is the same as that considered for the elastic problem. As for the elastic problem, if we assume X = A^ = 2L we may regard (5.1.4) as x representing the load due to the infilling of a basin of depth V m ( '" l e ) / T - T h e basin, assumed to be formed by some process such as thermal contraction, therfore increases in depth exponentially over a period of t m i l l i o n years. e W i t h an appropriate choice of X and X the boundary conditions (5.1.3) may be x y satisfied. Notice also, that the deflection at t — 0 can be made arbitrarily small by an appropriate choice of m although large values of m result i n larger rates of increase i n the forcing which could have serious repercussions when it comes to finding a numerical solution. We shall return to this problem later but in the meantime we shall, on the above grounds, consider (5.1.4) to be an acceptable analytic solution. T a k i n g X = riX = 2L as i n the elastic case we may again derive an expression y x for the degree, R to which the load is supported by the plate's rigidity. For the case } outlined above we have: R 1+ - .TTJ ^ 1+ \DJ V mJ V l + n (5.1.8) 2 Table I X gives appropriate values of R for various choices of L and D where values of n = 2 and m = 0.01 are assumed. If either m or n is increased, then R will increase for a given L and D. A s can be seen, the case L = 1.0 x 1 0 m is not particularly 6 interesting as the lithosphere is close to being i n isostatic adjustment for all choices of Chapter 5. The Viscoelastic 58 Problem L(m) 1 1 x 10 2.5 x 10 5.0 x 10 1.0 x 10 5 5 6 21 D (Nm) 1 x 10 1 x 10 11.2 0.78 0.05 92.6 44.1 4.7 23 0.12 0.0078 0.0005 L — L e n g t h of P l a t e D = F l e x u r a l Rigidity- 25 Table X : Degree of support (Viscoelastic M o d e l ) D. We therefore restrict ourselves to the two smaller domains for the remainder of this discussion since together they cover an extensive range of values of R. 5.2 Numerical Solution Numerically the problem is solved by first rewriting (5.1.1) i n the form (2.7.1) as follows: w (5.2.1) =u (5.2.2) Vu 2 (5.2.3) =v These equations may now be viewed as a differential algebraic equation ( D . A . E . ) i n time. T h e time component can be regarded as differential in w and algebraic i n u and v. For a more in-depth discussion on the numerical solution of D . A . E . ' s see Ascher [l]. We now choose a discretization for the time derivative. A s a first choice consider B a c k w a r d E u l e r ( B . E . ) which is a first order method. K n o w i n g quantities at time t , n the p r o b l e m at time t n+1 may be written: ,n+l _ (5.2.4) At w + U n+l T (5.2.5) Chapter 5. The Viscoelastic 59 Problem X7 u 2 Eliminating w n+1 = n+1 (5.2.6) n + 1 v from (5.2.5) using (5.2.4) we have: + 1) u + 7 ( ^ = v Vu 2 n+1 n+l = g (5.2.7) n+1 (5.2.8) n+l where g = - n+1 +p- w 1 r r is k n o w n at each time step. (5.2.9) n We now have two coupled elliptic equations for u and v. R e c a l l that i n section 2.6 we showed that our boundary conditions may be w r i t t e n equivalently as: v = 0 on dn u = 0, (5.2.10) T h u s , at each time step, our problem may be w r i t t e n in the form (2.7.6) w i t h boundary conditions (2.7.7) where ui = u a n d u = v . Hence, we may use a m u l t i g r i d n+1 algorithm to solve for u n + 1 n+1 2 at each time step and then use (5.2.4) to solve for w . n+1 A s an alternative, following Ascher :[1, p.12], we may consider computing w at time steps t while u and v w i l l be computed at n t n + 2. W i t h this formulation we may use a centred difference formula ( C . D . F . ) for the time discretization hence giving us a second order method i n time which requires no more work or storage than B a c k w a r d Euler. So we have w -w n+1 u ? n+ = n At ^-1 +"+§ u VV 2=i; § + n+ v = 1 +p i n+1 . (5.2.11) ' (5.2.12) (5.2.13) Chapter Now, 5. The Viscoelastic 60 Problem writing w = = n+ n+i = at each time step i n + 2 , -(w +w ) 2 l -(Mu * n (5.2.14) +2w ) n+ n we use m u l t i g r i d to solve the system: + + i) 7 V u 2 n + 2=u = n + (5.2.15) § (5.2.16) where (?"' = ^ - ^ (5.2.17) - Iiy" + T T is known at each time step. H a v i n g solved for we use (5.2.11) to solve for u i n+ w . n+1 Second order accuracy may also be obtained ( w i t h a little extra work) using the two-step backward difference formula ( B . D . F . ) : u 1 / 3 ^ , - — -w A t V2 n + 1 n+1 _ „ 1 - 2w + - u ; " ) 2 J 7 1 n 1 K (5.2.18) .' T h i s method also has greater storage requirements than B a c k w a r d Euler or the C . D . F . approach a n d is slightly more awkward to implement if a change i n step size, A t , is contemplated. U s i n g (5.2.18) we have the following system of equations: DV v 2 + n+1 + l) w 7 Vw = v 2 n+1 + pn + l _ _ n n+1 = .-g n+1 (5.2.19) (5.2.20) n+1 where now n+l g „~n+l n+l = t T A w 3 i + _ n-l w 3 (5.2.2T A g a i n , this system may be written in the form (2.7.6) and our m u l t i g r i d algorithm used once more. Chapter 5. The Viscoelastic 61 Problem We use the same 5-point operator to discretize the L a p l a c i a n and again use full weighting for the restriction and interpolation operators. Collective Gauss-Seidel is called u p o n once more for the relaxation procedure. 5.3 Numerical Results We begin by considering a case where an analytic solution is known thereby giving us a means of assessing the numerical results more objectively. Later we shall consider cases which are more meaningful geophysically but for which there is no analytic solution readily available. 5.3.1 Periodic Loading In the results which follow we assume a forcing function of the form (5.1.4). W i t h so many parameters (D, L , m, n , r ) contained i n the present formulation of the problem we must make some restrictions on their variability i n order to-make meaningful numerical •tests. We will allow n to take on only the value 2 since this results i n a forcing function which is both interesting and well represented on the fine grid. In other words, we avoid the uninteresting case of n = 1, which is similar i n form to Gaussian loading, as well as larger values of n which may result i n some error being introduced to the solution because of poor representation of the forcing function on the finest grid. Notice that i n restricting ourselves to the case of n = 2 we are in some sense losing geophysical applicability since the application to sedimentary basins requires the assumption that n = 1. T h e case n = 2 is applicable to topographic loading although it is hardly reasonable to assume that the topography grows exponentially w i t h time when in fact surface topography generally decreases i n size after formation due to erosion. A more meaningful forcing function, from the geophysical standpoint, w i l l be considered Chapter 5. The Viscoelastic 62 Problem in section 5.3.3. We must also provide appropriate choices for m . Recall that m not only controls the magnitude of the forcing at t = 0 but also the rate of increase i n the forcing. Increasing m yields a smaller forcing at t = 0 w h i c h is desirable, at least from a geophysical standpoint, but this also increases the rate of growth of the forcing w h i c h provides us w i t h a tougher problem to solve numerically. a happy balance. Values of m Clearly, we wish to find — 0.01 and 0.03 result in an i n i t i a l forcing having approximately 37% and 5% of its final value, respectively. T h e second of these is more reasonable geophysically but we expect it to be more difficult to solve numerically. Following N u n n & Sleep [7, p.607], we take the relaxation parameter, r , to have the value 1 x 10 years. 6 N u n n & Sleep [7, p.588] also indicate that any physical theory should allow continual subsidence for a period of at least lOOMa i n the modelling of sedimentary basins. Integration w i l l therefore be carried out from 0 to lOOMa, however, all equations have been nondimensionalized w i t h respect to a time scale T = lOOMa so that numerically, the equations are integrated from 0 to 1. We begin w i t h a comparison of the L M G and M F M G algorithms for time dependent problems as outlined in section 3.4. Tables X I , X I I and X I I I present the r.m.s. error, Ei so and Ed , calculated as i n (4.3.6), for the solution w, and its derivative u, for er various values of D. For this comparison we allow m to have only the value 0.01. A n r.m.s. estimate of the discretization error, S , 2h calculated as i n (3.4.9), and an r.m.s. representation of the exact difference in the derivative u, at time t = lOOMa are also presented. T h e L M G and M F M G algorithms were tested using B a c k w a r d Euler, the C . D . F . m e t h o d and the two-step B . D . F . method all w i t h a time step of l O M a . One of the first things that can be noticed from the results is that a better estimate of the derivative does not. necessarily lead to a better solution. A t first this may appear somewhat strange but we must bear in m i n d that we are dealing w i t h b o t h a time and Chapter 5. The Viscoelastic E i(%) MFMG LMG so D 1 x 10 1 x 10 1 x 10 2 1 2 3 2 5 63 Problem 0.029 0.470 2.839 0.081 0.259 1.420 E (%) LMG 2.506 2.323 0.088 ph MFMG (xlO- ) A (xlO- ) 2.553 2.127 1.669 0.0079 0.1794 0.3444 0.3406 0.3029 0.0251 der 4 2 2 Table X I : B a c k w a r d Euler, L M G vs. M F M G E (%) 1 x 10 1 x 10 1 x 10 MFMG LMG D 21 2 3 2 5 0.066 0.123 1.119 0.050 0.090 0.144 ph Ed sol LMG 0.063 0.071 0.180 MFMG (xlO" ) A* (xlO" ) 0.037 0.146 1.808 0.0120 0.2404 0.0330 0.3406 0.3029 0.0251 2 2 Table X I I : C . D . F method, L M G vs. M F M G E (%) 1 x 10 1 x 10 1 x 10 MFMG LMG D 21 2 3 25 0.034 0.413 0.248 ph Ed sol 0.002 0.197 1.153 LMG 0.139 1.624 0.176 MFMG (xlO- ) A* (xlO" ) 0.163 0.035 1.808 0.0113 0.2268 0.0347 0.3406 0.3029 0.0251 2 Table X I I I : Two-step B . D . F . , L M G vs. M F M G D = E l = Eder = ph = so = Flexural Rigidity Relative E r r o r i n Solution Relative E r r o r i n Derivative Discretization E r r o r E x a c t change i n derivative 2 Chapter 5. The Viscoelastic 64 Problem space discretization and hence some cancellation of the two discretization errors could well lead to some some spurious results. A very consistent pattern may be observed for results obtained for D = 1 x 1 0 N m . 2 5 In these cases, L M G consistently provides a better approximation to b o t h the solution a n d the derivative t h a n M F M G . It can also be noticed that i n these cases the discretization error on the next-to-finest level is significantly larger than the change i n the derivative and so we expect L M G to perform better t h a n M F M G as discussed i n section 3.4. For other values of D, things are not quite so clear cut. For D = 1 x 1 0 N m 2 1 where the discretization error is consistently smaller than the change i n the derivative, M F M G is expected to at least provide a better estimate of the derivative if not for the solution, but this occurs i n only one out of the three cases. O v e r a l l , L M G provided a better a p p r o x i m a t i o n t h a n M F M G to the solution i n 6 out of the 9 cases and a better a p p r o x i m a t i o n to the derivative also i n 6 out of 9 cases. We w i l l therefore conclude that L M G being the cheaper method should be employed for the remainder of our discussion on the viscoelastic problem. In the following discussion we examine the effect of altering m and attempt to determine how the degree of support, and therefore the importance of the biharmonic operator, affects the accuracy of the solution. We compare results using the B . E . , C . D . F . a n d two-step B . D . F . methods w i t h one V ( l , l ) cycle per time step except the first where five V ( l , 1) cycles were used to ensure the solution got off to a good start. The results are presented in Tables X I V , X V and X V I . We can see that as the degree of support increases, the accuracy of the solution declines consistently. Also, for a given D, the accuracy decreases when m is increased. T h e reason for this is two-fold. F i r s t , increasing m increases the degree of support and therefore the significance of the biharmonic operator and second, it increases the rate of growth of the forcing function Chapter 5. The Viscoelastic D m 1 1 1 1 1 1 0.01 0.03 x x x x x x 65 Problem E i(%) so 10 10 10 10 10 10 2 1 2 3 2 5 21 2 3 2 5 Eder{%) R {%) v 0.081 0.259 1.420 2.506 2.323 0.088 6.741 6.092 1.615 0.600 1.350 6.284 Table-'-XIV B a c k w a r d E u l e r , m D 0.01 1 x 10 1 x 10 1 x 10 2 1 1 x 10 1 x 10 1 x 10 2 1 0.03 E {%) soi 2 3 2 5 2 3 2 5 0.050 0.090 0.144 0.458 1.193 1.490 E „{%) d Re(%) 0.01 10.8 92.4 0.13 11.2 92.6 0.37 27.0 97.3 -0.68 25.2 97.1 A t = lOMa Rp(%) 0.063 0.071 0.180 0.25 11.3 92.6 0.239 0.997 1.423 1.36 28.6 97.4 R (%) e 0.13 11.2 92.6 0.37 27.0 97.3 Table X V C . D . F . m e t h o d , A t = 1 0 M a m D 1 x 10 1 x 10 1 x 10 1 x 10 •1 x 1 0 1 x 10 0.01 0.0? E (%) Ed {%) 0.034 0.023 0.248 0.139 0.183 0.176 1.047 1.684 1.447 sol 21 2 3 2 5 2 1 2 3 2 5 0.264 0.413 0.122 eT Table X V I : Two-step B . D . F . , m D E i Eder R R so p e = — = = = = &,(%) Re(%) 0.10 11.2 92.6 0.06 27.6 97.3 0.13 11.2 92.6 0.37 27.0 97.3 A t = lOMa Forcing function parameter, see (5.1.4) Flexural Rigidity Relative E r r o r i n Solution Relative E r r o r i n Derivative Predicted degree of support E x a c t degree of support Chapter 5. The Viscoelastic Problem 66 so that the solution from the previous time step is not as good an approximation to the current iterate. W e expect the latter reason is the more significant since if we do not start w i t h a good approximation to the solution at the current time step, we cannot expect to achieve reasonable accuracy from only one V ( l , l ) cycle per time step. We should expect this problem to be remedied by either selecting a smaller time step or choosing a higher order discretization for the time derivative. In comparing C . D . F . and B . D . F . to B . E . we see that the higher order discretization does much to reduce the error i n b o t h the solution and its derivative. Later we shall examine the effect, of reducing the time step. C . D . F . consistently provides better results t h a n B . E . for the same amount of work m a k i n g C . D . F . a very attractive method. C . D . F . also provides less spurious results in that as D increases, the errors i n b o t h the solution and the derivative are seen to increase as expected. T h i s is not observed for either B . E . or B . D . F . . A comparison of C . D . F . and B . D . F . does not lead to any conclusive results. T h e y should b o t h provide approximations which are second order i n time for the deflection w, however, C . D . F . provides only a first order approximation to the derivative u, whereas B . D . F . again supplies a second order approximation. There is a simple fix for this potential deficiency in C . D . F . , see Ascher [1, p.13], however we have not used it because the theory.is not borne out i n the results. It is in fact C . D . F . which i n general provides better estimates for the derivative but these estimates have not led to better approximations to the solution. B . D . F . generally provides the better approximation to the solution. However, in consideration of the additional work and storage involved i n using the two-step method, we may conclude that of the two second order methods, C . D . F . is preferable here. Before discarding B . D . F . completely, though, we w i l l examine the effect of reducing the time step. Tables X V I I , X V I I I and X I X give results for all three methods. We restrict ourselves to the case m = 0.03 which has proved to result i n less Chapter 5. The Viscoelastic 67 Problem m D 0.03 1 x 10 1 x 10 1 x 10 Esol(%) E {%) der 2 1 2 3 2 5 0.198 0.899 3.563 3.629 2.973 0.353 R (%) Re{%) 0.09 25.8 97.2 0.37 27.0 97.3 P Table X V I I : B a c k w a r d Euler, A t = 5 M a m D 0.03 1 x 10 1 x 10 1 x 10 E i(%) so 2 1 2 3 2 5 0.120 0.298 0.278 E „{%) d 0.074 0.249 2.257 R (%) Re(%) 0.63 27.4 97.4 0.37 27.0 97.3 P Table X V I I I : C . D . F metb.od,..At = 5 M a m D 0.03 1 x 10 1 x 10 1 x 10 E (%) sol 2 1 2 3 2 5 0.067 0.100 0.126 Eder{%) £P(%) 0.308 0.453 0.263 0.33 27.1 97.4 R {%) e 0.37 27.0 97.3 Table X I X : Two-step B . D . F . , A t =-5 M a m D E i Eder R R so p e = = — — = = Forcing function parameter, see (5.1.4) F l e x u r a l Rigidity Relative E r r o r i n Solution Relative E r r o r i n Derivative P r e d i c t e d degree of support E x a c t degree of support accurate solutions. A s expected, reducing the time-step has led to more accurate results for a l l three methods. B a c k w a r d Euler, being a.'hrst order method, provides an O ( A t ) approximation (in time) to b o t h the solution a n d its derivative. So, if there were no error arising from the spatial discretization and no significant roundoff error, reducing the time step by a factor of 2 should reduce the errors i n the solution and its derivative by 50%. T h e errors we have presented, however, account for errors arising from b o t h the time and spatial discretizations and so we should not necessarily expect a 50% improvement i n 68 Chapter 5. The Viscoelastic Problem our results. We expect that reducing the time step by a factor of 2 for the C . D . F . method should reduce the error i n the solution by a factor of 4 whereas the error i n the derivative should be reduced by only 50% since C . D . F . provides a second order a p p r o x i m a t i o n to the solution but only a first order a p p r o x i m a t i o n to its derivative. B . D . F . , on the other hand, provides a second order a p p r o x i m a t i o n to b o t h the solution and its derivative and so we should expect the errors i n b o t h to be reduced by 75%. Once again, however, the above arguments hold only for the error arising from the time discretization. We should not necessarily expect to observe the above improvements i n accuracy since we also have an error arising from the spatial discretization. B a c k w a r d E u l e r i n fact provides close to a 50% improvement i n the approximations to b o t h the solution and its derivative. T h e C . D . F . method i n some sense behaves better t h a n expected in that a reduction of close to 75% is observed not only for the error i n the solution, but also for the error in the derivative when only a 50% improvement i n the latter is expected. B . D . F . provides a near 75% improvement for b o t h the solution and its derivative, as expected, except for the case D = 1 x 1 0 N m where the error 2 5 in the solution actually gets worse! A n explanation for this strange observation is extremely hard to find. O n comparing the two second order methods, we see that B . D . F . once more provides better estimates of the solution whereas C . D . F . provides better estimates of the derivative. In view of the added difficulty and additional expense i n using B . D . F . we shall from hereon restrict ourselves to the C . D . F . method. W e shall also discard Backward E u l e r since we clearly gain nothing in using this m e t h o d from the point of view of either accuracy or efficiency. Chapter 5.3.2 5. The Viscoelastic 69 Problem Integration with steady forcing We now consider the effect of allowing the forcing term to reach a steady state and observe the resulting deflection. Geophysically, this is of interest since, i n consideration of the evolution of sedimentary basins, we may assume that after a given time, the sedimentary loading reaches a steady state. Hence, for t > t we look for a solution of e the equation: DV u) + 7 (5.3.1) + w^j = P 4 where now, p is independent of time. T h e solution to this equation is: (5.3.2) w=7 in other words, an isostatic adjustment of the lithosphere. Therefore, if we continue to integrate after the forcing reaches a steady state, the solution should converge to the isostatic solution. O w i n g to the expense of integrating for long periods of time, we have considered only two cases. We allow the forcing to have the form (5.1.4) as we integrate from 0 to t which we e take to be lOOMa. For times greater t h a n t , the forcing function has the form: e P = Pcgh 0 . /2"7rx\ sin [~^~J . s / 27ry \ J m (5.3.3) so that it maintains the value it h a d at lOOMa. W e consider two cases having different values of L, D a n d m. T h e first takes L = 2.5 x 1 0 m , D = 1 x 1 0 N m and m = 0.03 5 2 3 while the second has L = 5.0 x 1 0 m , D = 1 x 1 0 N m and m = 0.01. We take 5 2 5 At = l O M a for b o t h a n d use the C . D . F . method. In Table X X we present computed values for the m a x i m u m deflection, w , at intervals of lOOMa for case (i). T h e max m a x i m u m deflection expected when the lithosphere is i n isostatic adjustment is 16800m. It is very clear that the numerical solution behaves as predicted and converges to the Chapter 5. The Viscoelastic 70 Problem t(Ma) 100 200 300 400 500 11995.21 16797.98 16798.74 16799.10 16799.76 Table X X : Steady state integration : case 1 t(Ma) 100 . 200 300 400 500 600 700 800 900 1000 9408.81 14790.80 16253.75 16649.53 16756.77 16785.98 16794.06 16796.40 16797.17 16797.49 Table X X I : Steady state integration : case 2 Chapter 5. The Viscoelastic 71 Problem isostatic solution. Table X X I gives results for the second case, again at intervals of lOOMa. These again show that the plate continues to deform u n t i l it is in isostatic adjustment w i t h the loading. We expect that the solution took longer to converge in the second case because it was further from isostatic adjustment at t = lOOMa from which time the forcing function was held in steady state. It should be noted that integrating to lOOOMa is somewhat unreasonable geophysically but the results obtained are certainly of interest since they would seem to indicate that under a given constant loading, any viscoelastic plate (regardless of length or flexural rigidity) will deform w i t h time u n t i l it reaches isostatic adjustment. 5.3.3 Geophysical Application F i n a l l y , we t u r n our attention to a forcing function which is more meaningful geophysically but for which there is no analytic solution available. We w i l l assume a Gaussian distribution for the loading. So, allowing f(t) to represent the time dependence and t a k i n g hx to be as represented i n equation (4.3.9) we have: (5.3.4) p = gh (r)f(t) P c T Ideally, we would like to begin with zero forcing and have the load increase to reach a steady state asymptotically. We therefore consider a function of the form: f(t) = - arctan ( — ) 7r which has the required properties that / ( 0 ) \ r r i T (5.3.5) J = 0 and f(t) —> 1 as t —• oo. application here is to the formation of sedimentary basins. The A basin of shape KT is assumed to be formed by tectonic subsidence. T h e basin initially has zero depth but increases w i t h time to reach a steady state asymptotically. T h e relaxation parameter, r , w i l l again take on the value 1 x 10 years so that the forcing grows by 50% i n m M a . 6 Chapter 5. The Viscoelastic 72 Problem m /(ioo) /(500) /(1000) 50 75 100 0.844 0.674 0.500 0.994 0.986 0.975 0.998 0.996 0.994 Table X X I I : T i m e dependence i n terms of m t(Ma) m = 50 m = 100 100 ,'.200 s* 300 •• 400 10102 13783 14975 15537 '15872 16092 16247 16358 16442 16505 5592 11433 13752 14802 15381 15743 15987 16158 16284 16378 •• 5 0 0 600 700 800 900 1000 : Table X X I I I : M a x i m u m deflection (metres) for case (i) In Table X X I I we give the value of f(t) at t = 100, 500 and lOOOMa for various values of m. We propose to integrate from 0 to lOOOMa which we hope to be a long enough time period for the forcing to have reached a reasonably steady state - appropriate choices of m will ensure this. We wish to choose values of m such that the following two conditions are met: (i) the forcing function is close to being i n steady state towards the end of integration, (ii) the rate of increase at the beginning of integration is not too great. We w i l l consider the following two cases: (i) L = 2.5 x 10 m, D = 1 x 1 0 N m , 5 (ii) L — 5.0 x 10 m, D = 1 x 10 5 25 23 N m and allow m to have the values 50 or 100 for both. Since, p = 0 at. t = 0 we have the i n i t i a l condition w = 0 at t = 0. T h e C . D . F . m e t h o d was used w i t h one V ( l , l ) cycle per time step w i t h At = l O M a . are presented i n Tables X X I I I and X X I V . T h e results Chapter 5. The Viscoelastic 73 Problem t(Ma) 100 200 300 400 500 600 700 800 900 1000 m = 50 m = 100 6688 9271 10719 11759 12539 13135 13598 13962 14254 14491 3565 7585 9579 10904 11873 12607 13173 13617 13969 14254 Table X X I V : M a x i m u m deflection (metres) for case (ii) T h e results show that for case (i) the deflection is close to being in isostatic adjustment towards the end of integration - this is as expected since the derivative, p, is close to zero at this stage and, as was seen i n section 5.3.2, the steady state solution is one representing an adjustment to isostatic equilibrium. We also see that increasing m reduces the i n i t i a l rate of increase in the loading but towards the end of integration there is little difference in the two cases - again this is as expected owing to the asymptotic behaviour of the forcing function. For case (ii) the deflection is not as close to being i n isostatic adjustment at the end of integration, we expect this to be due to the fact that the degree of support for this case is larger and therefore a longer period of integration will be required before the plate reaches isostatic adjustment. T h i s k i n d of behaviour has already been observed in the previous section. T h e effect of increasing m is the same as was observed for case (i)- Chapter 6 Conclusion 6.1 Results A derivation of the governing equation for deflection of a viscoelastic plate under an applied load has been presented. T h e equation was found to be fourth order (in space) and time-dependent. For the elastic problem, the time-dependence is dropped but the governing equation remains fourth order. B o u n d a r y conditions were derived from the assumption that the plate is simply supported at the edges so that there is zero deflection and zero curvature on the boundary. M u l t i g r i d algorithms for the numerical solution of the elastic and viscoelastic problems were presented. T h e discretization was applied after reformulating the fourth order equations i n terms of two coupled second order equations. T h e m u l t i g r i d method has proved to be a very fast and efficient solver for the elastic, problem. It has dealt effectively w i t h a wide range of plate lengths and flexural rigidities representing degrees of support ranging from 0.05% to 99.9%. T h e numerical solution has therefore proven to be stable w i t h respect to varying degrees of importance of the b i h a r m o n i c operator (which represents the highest derivative in the governing equation). In cases where the degree of support is very low and we have a singularly perturbed problem the numerical scheme i n fact converges quickly to the solution of the reduced equation. A s the degree of support increases and the biharmonic operator becomes more significant, the numerical solution takes longer to converge and becomes 74 Chapter 6. 75 Conclusion less accurate, however, this deterioration cannot be considered to be serious as the convergence rate and accuracy are still w i t h i n acceptable limits. A l l i n all we can be very pleased w i t h the performance of the m u l t i g r i d m e t h o d as applied to the elastic problem. T h e viscoelastic or time-dependent problem has proved to be the more challenging of the two. M a n y possible avenues may be travelled on b o t h the algorithmic level and on the level of choosing appropriate time discretizations. We examined two multigrid algorithms, L M G and M F M G , and compared results obtained from each of them against an analytic solution. A n estimate of the 0(2/i) discretization error was obtained from 2 the M F M G algorithm and compared to the exact change i n the derivative over the last time step. It was found that i n those cases where the discretization error was higher t h a n the change i n the derivative, L M G consistently provided better approximations to the solution and its derivative than M F M G - as expected. .Tfowever, M F M G did not consistently provide better results when the discretization error was less than the change i n the derivative. O n the whole, L M G provided better results at a lower cost and so the M F M G a l g o r i t h m was discarded. The time dependence was dealt w i t h i n three different ways: (i) B a c k w a r d E u l e r - a first order m e t h o d , (ii) A centred difference scheme which provides a second order approximation w i t h the same amount of work as B a c k w a r d E u l e r and (iii) a two-step B . D . F . scheme which again provides a second order approximation but at greater cost than the centred difference scheme. The two second order methods consistently provided better approximations than the first order B a c k w a r d E u l e r scheme. A comparison of the two second order methods d i d not reveal any conclusive results. However, it was decided, from the point of view of efficiency and ease of implementation to use the centred difference scheme for future numerical experiments. W e also examined the effect of m a i n t a i n i n g the forcing function at a steady state Chapter 76 6. Conclusion after a given period of time and observing the subsequent deflection. A n a l y t i c a l l y we expect the plate to relax, asymptotically, u n t i l it is in isostatic adjustment regardless of plate size or elastic properties - this behaviour was also borne out i n the numerical results. It appears that a very good degree of accuracy may be obtained w i t h only one V ( l , 1) L M G cycle per time step. T h i s represents a fast and highly efficient method. Perhaps the greatest advantage of the solution method we have developed lies i n its adaptability i n that a m i n i m u m amount.of work is involved i n solving the problem w i t h different forcing functions once the program has been set up. T h i s cannot be said of methods such as integral transform techniques and Fourier Series expansion methods which have been employed previously. T h e finite element technique also involves a lot of work in setting up (rather than solving) the problem when different forcing functions are considered. 6.2 Future Work It is fairly clear that most of our interest for future work should lie i n dealing w i t h timedependent m u l t i g r i d . For this application it was shown that M F M G did not represent a viable alternative to the very straightforward L M G algorithm. T h i s appears to be due to the high discretization error on the next-to-finest level relative to the change i n the solution from one time step to the next. However, further investigation into the applicability of M F M G against L M G is certainly required. Further investigation into the number of pre and post C G C relaxation iterations and the number of V-cycles per time step is also called for. We have consistently used one V ( l , 1) cycle per time step without any justification as to whether this is an o p t i m a l choice. T h e optimal number of V-cycles per time step is of course related to Chapter 6. Conclusion 77 the choice and control of the time step, A t , i n that integrating w i t h a time step A t using, for example, one V ( 2 , 2) cycle per time step involves the same amount of work as integrating w i t h a time step A t / 2 using one V ( l , 1) cycle per time step. T h i s question of a t t a i n i n g o p t i m a l efficiency at the algorithmic level is certainly worthy of further attention. In the geophysical context there is still room for improvement i n the development of forcing functions which are more realistic. A s regards the application to the formation of sedimentary basins, we have assumed that the density of the infilling sediment is constant. In fact, the body of sediment infilling the deflection may consist of several layers each having a different density so that for the elastic model we might have a forcing function of the form P = ^2Pi9hi(x,y) (6.2.1) i=i where p is the density of the z l t h layer having a thickness h,(x,y). A comparison of the elastic and viscoelastic models would also be of interest. T h e deflection of an elastic lithosphere w i t h time-dependent loading is governed by the equation: DV w(t) 4 +jw(t) = p{t) (6.2.2) H a v i n g solved the viscoelastic problem, the solution of (6.2.2) should not present any further difficulty numerically since it represents a simple continuation problem. T h e L M G algorithm w i t h one V ( l , 1) cycle per time step should prove successful provided p(t) does not change rapidly. A s regards the formation of sedimentary basins, a comparison of the stratigraphy predicted by the elastic and viscoelastic models for a given p(t) would be a valuable exercise since, to the author's knowledge, no definitive statement has yet been made as to which is the more realistic model. Chapter 6. Conclusion 78 A further necessary improvement if a true geophysical situation were to be modelled, is the extension of the m u l t i g r i d method to deal w i t h irregular domains. T h i s involves much computational work since the discretization, interpolation a n d restrict i o n operators become non-standard at points near the boundary. T h e grid coarsening process also involves further difficulty since a regular coarsening may not be possible near the boundaries. It is highly unlikely that grid coarsening could be continued u n t i l only one interior grid point remains thus m a k i n g it necessary to develop an exact solver for the problem on the coarsest grid, however, this is not a serious drawback. To conclude, we suggest that many more hours of research into b o t h the numerical and geophysical aspects of this topic are in order. We content ourselves w i t h having provided what we hope to be a sound basis for future investigation i n this field of study. Bibliography [l] Ascher, U . conductor On numerical device [2] B e a u m o n t , C . differential algebraic problems with application to semi- S . I . A . M . J o u r n a l of N u m e r i c a l Analysis (to appear). simulation The evolution of Sedimentary .Basins on a Viscoelastic Lithosphere Geophysical J o u r n a l of the R o y a l A s t r o n o m i c a l Society 55 p.471, 1978. [3] Drucker, D . C . Introduction to the mechanics of deformable bodies McGraw-Hill 1967. [4] G e n d l e r , E . Multigrid methods for time-dependent parabolic equations Master's Thesis, W e i z m a n n Institute of Science, Rehovot, Israel 1986. [5] Hackbusch, W . [6] Hunter, S . C . Multi-Grid Mechanics methods of Continuous [7] N u n n , J . A . , a n d N . H . Sleep ' Basins: a three and Applications dimensional Media Thermal study Springer Verlag 1985. E l l i s H o r w o o d 1983. Contraction and Flexure of the Michigan Basin of Iniracratonal Geophysical J o u r n a l of the R o y a l A s t r o n o m i c a l Society 76 p.587, 1984. [8] N u n n , J . A . , N . H . Sleep a n d W . E . M o o r e hydrocarbons in the Michigan Basin Thermal Subsidence Flexure Advances i n Geophysics, 21 1979. [10] Turcotte, D . L . a n d G . Schubert of A m e r i c a n Association of Petroleum Geologists B u l l e t i n 68 N o . 3 p.296, 1984. [9] Turcot te, D . L . and generation Geodynamics 79 J . W i l e y 1982. 80 Bibliography [11] V i n s o n , J . R . Structural Mechanics [12] W a t t s , K a r n e r and Steckler Basins Philosophical J.Wiley Lithospheric 1974. Flexure and the evolution of Sedimentary Transcripts of the R o y a l Society of L o n d o n A 3 0 5 p.249, 1982. [13] Wesseling, P. in Multigrid Methods M a t h e m a t i c s , S . I . A . M . 1987. S.F. M c C o r m i c k ed., Frontiers i n A p p l i e d
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A multigrid method for determining the deflection of...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A multigrid method for determining the deflection of lithospheric plates Carter, Paul M. 1988
pdf
Page Metadata
Item Metadata
Title | A multigrid method for determining the deflection of lithospheric plates |
Creator |
Carter, Paul M. |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | Various models are currently in existence for determining the deflection of lithospheric plates under an applied transverse load. The most popular models treat lithospheric plates as thin elastic or thin viscoelastic plates. The equations governing the deflection of such plates have been solved successfully in two dimensions using integral transform techniques. Three dimensional models have been solved using Fourier Series expansions assuming a sinusoidal variation for the load and deflection. In the engineering context, the finite element technique has also been employed. The current aim, however, is to develop an efficient solver for the three dimensional elastic and viscoelastic problems using finite difference techniques. A variety of loading functions may therefore be considered with minimum work involved in obtaining a solution for different forcing functions once the main program has been developed. The proposed method would therefore provide a valuable technique for assessing new models for the loading of lithospheric plates as well as a useful educational tool for use in geophysics laboratories. The multigrid method, which has proved to be a fast, efficient solver for elliptic partial differential equations, is examined as the basis for a solver of both the elastic and viscoelastic problems. The viscoelastic problem, being explicitly time-dependent, is the more challenging of the two and will receive particular attention. Multigrid proves to be a very effective method applicable to the solution of both the elastic and viscoelastic problems. |
Subject |
Multigrid methods (Numerical analysis) Viscoelasticity Elastic plates and shells |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079633 |
URI | http://hdl.handle.net/2429/27854 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1988_A6_7 C37_7.pdf [ 4.3MB ]
- Metadata
- JSON: 831-1.0079633.json
- JSON-LD: 831-1.0079633-ld.json
- RDF/XML (Pretty): 831-1.0079633-rdf.xml
- RDF/JSON: 831-1.0079633-rdf.json
- Turtle: 831-1.0079633-turtle.txt
- N-Triples: 831-1.0079633-rdf-ntriples.txt
- Original Record: 831-1.0079633-source.json
- Full Text
- 831-1.0079633-fulltext.txt
- Citation
- 831-1.0079633.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079633/manifest