A M U L T I G R I D M E T H O D F O R D E T E R M I N I N G T H E D E F L E C T I O N O F L I T H O S P H E R I C P L A T E S B y Paul M . Carter B . Sc., The University of Sheffield, 1985 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SC IENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES INSTITUTE OF APPL IED MATHEMAT ICS W e - accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA June 1988 © Paul M . Carter, 1988 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of MATHEMATICS The University of British Columbia Vancouver, Canada Date 24 th JUNE 1988 DE-6 (2788) Abs t r ac t 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 wi th min imum 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 mult igr id method, which has proved to be a fast, efficient solver for ell iptic part ial differential equations, is examined as the basis for a solver of both the elastic and viscoelastic problems. The viscoelastic problem, being explici t ly time-dependent, is the more challenging of the two and wi l l receive particular attention. Mu l t i g r i d proves to be a very effective method applicable to the solution of both 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 Elast ic i ty and Viscoelastic Theory 10 2.3 T h i n Elast ic Plates 13 2.4 Transverse Loading 19 2.5 Nondimensionalization 20 2.6 Boundary and Init ial Condit ions 21 2.7 Reformulation of the governing equations 22 3 The Multigrid Method 25 3.1 The basic method 25 3.1.1 Relaxat ion 28 3.1.2 Restr ict ion 31 3.1.3 Prolongation 31 ii i 3.1.4 The Coarse G r i d Operator 32 3.2 Basic Mul t i g r id Algor i thms 32 3.3 Mul 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 4 The Elas t ic P r o b l e m 41 4.1 Ana ly t i c Solution 41 4.2 Numerical Solution 46 4.3 Numerical Results 48 4.3.1 Periodic Loading 48 4.3.2 Mul t i g r i d vs. Gauss-Seidel 52 4.3.3 Gaussian Loading 53 5 The Viscoelas t ic P r o b l e m 56 5.1 Ana ly t i c Solution 56 5.2 Numerical Solution 58 5.3 Numerical Results 61 5.3.1 Periodic Loading 61 5.3.2 Integration wi th steady forcing 69 5.3.3 Geophysical Appl ica t ion 71 6 Conc lus ion 74 6.1 Results 74 6.2 Future Work 76 Bib l iog raphy 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 10 5 m, Periodic Loading 50 I V L = 5.0 x 10 5 m, Periodic Loading 50 V L = 1.0 x 10 6 m, Periodic Loading 50 V I M u l t i g r i d vs. Gauss-Seidel 52 V I I L = 2.5 x 10 5 m , Gaussian Loading 54 VI I I L — 5.0 x 10 5 m , Gaussian Loading 54 I X L = 1.0 x 10 6 m , Gaussian Loading 54 X Degree of support (Viscoelastic Model) 58 X I Backward Euler, L M G vs. M F M G 63 X I I C . D . F . method, L M G vs. M F M G 63 X I I I Two-step B . D . F . , L M G vs. M F M G 63 X I V Backward Euler, A t = l O M a -. . 65 X V C . D . F . method, A t = l O M a ". . 65 X V I Two-step B . D . F . , A t = l O M a 65 X V I I Backward Euler, At = 5 M a 67 X V I I I C . D . F . method, A< = 5 M a 67 X I X Two-step B . D . F . , A t = 5 M a 67 X X Steady state integration : case 1 70 X X I Steady state integration : case 2 70 v X X I I T ime dependence in terms of m X X I I I M a x i m u m deflection (metres) for case (i) X X I V M a x i m u m deflection (metres) for case (ii) v i L i s t o f F i g u r e s 2.1 Elements of Stress 11 2.2 T h i n plate deformation 14 2.3 Hydrostat ic Restoring Force 19 3.4 The mult igr id 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 Dr . U r i Ascher and Dr . Ma t t Yedl in during the research and wri t ing of this thesis. vm C h a p t e r 1 In t roduc t ion The aim of this dissertation is to provide a fast, efficient, numerical method for deter-mining the deflection of lithospheric plates under an applied load. In particular, the mul t igr id method is employed in this capacity. We wi l l examine plates having hori-zontal dimensions ranging from those of sedimentary basins to continental plates. We wi l l also consider a wide range of other physical parameters. Therefore, our atten-t ion 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 mainly on the numerical aspects of the topic while still maintaining a geophysical framework wi th in which our forcing functions and model parameters wi l l lie. Our discussion wi l l begin wi th a justification of the treatment of lithospheric plates as thin 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 thin plate deflection to which much attention has been given in the engineering literature. The equations governing thin 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. The application of the current work to the evolution of sedimentary basins is of particular interest since such areas tend to be rich in 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 in i t ia l deflection in the plate and incidentally, the thermal anomaly causing this contraction also plays an important role in the maturat ion of petroleum, see Nunn , Sleep and Moore [8]. Various models for the evolution of sedimentary basins have been developed for both elastic and viscoelastic plates. Beaumont [2] examined the formation of sedimen-tary basins due to the loading of elastic and viscoelastic lithospheres. The problem was considered in cyl indrical co-ordinates wi th 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 wi th the heavy-side-step function to represent the time dependence. Model l ing of the plate loading was therefore very simple, however, math-ematical tractabili ty was maintained in the process. B o t h the elastic and viscoelastic problems were solved using integral transform techniques. Beaumont also addressed the problem of basin ini t iat ion, i.e. that mechanism by which an in i t ia l deflection is created in the lithosphere in which sediment subsequently accumulates. Six possible mechanisms were proposed including 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. Introduction 3 instantaneously formed graben was examined. The 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 in best agreement wi th seismic data obtained from the Nor th Sea basin. Watts , Karner & Steckler [12] examined a two-dimensional model where thermal contraction was assumed to control the tectonic subsidence of the plate, i.e. that subsi-dence which is independent of subsequent deflection due to sedimentary infil l . 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. Two models were considered for tectonic subsidence the latter being attr ibuted to thermal contraction in both cases. The first assumed that the depth of the basin increases as the square root of t ime since basin ini t ia t ion. Aga in , formation of the basin was assumed to occur over a period of lOOMa. The second model was based on the assumption that the lithosphere is stretched laterally at the time of basin ini t ia t ion resulting in plate thinning and heating which in turn resul t / in plate subsidence on cooling. A variation in elastic thickness, and therefore flexural rigidity, w i th time was allowed for in both the elastic and viscoelastic. models between time periods. Spatial variation of the flexural rigidity was also considered. The 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 in two-dimensions, the extension to three dimensions being non-tr ivial . It was concluded that the best model fit to overall basin geometry arises from an elastic model in which the flexural rigidity increases with Chapter 1. Introduction 4 t ime after basin ini t iat ion. This increase is due to an increase in elastic thickness wi th age as the plate cools. It should be noted that the conclusion here is at variance wi th 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 in the formation of sedimentary basins are thermal contraction, which controls the overall shape of the basin, and sedimentary loading. Nunn & Sleep [7] considered the load due to both thermal contraction and sedimen-tary loading for an elastic and viscoelastic lithosphere. A sinusoidal spatial variation was assumed for both 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 in the model. A n estimate of the thickness of sedimentary layers obtained from lithostratigraphic data for the Michigan Basin was used to compute the Fourier coefficients for the ex-pansion representing the load due to sedimentary infi l l . 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, tmax, after basin ini t ia t ion. These coefficients were therefore re-computed for each set of rheological parameters. It was estimated that the load due to sedimentary infill constitutes approximately 75% of the total driving force. The best model fit to observations was obtained for a viscoelastic lithosphere wi th low flexural rigidity. However, it was pointed out that wi th the data available at present, either model yields acceptable results, and a better un-derstanding of basin ini t iat ion, sediment budget and sediment compaction is necessary before any definitive statements can be made. We wi l l examine the deflection due to the loading of both elastic and viscoelastic Chapter 1. Introduction 5 plates, al l models being three dimensional. Simple models which account only for sedimentary loading wi 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 in a continuous fashion. The infilling sediment is assumed to have constant density. We now present an overview of material to be found in the following chapters. In Chapter 2 we wi l l call upon material found in the engineering literature to provide a derivation of the governing equations for the deflection of a viscoelastic plate. In doing so we wi l l closely follow a derivation of the equation governing deflection of an elastic plate replacing the elastic constitutive equations wi th those for a viscoelastic solid. We wi l l also present appropriate boundary and ini t ia l conditions along wi th a nondimensionalized form of the governing equations. To conclude the chapter we wi l l consider an appropriate reformulation of the governing equations which wi l l make their solution easier numerically. For both the elastic and viscoelastic problems, we are faced wi th 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 mult igrid method and consider its application to the two problems in hand. The presentation includes a discussion of the main components of mult igr id as well as some basic algorithms. The 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. Chapter 3 is therefore concluded wi th a discussion of mult igrid for time-dependent problems. Two algorithms are considered, one of which follows immediately from the time-independent, elastic problem. The applicabili ty of the two algorithms is discussed and some attempt is made to ascertain the conditions under which one is expected Chapter 1. Introduction 6 to perform more efficiently than the other. The performance of the two algorithms is tested numerically in Chapter 5. In Chapter 4 we examine the elastic problem on a square domain. The assumption of a square domain places some restriction on the problem geophysically, but it s t i l l maintains a definite advantage over the assumption of two-dimensionality which has been used by other authors. However, the regular domain 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 with the presentation of an analytic solution to which later numerical results wi l l be compared. This 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 wi l l be derived. It wi 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 mul t igr id solution wi l l be examined for various plate dimensions and elastic constants representing a wide spectrum of degrees of support. A comparison between mul t igr id and a simple Gauss-Seidel iteration wi l l be carried out showing that mul t igr id performs relatively more efficiently except in cases where the degree of support is very low in which case the numerical solution becomes t r iv ia 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 with the presentation of an analytic solution representing an exponentially increasing forcing function wi th periodic spatial dependence. This representation has some major shortcomings in the Chapter 1. Introduction 7 geophysical context but it is nonetheless considered acceptable in so far as it provides a means for assessing the numerical results more objectively. Two mult igr id algorithms for time-dependent problems wi l l be contrasted and compared. After an examination of the results and a consideration of the work involved in each of the algorithms, one wi l l be eliminated on the grounds of both 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. Having chosen a suitable algorithm we then examine the performance of three dif-ferent approaches to dealing wi th the time discretization, namely, Backward Euler , a centred difference scheme and a two-step backward difference formula. The first of these provides a first order approximation in time to the solution while the other two provide second order approximations. It wi l l be seen that the centred difference scheme, while providing a second order approximation, involves no more work or storage than Back-ward Euler, the first order method, making it an extremely attractive method from the point of view of both 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 in future numerical experiments. Next we consider the effect of allowing the forcing to reach a steady state. Ana ly t -ically we expect all plates, regardless of size or elastic properties, to relax unt i l they are in hydrostatic equil ibrium wi th 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 equil ibr ium. Final ly , we consider a forcing function which is more reasonable geophysically. The spatial variation is represented by a Gaussian distr ibution and the time dependence is Chapter 1. Introduction 8 such that forcing is ini t ia l ly zero and reaches a steady state asymptotically. T w o sizes of plate wi th different elastic properties wi 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 equil ibrium. Conclusions arising from the numerical experiments are presented in Chapter 6. Suggestions as to the direction of future research are also proposed. C h a p t e r 2 D e r i v a t i o n o f t h e G o v e r n i n g E q u a t i o n 2.1 L i t h o s p h e r i c P l a t e s In considering the deflection of lithospheric plates, we wi l l be encountering horizontal dimensions ranging from 250km to 1000km. This covers anything from the formation of sedimentary basins to the deflection of continental plates due to topographic load-ing. In order to derive an equation governing plate deflection we make the following observations. Firs t , that the horizontal extent of the plate is not so large as to make the effects of the earth's curvature of importance. We therefore assume that in their undisturbed state, the plates are completely flat. Second, the thickness of the plate is very small in 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. This is considerably less than the seismic thickness which is of the order of 100km, see Nunn & Sleep [7, p.590]. It is the mechanical thick-ness which is of importance in 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 Equation 10 wi th an adequate model for achieving this. However, if we wish to study the evolution of the plate's deflection then a viscoelastic model is more appropriate. 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 t ime goes on. Obviously, from a geophysical standpoint, we wi l l not be considering the complete removal of the load but, in 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 Having justified the treatment of lithospheric plates as th in elastic or viscoelastic plates we now derive an expression governing the deflection of thin 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 in Figure 2.1 with Cartesian axes-in the given orientation. This volume element is subject to normal and shear stresses due to its interaction wi th neighbouring volume elements. Figure 2.1 indicates the direction in which these stresses are taken as positive. Note that stresses acting on surfaces normal to the x 2 -di rect ion 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. We adopt a suffix notation wi th respect to Cartesian axes (x1, x2, x3) and assume the Einstein Summation Convention to be in operation over repeated suffices. The notation (. )_j Chapter 2. Derivation of the Governing Equation 11 ^33 Figure 2.1: Elements of Stress represents differentiation wi th respect to xr Following Vinson [11], we now present the governing equations for the deformation of an elastic, body. A more in-depth discussion than that given in V inson may be found in Hunter [6]. Consideration of conservation of momentum for a static, elastic body in the absence of body forces (such as gravity) yields the following equil ibr ium equations °~H,j = 0 (2-2.1) Conservation of angular momentum requires that a,? = on (2.2.2) in other words that the stress tensor be symmetric. A deformed elastic body has associated with it not only stresses but also strains and a displacement field. Assuming small displacements, it can be shown that the strains e,7 (suffix convention being the same as for the stresses at]) are related to the Chapter 2. Derivation of the Governing Equation 12 displacement field in the following way £H = \(u*J + U3,i) (2.2.3) F ina l ly we have the constitutive or strain-stress relations. Assuming the material to be isothermal, isotropic and homogeneous we have: 1 £ij E (1 + v)uij - v ukk 8{ (2.2.4) where the superscript e refers to an elastic solid, E = Young's Modulus , v = Poisson's Rat io and 6{j = Kronecker Del ta 1 if % = j, • 0 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 wi th the replacement of the strains e\- by the rate of strain zv- where the dot denotes differentiation wi th respect to time and the superscript v refers to a viscous solid. For an isotropic material we replace Young's Modulus E, by the coefficient of viscosity, p. Hence, for a linear viscous solid we have: £ii = ~~ K1 + v) ~ v t>ij} (2.2.5) A4 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 Maxwel l model and the K e l v i n model. The Maxwel 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. The 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. As pointed out by Turcotte [9, p.70], the usual approach, at least in the geophysical context, is to deal wi th the problem using the Maxwel l model. Chapter 2. Derivation of the Governing Equation 13 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 wi l l be of use to have a form of the constitu-tive equations where summation from only 1 to 2 and not to 3 is implied over repeated suffices. Inverting equation (2.2.6) wi th this in mind we have: <7„ + o~ij _ E T (1 + 1/ ^ + (1 - ") - kk^ij where T _ E [2.2.7) (2.2.8) is known as the relaxation constant. 2.3 T h i n Elas t ic Plates Vinson [ l l ] derives an equation for the deflection of a thin elastic plate. We shall follow along the same lines using the constitutive equations for a viscoelastic solid in place of those for an elastic solid, hence deriving an equation governing the deflection of a thin viscoelastic. plate. Consider a thin, rectangular, isotropic plate with horizontal dimensions x\, x\ and thickness h. We assume that the conditions h -C x\ and h <C x% define a thin plate. We now position a set of Cartesian axes at the centroid of the plate as shown in Figure 2.2 so that the mid-plane of the plate is the plane x3 = 0. Consider a cuboid element of the plate extending through the plate thickness as shown. We assume that on the application of a transverse load this element remains Chapter 2. Derivation of the Governing Equation 14 C U B O I D E L E M E N T M I D - P L A N E 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 wi th respect to the co-ordinate system. The element therefore remains straight upon load application. This means that sides of the cuboid init ia l ly ly ing in the X ] and x 2 -directions remain perpendicular to sides ly ing in the x 3 -direct ion after deformation. Hunter [6, p . I l l ] interprets the strains etJ (i ^ j) as representing a measure of the deviation from perpendicularity, after deformation, of sides init ia l ly pointing i n the i and j-direct ions. Since in the current problem sides in the X j and x 2 -directions remain perpendicular to sides in the x 3 -direct ion, we may deduce that e 1 3 — £ 2 3 = 0-We also assume that an applied transverse load results in plate bending but does not lead to compression or tension in the plate thickness direction. This is in contrast to the behaviour of a sponge where a transverse load is absorbed by compression in 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 Equation 15 the elongation or compression of an element in the z-direction. Since in the current application there is no elongation or compression in the ^ - d i r e c t i o n we therefore deduce that e 3 3 = 0. Equat ion (2.2.3) then implies that u3 (the displacement in the x3-direction) is independent of x3. Hence, the admissible form of the subsequent displacement field for a thin elastic plate is: ux = u1(x1,x2,x3,t) = u<t(x1,x2,t) + x3a1(x1,x2,t) u2 = u2(x1,x2,x3,t) = u°2(x1,x2,t) + x3a2(xux2,t) > (2.3.1) u3 = u3(x1,x2,f) where u°(xi,x2,t) — u , ( x i , x 2 , 0, t) for i = l , 2 and a1(xi,x2,t) represent rotations to be defined later. From here we shall adopt modified forms of the suffix notation and summation convention which have been used to this point. We shall assume that al l suffices run only from 1 to 2 and not to 3. A n y variable wi th a suffix 3 wi l l be writ ten explicitly. We now derive expressions for the net moments and forces acting on an element across the plate thickness. The moments, per unit area about the mid-plane are given by: = cr%jX3 Integrating these across the plate thickness we have: h Mij — J\ azjx3 dx3 (2.3.2) The transverse shears, ql, in the x 3 -direct ion are given by: and so integrating over the plate thickness we have: h Q% = f\ al3 dx3 (2.3.3) Chapter 2. Derivation of the Governing Equation 16 We now derive an integrated form of the equil ibrium equations (2.2.1). Mu l t i p ly ing (2.2.1) by x3 and integrating over the plate thickness we have 0 - ^ 2 : 3 dx3 + / c r l 3 ,3X 3 dx3 = 0 Mijtj + Ol3X3 o%3 dx3 = 0 =*> Mijj - Q, .= 0 (2.3.4) assuming there are no shear stresses acting on the surface of the plate. Integrating (2.2.1) for t = 3 we have 2 12 h o3h] dx3 + I 0-33,3 dx3 = 0 Qj,i + o-33 =4- Qj,j +p = 0 (2.3.5) where p = p(x1,x2,t) = o-33 represents the net force per unit area acting normal to, the plate surface. Combining equations (2.3.4) and (2.3.5) we have Mijjti + p = 0 (2.3.6) We now relate the moments Ml3 to the vertical displacement of the plate. Combining equations (2.2.3) and (2.2.7) we have Chapter 2. Derivation of the Governing Equation 17 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 making use of (2.3.1) we have: h 2 2=3 E [l + u) / • 0 , - 0 \ , V • 0 r + E (1 + 1 -[Ct .2 ( 1 - ^ ) We now relate the rotations a%(xi,x2,t) to the vertical deflection u3(x1,x2,t) which from here wi l l be denoted by w = w(x1, x2,t). Equations (2.2.3) and the assumption that e 1 3 = e23 = 0 imply: a h l ) X3 + 2 c dx-. dx-- (2.3.7) £23 = - ' dxi dw 0 =±> 0 d w a, 9xi dw dx2J dx2 Notice we have dropped the comma notation for derivatives and wi l l continue to do so for the remainder of this discussion. Equat ion (2.3.7) may now be writ ten Mi. -D 1 - V d2w d2w dxjdx., + v dxf.dxk ZJ (2.3.8) where D , the flexural r igidity of the plate, is given by D = E h 3 12(1 - v2) Equations (2.3.6) and (2.3.8) represent a closed system of equations for the vertical deflection w, in terms of the applied transverse load p. E l imina t ing Ml3 from these equations we have dx\ I D d2w | d2w dx2\ + 2 dxidx2 D(l -d2 w dx-t dx* Chapter 2. Derivation of the Governing Equation 18 dx\ d2w d2w .dx\ dx\ which may be written p , • = - +P T V where and V2{D(V2w)) - (1 - u)0\D,w) = -+p T (2.3.9) 04(D,w) = d2Dd: w dx\ dx\ d2D d2w d2D d2w + dx\ dx\ dx1dx2 dxidx2 dx\ dx\ Recal 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 maintain some time dependence if h changes with time) and so we have: D V 4 P w — hp r (2.3.10) wnere V \ i ; = V2(V2w) Inverting equations (2.2.4) and using the resulting equations in place of (2.2.7) in the above derivation leads to an equation relating the vertical deflection of an elastic, rather than viscoelastic, plate to transverse loading. Vinson [11, p.16] shows that the equation thus obtained is DV4w=p (2.3.11) instead of (2.3.10). Chapter 2. Derivation of the Governing Equation 19 u Surface Rock Pa Pc w h. Mantle •m B F l u i d Mant le A Pm Figure 2.3: Hydrostatic Restoring Force 2.4 T r a n s v e r s e L o a d i n g When 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, wi l l infill the cavity thus formed, resulting in a net hydrostatic restoring force. Following Turcotte & Schubert [10, p.121] we consider the case shown in figure 2.3 where the continental lithosphere of density pm and thickness hm rides on the fluid mantle, also of density pm. Overlying the lithosphere is a layer of sediment of thickness hc and density pc. A n applied transverse load pa causes a deflection w in the lithospheric plate. We shall now calculate the hydrostatic restoring force and hence the net transverse load. The pressure acting at the point A in Figure 2.3 below the deflected plate is given by Chapter 2. Derivation of the Governing Equation 20 whereas the pressure at the point B under an undisturbed portion of the plate is given by: p c g h c + p m g ( h m + w) (2.4.2) Subtracting (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 = p a + {pc - p m ) g W = p a - JW (2.4.4) where 7 = ( p m — p c ) g . Hence, the governing equation for deflection of a viscoelastic plate under an applied load p a may now be written: DVAw + 7 ^ - + iu) = — +pa (2.4.5) and for an elastic plate we have: DV4w + j w = p a (2.4.6) From hereon, the subscript 'a ' on the loading, p , wi l l be dropped so that p w i l l be understood to mean the applied load. 2.5 N o n d i m e n s i o n a l i z a t i o n For the purposes of numerical computation we.wil l now express equations (2.4.5) and (2.4.6) in nondimensional form. As our solution domains wi l l be square, we scale all lengths by a factor L . A l l times wi l l be scaled by a factor T. Using primes to denote nondimensional parameters, we have: w = w'L, p = p'/L, D = D'L2 Chapter 2. Derivation of the Governing Equation 21 7 = 7'/L 2 , T — T'T d I d d I d d i d dx! Ldx\\ dx2 Ldx'2' dt T dt' Hence, equation (2.2.5) becomes Cancel l ing the common factor 1/LT from each of the terms and dropping the primes we have: 7 4 - . , (w , -'\ P DX7qw + 7 ^ - +wj= - + p (2.5.1) where now, all variables and parameters are nondimensionalized wi th respect to a length scale L and a time scale T . Similarly, equation (2.4.6) remains essentially unaltered in nondimensional form. 2.6 B o u n d a r y and In i t i a l Condi t ions Consider a solution domain Q, wi th boundary d£l. Following Nunn & 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, <9fi. Hence, .w = 0, V2u> = 0 on dn ' (2.6.1) For the time-dependent (viscoelastic) problem, we also need some ini t ia l condition. A n y condition consistent wi th (2.6.1) at t = 0 wi l l suffice. In cases where the exact solution, we, is known we wi l l take: w = we{t = 0) on Q, (2.6.2) as our in i t ia l condition where it is assumed that we satisfies (2.6.1). In cases where the exact solution is not known, we follow Nunn & Sleep [7, p.630] in assuming that the Chapter 2. Derivation of the Governing Equation 22 deflection is ini t ia l ly 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 al l time they may be wri t ten equivalently as: w = 0, V2w = 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) wi th the above boundary and in i t ia l conditions. 2.7 Reformulation of the governing equations Part icular analytic solutions are readily available for certain choices of the forcing function. Examples of these wi l l be encountered later where an analytic solution wi l l be used to assess the accuracy of a numerical solution. Nunn & Sleep [7] provide an in-depth discussion of the solution for periodic loading using the elastic and viscoelastic models. The current aim however, is to provide a means for solving the elastic or viscoelastic problem numerically in as efficient a manner as possible. We propose to solve the equations using finite difference discretizations, however, as both problems have 4th order derivatives in space, we are presented wi th several computational difficulties if they are left in their current formulation. We wi 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. Equat ion Chapter 2. Derivation of the Governing Equation 23 (2.4.5) for the viscoelastic problem becomes: DV2v + 7 ^ — + uj = - + p V2u = v u = w (2.7.1) wi th boundary conditions: u = 0, v = 0 on dft (2.7.2) and ini t ia l condition: w(t = 0) = 0 on 9, (2.7.3) assuming the exact solution is unknown. Equat ion (2.4.6) for the elastic problem may (2.7.4) be written: D V 2 u + jw = p V2u> = v with boundary conditions: w = 0, v = 0 on dn (2.7.5) Notice that the boundary and ini t ia l conditions arise from a reformulation of those presented in section 2.6. It is clear that the elastic, rather than viscoelastic, problem is the easier of the two to solve. The solution of the elastic problem requires an efficient solver for coupled elliptic p.d.e.'s of the form: £ 1(u 1,« 2) = / 1 C2(u\u2) = 0 where u1 = u2 = 0 on dn (2.7.7) Later it wi l l be seen that having chosen a discretization for the time derivative, the viscoelastic problem may, in part, be written in the same form. So, once an efficient on n (2.7.6) 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. C h a p t e r 3 T h e M u l t i g r i d M e t h o d 3.1 The basic method 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. Then, u^, the exact solution of the discrete ^problem, satisfies: where C]x is a. discretization of the differential operator C. and fh is a discrete approx-imation to / . The 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. This fill-in is part icularly taxing as regards storage ^requirements. More efficient methods, in terms of storage, may be found amongst simple iterative schemes such as Gauss-Seidel. However, the asymptotic convergence of such methods tends to be slow and • Chuh = fh on Q (3.1.3) and Uh = 0 on dn 25 Chapter 3. The Multigrid Method 26 for very large systems such methods are extremely slow and therefore expensive. The mul t igr id method, about to be described, proves to be a fast, efficient solver for such a system. The mult igrid method can be viewed as a defect-correction method. Given an approximation Uh to Uh, the solution of (3.1.3), we define a defect, d^, by dh = h- £hUh (3.1.5) and a correction Vh satisfying: Ch{Uh + vh) = dh + ChUh (3.1.6) which, in the case of a linear operator Ch, reduces to Chvh = dh (3.1.7) where vh = 0 on dfl (3.1.8) Hence, our original problem (3.1.3) can be writ ten in terms of a defect, dh. and a correction, Vh- Having computed dh and solved exactly for Vh, our corrected solution becomes uh:=Ui + vh (3.1.9) so that uh now satisfies (3.1.3) exactly. As things are right now, we stand to gain nothing from our reformulation of the problem in 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 mult igridding. In the mult igr id method, we propose to construct a hierarchy of successively coarser grids on which the defect for the problem on the previous finer grid may be represented. For simplicity, in this ini t ia l discussion we shall confine ourselves to only two grids, a fine grid 0^ and a coarse one £7H, where, for example, H = 2h and C fit,,. Chapter 3. The Multigrid Method 27 R E L A X 1 C-h.Uk = fh # 4 R E S T R I C T S O L V E C-HVH p R O L O N G A T •E R E L A X ^ 2 £h.Uh = fh Figure 3.4: The mult igr id two-level cycle W i t h a coarser grid defined, it is proposed to solve (3 .1.6) on this coarser grid rather than the finer one as this w i l l , of course, be computationally 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. The elements of a mul t igr id two-level cycle are presented in Figure 3.4 where it has been assumed that is a linear operator. A s can be seen, the method consists of four main elements: 1. - Relax": generally consists of v sweeps of a simple iterative method such as Gauss-Seidel 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 Method 28 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 in more detail we must decide on a grid coars-ening process. We wi l l restrict ourselves to square domains on which we have a uniform grid. Assuming that the finest grid contains an odd number of grid points in each di-rection then a coarser grid can be defined by doubling the grid spacing in each direction so that al l points on the coarse grid are common to the fine grid and the grid remains uniform. Many other grid coarsening processes exist but the one described above is the easiest, to implement and wil l suit our purposes adequately. W i t h this assumed, we wi l l now consider the main components of the mul t igr id cycle in more detail. 3.1.1 Relaxation The purpose of the relaxation routine is to smooth out high frequency components of the error: eh = Uh-uh (3.1.10) The generation of other error components is also undesirable and so the relaxation method, in 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 = Ck, x = uh, b= fh so we are interested in solving the system Ax — b. Let A be an approximation to A such that B = A~x is easily computed. We may now write an iterative scheme for Chapter 3. The Multigrid Method 29 solving the given system as follows: n O + i ) : = (I - BA)xW + Bb (3.1.11; The matr ix M := I — BA is called the amplification matr ix. The methods we are about to describe are determined by the choice of A. We begin by considering a split t ing of the matr ix A into three components: A = D — E (3.1.12) 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'1(E + F) so that pointwise the iteration may b.e written: .(j'+i) 1 Q>kk X] akl- Si) 1=1 i^k (3.1.13) 2. The Gauss-Seidel Method where we choose A and the pointwise iteration may be written: D - E so that M = (D- E)~lF 1 o-kk Yl a k i x i i<k (j+i) l>k 3) (3.1.141 B o t h the Jacobi and Gauss-Seidel iterations can be shown to converge provided the matr ix 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 matr ix wi th m < n. The kernel of this matr ix, ker(l jf) , is therefore non-tr ivial . Hence, if dn •= iff dh represents Chapter 3. The Multigrid Method 30 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. This has disastrous consequences when solving the equation Cjjvjj = dn since if djj — 0 then VJJ = 0 and the iteration wi l l 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 wi th 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 wi th 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: Cheh = dh (3.1.15) Since C-h arises from the discretization of a differential operator, it does not, as noted by 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 high frequency components of dh which as we have seen are undesirable. We therefore require some means of smoothing out high frequency components of eh-It 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 in this capacity. However, it should be rioted that variations of Jacobi, i.e. damped Jacobi, can be used as effective smoothers wi th an appropriate choice of the damping 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, Al ternat ing Direction Implicit and Conjugate Gradient but we wi l l not discuss these methods here as Gauss-Seidel wi l l prove to be adequate for our needs. Chapter 3. The Multigrid Method 31 3.1.2 R e s t r i c t i o n The restriction procedure serves to transfer the defect from the fine to the coarser grid. The two most popular methods are: 1. Straight Injection where the defect is injected through to the coarser grid at points common to both 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 method may be written: (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 both grids. 3.1.3 P ro longa t ion The purpose of the prolongation routine is to transfer the correction from the coarse grid to the finer one. This may be achieved by injecting the corrections up to the finer grid at points common to both grids and then performing a bilinear interpolation. So, in stencil form we have: (3.1.17) Higher order interpolation may of course be carried out but, for our purposes, bilinear interpolation wi l l be sufficient. IH — 16 1 2 1 2 4 2 1 2 1 1 4 1 2 1 2 4 2 1 -2 1 Chapter 3. The Multigrid Method 32 3.1.4 T h e Coarse G r i d Opera to r Fina l ly we must define the coarse grid operator CH- There are two methods available to us for achieving this. 1. We 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 CH-2. The Gale rk in method defines the coarse grid operator in terms of the fine grid operator and the prolongation and restriction operators as follows: CH •= I K C J H • (3.1.18) The first method is the most popular for use with finite difference discretizations owing to its shear simplicity. However, for finite element methods the second method is the natural choice. It should also be noted that in 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 wi th finite differences and so the first of these options wi l l be employed. 3.2 Bas ic M u l t i g r i d A l g o r i t h m s Extending the concept of the two-grid method outlined in 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 n — 1 interior grid points in each direction then we could continue coarsening the grids unti l the coarsest grid contains only 1 interior grid point provided that the correct physical problem is st i l l represented on the coarsest grid. We assume this to be the case so that the exact Chapter 3. The Multigrid Method 33 solution of the discrete problem on the coarsest grid is easily obtained by one sweep of Gauss-Seidel. Even though this is generally considered to be an iterative method, when used in this capacity wi th only one grid point where the solution is unknown, it in 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 problem on that grid. Assuming 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 notation we may now define a basic mul t igr id algorithm, L M G , for the problem Cpup = fp as-follows: Algorithm 3.1 (LMG) Procedure L M G ( / ,7 , ui, u2)\ Begin if / = 1 then S O L V E ^ U i = e^ ) else begin for i:= 1 to ux R E L A X ( £ ^ = d\); RESTRICT; {d^x <- j / - 1 ( d j - £ / U j ) } for j : = 1 to 7 L M G ( / - 1, 7, uz, u2)\ PROLONGATE; {vt <- vt + I ^ v i ^ } for i : = l to u2 R E L A X ( £ ; u , = di); 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. The case 7 = 1 Chapter 3. The Multigrid Method 34 is known as the V(u1}u2) cycle while 7 = 2 is known as the W(ui,u2) cycle where vx and v2 are the number of pre and post C G C relaxations. A n alternative approach is to begin not on the finest level as in L M G but on the coarsest level. This method is commonly known as full mult igr id or F M G . Having obtained an approximation to the solution on the coarsest level this solution is inter-polated, 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, vx, u2); Begin S O L V E ( £ 1 u 1 = / a ) ; for I := 2 to p do begin I N T E R P O L A T E ; {ut <— f j ^ u ^ j } L M G ( / , 7, uu u2); end; End ; Final ly , it should be noted that for nonlinear operators, (3.1.6) may have to be used in place of (3.1.7) thus requiring that the solution-as. .well as the defect be transfered to each coarser level. This 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 wi l l not be considered further. Chapter 3. The Multigrid Method 35 3.3 M u l t i g r i d for coupled equations We now consider the application of the mult igrid method outlined above to the following problem. £V.«2) = / 1 (3- 3 - 1 ) £2(u\u2) = 0 (3.3.2) u1 = u2 = 0 on dQ (3.3.3) Recall that both the elastic and, i n part, the viscoelastic problem may be wri t ten in 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 d2h = 0 so that d\-.= dl-C\(ulu2h) (3.3.4) d2h:=d2h-C2h(ulu2h) (3.3.5) and a coupled system of equations for the corrections on n^ is also obtained as follows: 4 ( « ) = 4 (3-3.6) C l ( v l v 2 h ) = dl (3.3.7) where vl = v2h = 0 on dnh (3.3.8) The transfer of defects and corrections between coarse and fine grids remains standard since in 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 room for manoeuvre in this area. Basically we have to choose from the following: • Sequential Smoothing Chapter 3. The Multigrid Method 36 • Collective Smoothing We now consider each of these in more detail. Sequential smoothing is performed by successively smoothing one component of the solution while keeping the other fixed and then reversing the roles. So, in consideration of (3.3.6) and (3.3.7) for example, at each grid point we might first update vx using (3.3.6) while keeping v2 fixed and then update v2 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) in place of (3.3.6). We should expect these methods to work well when the variables v1 and v2 are not strongly coupled. A t any grid point equations (3.3.6) and (3.3.7) may be written: A V i \ C j \c2 ) (3.3.9) \ v2 where the 2x2 matr ix A and the right hand side ( c 1 , c 2 ) T are particular to the chosen discretization. If the matr ix A is strongly diagonally dominant then we may deduce that the variables v\ and v2 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. Computat ional ly this is more expensive than sequential smoothing but in cases where the variables are strongly coupled in one or other of the equations, this method is a necessary alternative to sequential smoothing. In our case, the added expense in 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. The Multigrid Method 37 3.4 Multigrid for time-dependent problems Consider the time-dependent problem Cu(t) = f(t) on ft (3.4.1) wi th the time-independent boundary condition: u = 0 on dn (3.4.2) Assuming a uniform discretization in time wi th time-step At, the discrete problem at the (n + l ) s t time step may be written: ^W + 1 = A n + 1 on nh (3.4.3) where < + 1 = 0 on dtth (3.4.4) A t each time step we may consider using L M G to solve the problem for u^. A t the first t ime step we take u°h = 0 on the entire domain as an in i t ia l guess to the solution. Clear ly at each successive time step it would be pointless to use zero as an ini t ia l guess since un 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 t ime step while st i l l maintaining reasonable accuracy. A n alternative to L M G is the F M G algorithm. Gendler [4] (following Brandt) and Hackbusch [5, p.272] suggest modified F M G F A S algorithms for dealing wi th time de-pendent problems. Hackbusch [5, p.273] also presents criteria under which the proposed algori thm wi l l converge. In our case, the F A S component is not required since we are dealing with a linear problem, however, the procedure outlined below follows the same philosophy. Chapter 3. The Multigrid Method 38 A g a i n the question uppermost in our minds is how information from the previous time step can best be used as an in i t ia l guess for the current time step. Consider the following: 1. Inject the solution un down to every level. 2. Solve exactly on the coarsest level to obtain u n + 1 3. Perform a high order interpolation of the correction, u n + 1 — un, to this solution and add it to the solution on the next finest level. 4. Perform one L M G cycle to obtain un+1. 5. Repeat 3. and 4. unti l 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. This prompts us to define the following algorithm: A l g o r i t h m 3.3 ( M F M G ) Procedure M F M G ( / , 7 , U2); Begin for / :=p - 1 to 1 do u n I N J « + 1 ) ; R E L A X ( £ 1 U 1 = / i ) ; u f + 1 ] for / := 2 to p do begin u n+l ;= u? + INT(u n + l 7-1 Chapter 3. The Multigrid Method 39 L M G ( / , 7, uu u2); {-> end; End ; where INJ represents the straight injection operator and INT is a high order interpo-lation. This algorithm involves more work than the L M G algori thm but only on the coarser grids so the cost is not significantly higher. It is obviously hoped that the additional expense results in a more accurate solution - this is yet to be established. In the M F M G algorithm we are making the assumption that the 0(2h)2 discretiza-t ion error arising from our approximation to u n + 1 on the next to finest level is a better ini t ia l guess to u n + 1 than the 0(h)2 approximation to un obtained at the previous time step. It is not clear why this should be the case. In the l imit as At —> 0 we can be certain that L M G should provide the more accurate solution since in this case the O ( A i ) error becomes insignificant in comparison to the 0(2h)2 error. A n estimate of the 0(2/i)2 discretization error would therefore be of use. The discretization errors, 8^-and Sfj1, at each grid point on the finest and next-to-finest levels respectively, are given by % := - ui: ^ ClJh2 (3.4.5) 6}} ':= u% - U i j ~ ClJ(2h)2 (3.4.6) Together these equations lead to the following: Chapter 3. The Multigrid Method 40 We wi 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 approximation is easily obtained using the M F M G algori thm 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 in the solution from one time step to the next to the discretization error, 82h, on the next-to-finest level. In cases where the discretization error is larger than the change in the solution we expect L M G to provide a better solution than M F M G , in Chapter 5 we wi l l show this to be the case., C h a p t e r 4 The Elas t ic P r o b l e m 4.1 A n a l y t i c Solu t ion First let us recall the governing equation wi th its boundary conditions: D V w + jw = p (4.1.1) where w ~ V"tt> = 0 on dfl (4.1.2) 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 domain 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 in the re-direction and Xy in the y-direction, so that the topographic displacement, hr, at any point (x,y) is given by: where ho is the maximum topographic amplitude. Hence, if pc represents the density of the rock forming the surface topography we have: (4.1.3) We now assume that the plate deflection responds in such a way that: (4.1.5) 41 Chapter 4. The Elastic Problem 42 where w0 is a constant to be determined. In other words, it is assumed that the plate deflection is also periodic wi th a period equal to that of the loading function. Notice that in order to satisfy the boundary conditions (4.1.2) we must insist that 2L/XX and 2L/Xy be positive integers where L represents the length of the domain. Now, substi tuting (4.1.4) and (4.1.5) into (4.1.1) we obtain: 2 - i - i (4.1.6) which represents the amplitude of the deflection. Whi l e the application here is obviously to topographic loading such as mountain ranges and valleys, equation (4.1.4) may also be seen to represent a simple model for the formation of sedimentary basins. This requires the assumption that A x = Xy = 2L 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 wi th the assumption that we are attempting to model deflection due to topographic loading. For most practical applications we wi l l be considering the case of mountain ranges which extend the length of the domain in the y-direction and vary periodically i n the x-direction so that Xy = 2L and nXx = 2L hence, Xy = nXx where n is a positive integer. We may therefore write (4.1.6) in the form: w0 ho 2 T T \ 4 D \ ( 1 y / 7 K.) \ Pc9J V1 n2 / \pcg (4.1.7) where now the two terms in the denominator are nondimensional. Turcotte & Schubert [10, p.123] point out that the quantity (D/pcg)1/4 has dimensions of length and is 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 than Chapter 4. The Elastic Problem 43 the natural wavelength so that: ( 1 + (4.1.8) then we have: Peg h{ •0 (4.1.9) 7 where wl represents the amplitude of deflection resulting from an isostatic adjustment of the lithosphere to the loading. In other words, the lithosphere may be considered to have no rigidity and to be in hydrostatic equil ibrium. Alternatively, if the wavelength of the forcing is much less than the natural wavelength, we have: and the deflection of the lithosphere is negligible. In other words, the r igidity 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 r igidity 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. At the opposite end of the spectrum, in cases where the plate's rigidity is not significant, we have the equivalent of a small parameter mul t ip lying the highest derivative in the equation. This yields a singularly perturbed problem which, although a very simple one, could again give rise to difficulties in obtaining a numerical solution. A s we proceed, we shall therefore bear in mind 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, so that: < h0 Cha.pt er 4. The Elastic Problem 44 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 wi l l compute the degree to which the plate's rigidity supports the load as this wi l l provide us wi th some measure of the importance of the biharmonic operator in the solution of the equation. Let R = (u>, — w0)/w.i so that for iy0 — 0, R. ~ 1 and we have 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 w0 — 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 in the solution of the governing equation. Substi tuting 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 Problem 45 L(m) D (Nm) I 1 x 102 1 1 x 10 2 3 1 x 10 2 5 2.5 x 10 5 11.3 92.7 99.9 5.0 x 10 5 0.8 44.3 98.8 1.0 x IO 6 0.05 4.7 83.2 L = Length of Plate D = Flexura l Rigidi ty Table I: Degree of support (Elastic Model) topographic wavelength. Hence, for a 50% support of the loading by the plate's r igidity we have: = 2 * 1 - 1 (4.1.121 1 + rt2 j \ 7 We may also express R as a function of the plate length, L, as follows - l R = r VD/ V l 1 n ' (4.1.13) Table I gives appropriate values of R for various values of D and L where a value of n — 2 has been assumed. The range of values of D covers most of those found in the literature for geophysical applications. The range of values of L covers everything from sedimentary basins to continental plates. The 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 m , the density of the mantle, to be 3300 kg m ~ 3 and p c , the density of the surface rock to be 2800 kg m ~ 3 . Hence, 7 has the value 4900 k g 2 m ~ 2 . Chapter 4. The Elastic Problem 46 4.2 N u m e r i c a l S o l u t i o n Recall that we propose to solve the problem numerically by reducing (4.1.1) to two coupled equations so that at the grid point we have: D^Vij+'ywij =Pij (4.2.1) V2Wij = Vij (4.2.2) with boundary conditions to = w = 0 o n 80, (4.2.3) We discretize the Laplacian operator using the standard five-point formula: V ^ u ' t ; = —(Wij-x + wlJ+1 + wz_ltj + wl+li3 - 4wli0) (4.2.4) where h represents the grid-spacing which we have assumed to be uniform over the entire domain. We use full weighting for restriction and bilinear interpolation for prolongation. Relaxat ion is carried out using Gauss-Seidel wi th 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. Final ly, 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 in each of the equations. Having chosen an appropriate discretization of the Laplacian, we may now write our equations in the form (4.2.4) where the 2x2 matr ix, A, becomes: ' 4/h2 - i j U \ (4.2.5) 1 4 /A 2 / Notice that the primes have been brought back into use to stress that numerically we deal wi th nondimensional quantities, see section 2.5. Sequential smoothing should be Chapter 4. The Elastic Problem 47 L ( m ) 2.5 x 10 5 5.0 x 10 5 1.0 x IO 6 5.0 x IO 6 ijU 1.9 31.0 4.9 x 10 2 3.1 x IO 6 L = Length of Plate 7'/D' = Nondimensionalized j/D. Table II: j' / D' as a function of plate length effective when A is diagonally dominant, this requires that: h2 D' anc J ? > 1 The 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 2 = 16. For sequential smoothing to be effective 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 10 2 5 (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 km. Notice that smaller values of D would result in an even stronger coupling for a given L. It is there-fore clear that sequential smoothing should be avoided and hence we choose collective Gauss-Seidel for the relaxation procedure throughout. As was pointed out earlier, this does not lead to any significant additional expense since we are dealing only wi th a 2x2 system. There is also the added benefit of obtaining an exact solution of the 2x2 system at each grid point. Chapter 4. The Elastic Problem 48 4.3 Numerical Results We now examine the performance of the suggested mult igrid 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 Loading In examining the performance of a numerical procedure there are two things to con-sider. Fi rs 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 wi th a good combination of both - clearly one is not much use without the other. As 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 maximum displacement wi l l be used in this capacity so that: 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. The easiest (4.3.6) where wl] is the numerical solution, u>,7 is the exact solution and w0 is the amplitude Chapter 4. The Elastic Problem 49 method is to examine convergence at one particular point and take this to be represen-tative of the overall convergence. Hence, it is imperative that we choose a point where maximum deflection, rather than zero deflection, is expected - knowing the exact solu-tion, this wi l l be possible. We wi l l terminate the iterative procedure at the kth iterate when (k) (k-(k) where the grid point is a location where max imum deflection is expected. In all cases TOL is chosen to be 1 x I O - 4 so that the iterative procedure is terminated when the solution at a point of maximum deflection has converged to wi th in one ten thousandth of its computed value. In Tables III, IV and V , we present results for various sizes of plates, flexural rigidities and topographic wavelengths. The number of iterations, I T N , required for convergence, the r.m.s. error, E , and the predicted, Rp, and exact, Re, degrees of support, are given in each case. The mul t igr id 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. The L M G algorithm wi th 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 mult igr id iteration generally converges quickly to a solution which is a very good approximation to the exact solution. This is worthy of note since for a small degree of support, the nondimensionalized equations have a small parameter mul t ip ly ing the highest derivative i.e. the biharmonic operator. We therefore have a singularly per-turbed problem for which our mult igr id method provides a good solution. However, in consideration of the discretized equations we see that in the l imit as the degree of < TOL (4.3.7) Chapter 4. The Elastic Problem 50 n D ITN E(%) Rp(%) Re(%) 1 x 10 2 1 4 0.13 11.1 11.3 2 1 x 10 2 3 5 1.09 92.6 92.7 1 x 10 2 5 5 1.18 99.9 99.9 1 x 10 2 1 4 3.12 57.2 59.5 4 1 x 10 2 3 4 5.42 99.2 99.3 1 x 10 2 5 4 5.46 99.9 99.9 Table III: L = 2.5 x 10 5 m, Periodic Loading n D ITN E{%) RP(%) Re(%) 1 x 10 2 1 5 0.01 0.8 0.8 2 1 x 10 2 3 3 0.52 43.7 44.3 1 x 10 2 5 5 1.16 98.7 98.8 1 x 10 2 1 4 0.42 7.7 8.4 4 1 x 10 2 3 4 4.87 89.3 90.2 1 x 10 2 5 4 5.46 99.9 99.9 Table IV: L = 5.0 x 10 5 m, Periodic Loading n D ITN E(%) Rp(%) Re(%) 1 x 10 2 1 3 0.001 0.05 0.05 2 1 x 10 2 3 3 0.04 4.7 4.7 1 x 10 2 5 ' 5 0.98 82.9 83.2 1 x 10 2 1 3 0.03 0.52 0.57 4 1 x 10 2 3 4 1.87 34.3 36.5 1 x 10 2 5 4 . 5.36 98.1 98.3 Table V : L = 1.0 x 10 6 m, Periodic Loading n = 2L/XX D — F lexura l Rigidi ty ITN = Number of iterations for convergence E = Relative error Rp = Predicted degree of support Re = Exact degree of support Chapter 4. The Elastic Problem 51 support tends to zero, equation (4.2.2) plays no role in the computation and equation (4.2.1) reduces to: wtj = ^ (4.3.8) 7 which is the solution of the reduced continuous equation representing an isostatic ad-justment of the lithosphere to the applied load. The 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 with 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 in the solution. This is likely caused by a combination of two factors. Firs 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. Overal l , we can be pleased wi th 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. Also , the r.m.s. error in the solution is at worst 5% of the maximum 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 Problem 52 L D G S - i t n M G - i t n Re{%) 1.0 X 10 6 1 x 10 2 1 5 3 0.05 5.0 X 10 5 1 x 1 0 2 3 58 3 44.3 2.5 X 10 5 1 x 10 2 5 94 5 99.9 L — Length of Plate D = Flexural Rigidi ty G S - i t n = Number of Gauss-Seidel iterations M G - i t n = Number of Mul t i g r i d iterations Re — Exact degree of support Table V I : Mu l t i g r i d vs. Gauss-Seidel 4 .3.2 M u l t i g r i d v s . G a u s s - S e i d e l In order to appreciate the efficiency of the mult igrid 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 maintain TOL at 1 x 10~ 4 . The results are presented in Table V I where GS- i tn is the number of straight Gauss-Seidel iterations and M G - i t n is the number of mult igr id iterations required for convergence. The work involved in relaxing the 4 levels of one mult igr id 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 mult igrid. However, we suspect that the degree of support does not have to become too large before mult igrid performs"relatively more efficiently. Mu l t i g r i d certainly proves to be more efficient when the degree of support is at 44% or higher. We can therefore claim with some confidence that in cases where the flexural r igidity is significant, the mult igr id method described proves to be a very efficient solver for the problem. Chapter 4. The Elastic Problem 53 4.3.3 Gaussian Loading We now turn our attention to a loading function for which an analytic solution may not be obtained. Consider a Gaussian distribution for the forcing so that The application here is to the formation of sedimentary basins where a basin of depth hr is formed by thermal contraction. We therefore attempt to predict the subsequent deflection, w, due to the sedimentary infill of this basin. The 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. The parameter a in some sense represents the wavelength of the loading. Notice that maximum deflection wi 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. We again expect that smaller wavelengths in the loading result in greater support from the plate's rigidity. We therefore consider two choices for a, namely 0.04 and 0.08. Bo th 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. The results are presented in Tables V I I , VI 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. This is as expected since a larger value of a results in a larger wavelength in the forcing which in turn results in a lower degree of support as was seen for periodic loading. Aga in we can be satisfied wi th the convergence rate but of course we no longer have any means of assessing the accuracy of the solution, however, the fact (4.3.9) where (x - 0.5) 2 + (y- 0.5) 2 Chapter 4. The Elastic Problem 54 a ITN D Rp{%) 4 1 x 10 2 1 28.3 0.04 5 1 x 1 0 2 3 85.5 5 1 x 10 2 5 99.8 4 1 x 10 2 1 13.3 0.08 5 1 x I O 2 5 76.9 5 1 x 10 2 5 99.7 Table V I I : L = 2.5 x 10 5 m, Gaussian Loading a ITN D Rv(%) 5 1 x 10 2 1 4.6 0.04 5 1 x 10 2 3 52.7 5 1 x 10 2 5 96.8 5 1 x 10 2 1 1.4 0.08 5 1 x 10 2 5 32.4 5 1 x 10 2 5 94.9 T I L L = 5.0sx 10 5 m , G aussian . a ITN D Rp(%) 3 1 x 10 2 1 0.3 0.04 4 1 x 1 0 2 3 17.2 4 1 x 10 2 5 74.8 3 1 x 10 2 1 0.1 0.08 5 1 x 10 2 5 6.9 4 1 x 10 2 5 60.6 Table I X : L = 1.0 x IO 6 m, Gaussian Loading a = Forcing function parameter, see (4.3.9) ITN — Number of iterations for convergence D = F lexural Rigidi ty Rp = 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. Chap te r 5 The Viscoelas t ic P r o b l e m 5.1 A n a l y t i c So lu t ion We begin by restating the governing equation wi th its ini t ia l and boundary conditions: 74 • , (w , _-.\ P DV*w + 7 ^ - +wj = ^ + p (5.1.1) where, assuming the exact solution is known, w = We(t = 0) on (1 (5.1.2) and w = V 2 u ; = 0 on dn (5.1.3) As for the elastic problem, we seek an analytic solution to which we may compare future numerical results. We 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 with time. As an approximation to this situation, bearing in mind that we need an analytic solution, we assume a loading function of the form: p = Pcgh0e^-^ s i n s m ^ ( 5 . L 4 ) for which we have a corresponding deflection: w = sin ( ~ ) sin (5.1.5) 56 Chapter 5. The Viscoelastic Problem 57 where . o -I -1 1 1 w0 = (m + l)pcghQ 1 6 Z W 4 ( — + — ) + ( m + l ) 7 and t e is the time at which the integration is terminated. Notice that at t = te, the (5.1.6) forcing reaches its maximum value of . (2irx \ . (2-Ky\ Pmai = Pcgho sin J s m \ j (5.1.7) which is the same as that considered for the elastic problem. As for the elastic problem, if we assume Xx = A^ = 2L we may regard (5.1.4) as representing the load due to the infilling of a basin of depth V m ( ' " l e ) / T - The basin, assumed to be formed by some process such as thermal contraction, therfore increases in depth exponentially over a period of te mi l l ion years. W i t h an appropriate choice of Xx and Xy the boundary conditions (5.1.3) may be 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 in larger rates of increase in 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. Taking Xy = riXx = 2L as in the elastic case we may again derive an expression for the degree, R} to which the load is supported by the plate's rigidity. For the case outlined above we have: R 1 + - ^ 1 + - (5.1.8) .TTJ \DJ V mJ V l + n 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 w i l l increase for a given L and D. As can be seen, the case L = 1.0 x 10 6 m is not particularly interesting as the lithosphere is close to being in isostatic adjustment for all choices of Chapter 5. The Viscoelastic Problem 58 L ( m ) D (Nm) 1 1 x 10 2 1 1 x 1 0 2 3 1 x 10 2 5 2.5 x 10 5 0.12 11.2 92.6 5.0 x 10 5 0.0078 0.78 44.1 1.0 x 10 6 0.0005 0.05 4.7 L — Length of Plate D = Flexural Rigidity-Table X : Degree of support (Viscoelastic Model) 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. w = u V2u = v 5.2 N u m e r i c a l S o l u t i o n Numerical ly the problem is solved by first rewrit ing (5.1.1) in the form (2.7.1) as follows: (5.2.1) (5.2.2) (5.2.3) These equations may now be viewed as a differential algebraic equation ( D . A . E . ) in time. The time component can be regarded as differential in w and algebraic in 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 Backward Euler (B .E . ) which is a first order method. Knowing quantities at time tn, the problem at time tn+1 may be written: (5.2.4) ,n+l _ w + U At n+l T (5.2.5) Chapter 5. The Viscoelastic Problem 59 X72un+1 = v n + 1 (5.2.6) E l imina t ing wn+1 from (5.2.5) using (5.2.4) we have: + 7 ( ^ + 1) un+l = gn+1 (5.2.7) V2un+1 = vn+l (5.2.8) where gn+1 = - +p-1wn (5.2.9) r r is known at each time step. We now have two coupled elliptic equations for u and v. Recal l that in section 2.6 we showed that our boundary conditions may be writ ten equivalently as: u = 0, v = 0 on dn (5.2.10) Thus, at each time step, our problem may be writ ten in the form (2.7.6) wi th boundary conditions (2.7.7) where ui = un+1 and u2 = vn+1. Hence, we may use a mul t igr id algori thm to solve for u n + 1 at each time step and then use (5.2.4) to solve for wn+1. A s an alternative, following Ascher :[1, p.12], we may consider computing w at time steps tn while u and v w i l l be computed at 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 in time which requires no more work or storage than Backward Euler. So we have wn+1 -wn . un+? = (5.2.11) At v ' ^-1 + u"+§ = 1 +pn+1i (5.2.12) VV +2=i; n +§ (5.2.13) Chapter 5. The Viscoelastic Problem 60 Now, wr i t ing w n+= = -(wn+i+wn) 2 = l-(Mun+* +2wn) (5.2.14) at each time step i n + 2 , we use mult igr id to solve the system: + 7 + i ) = (5.2.15) V 2 u n + 2 = u n + § (5.2.16) where (?"+' = ^ - ^ - Iiy" (5.2.17) T T is known at each time step. Having solved for un+i we use (5.2.11) to solve for wn+1. Second order accuracy may also be obtained (with a little extra work) using the two-step backward difference formula ( B . D . F . ) : 1 / 3 ^ , _ „ 1 u n + 1 - — -wn+1 - 2wn + - u ; 7 1 " 1 ) (5.2.18) A t V2 2 J K .' This method also has greater storage requirements than Backward Euler or the C . D . F . approach and is slightly more awkward to implement if a change in step size, A t , is contemplated. Using (5.2.18) we have the following system of equations: DV2vn+1 + 7 + l ) wn+1 = .-gn+1 (5.2.19) V2wn+1 = vn+1 (5.2.20) where now ~n+l „ n + l A i gn+l = t + pn + l _ _wn + _wn-l (5.2.2T T 3 3 Again , this system may be writ ten in the form (2.7.6) and our mult igr id algorithm used once more. Chapter 5. The Viscoelastic Problem 61 We use the same 5-point operator to discretize the Laplacian and again use full weighting for the restriction and interpolation operators. Collective Gauss-Seidel is called upon 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 in the present formulation of the problem we must make some restrictions on their variabili ty in order to-make meaningful numerical •tests. We wi l l allow n to take on only the value 2 since this results in 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 in form to Gaussian loading, as well as larger values of n which may result in some error being introduced to the solution because of poor representation of the forcing function on the finest grid. Notice that in restricting ourselves to the case of n = 2 we are in some sense losing geophysical applicabil i ty since the application to sedimentary basins requires the assumption that n = 1. The case n = 2 is applicable to topographic loading although it is hardly reasonable to assume that the topography grows exponentially wi th time when in fact surface topography generally decreases in size after formation due to erosion. A more meaningful forcing function, from the geophysical standpoint, wi l l be considered Chapter 5. The Viscoelastic Problem 62 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 in the forcing. Increasing m yields a smaller forcing at t = 0 which is desirable, at least from a geophysical standpoint, but this also increases the rate of growth of the forcing which provides us wi th a tougher problem to solve numerically. Clearly, we wish to find a happy balance. Values of m — 0.01 and 0.03 result in an in i t ia l forcing having approximately 37% and 5% of its final value, respectively. The second of these is more reasonable geophysically but we expect it to be more difficult to solve numerically. Following Nunn & Sleep [7, p.607], we take the relaxation parameter, r , to have the value 1 x 10 6 years. Nunn & Sleep [7, p.588] also indicate that any physical theory should allow continual subsidence for a period of at least lOOMa in the modell ing of sedimentary basins. Integration wi l l therefore be carried out from 0 to lOOMa, however, all equations have been nondimensionalized wi th respect to a time scale T = lOOMa so that numerically, the equations are integrated from 0 to 1. We begin wi th 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, Esoi and Eder, calculated as in (4.3.6), for the solution w, and its derivative u, for 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, S2h, calculated as in (3.4.9), and an r.m.s. representation of the exact difference in the derivative u, at time t = lOOMa are also presented. The L M G and M F M G algorithms were tested using Backward Euler, the C . D . F . method and the two-step B . D . F . method all wi th 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 mind that we are dealing wi th both a time and Chapter 5. The Viscoelastic Problem 63 Esoi(%) Eder(%) ph A4 D L M G M F M G L M G M F M G ( x l O - 2 ) ( x l O - 2 ) 1 x 10 2 1 0.081 0.029 2.506 2.553 0.0079 0.3406 1 x 1 0 2 3 0.259 0.470 2.323 2.127 0.1794 0.3029 1 x 10 2 5 1.420 2.839 0.088 1.669 0.3444 0.0251 Table X I : Backward Euler , L M G vs. M F M G Esol(%) Ed ph A* D L M G M F M G L M G M F M G ( x l O " 2 ) ( x l O " 2 ) 1 x 10 2 1 0.050 0.066 0.063 0.037 0.0120 0.3406 1 x 10 2 3 0.090 0.123 0.071 0.146 0.2404 0.3029 1 x 10 2 5 0.144 1.119 0.180 1.808 0.0330 0.0251 Table X I I : C . D . F method, L M G vs. M F M G Esol(%) Ed ph A* D L M G M F M G L M G M F M G ( x l O - 2 ) ( x l O " 2 ) 1 x 10 2 1 0.034 0.002 0.139 0.163 0.0113 0.3406 1 x 10 2 3 0.413 0.197 1.624 0.035 0.2268 0.3029 1 x 10 2 5 0.248 1.153 0.176 1.808 0.0347 0.0251 Table X I I I : Two-step B . D . F . , L M G vs. M F M G D = Flexura l Rig id i ty Esol = Relative Error in Solution Eder = Relative Error in Derivative ph = Discretization Error = Exact change in derivative Chapter 5. The Viscoelastic Problem 64 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 2 5 N m . In these cases, L M G consistently provides a better approximation to both the solution and the derivative than M F M G . It can also be noticed that in these cases the dis-cretization error on the next-to-finest level is significantly larger than the change in the derivative and so we expect L M G to perform better than M F M G as discussed in section 3.4. For other values of D, things are not quite so clear cut. For D = 1 x 1 0 2 1 N m where the discretization error is consistently smaller than the change in 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 in only one out of the three cases. Overal l , L M G provided a better approximation than M F M G to the solution in 6 out of the 9 cases and a better approximation to the derivative also in 6 out of 9 cases. We wi 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 . and two-step B . D . F . methods wi th 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. The reason for this is two-fold. Fi rs 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 Problem 65 m D Esoi(%) Eder{%) Rv{%) Re(%) 1 x 10 2 1 0.081 2.506 0.01 0.13 0.01 1 x 1 0 2 3 0.259 2.323 10.8 11.2 1 x 10 2 5 1.420 0.088 92.4 92.6 1 x 10 2 1 0.600 6.741 -0.68 0.37 0.03 1 x 1 0 2 3 1.350 6.092 25.2 27.0 1 x 10 2 5 6.284 1.615 97.1 97.3 Table-'-XIV Backward Euler , A t = l O M a m D Esoi{%) Ed„{%) Rp(%) Re(%) 1 x 10 2 1 0.050 0.063 0.25 0.13 0.01 1 x 1 0 2 3 0.090 0.071 11.3 11.2 1 x 10 2 5 0.144 0.180 92.6 92.6 1 x 10 2 1 0.458 0.239 1.36 0.37 0.03 1 x 1 0 2 3 1.193 0.997 28.6 27.0 1 x 10 2 5 1.490 1.423 97.4 97.3 Table X V C . D . F . method, A t = 10M a m D Esol(%) EdeT{%) &,(%) Re(%) 1 x 10 2 1 0.034 0.139 0.10 0.13 0.01 1 x 1 0 2 3 0.023 0.183 11.2 11.2 1 x 10 2 5 0.248 0.176 92.6 92.6 1 x 10 2 1 0.264 1.047 0.06 0.37 0.0? •1 x 1 0 2 3 0.413 1.684 27.6 27.0 1 x 10 2 5 0.122 1.447 97.3 97.3 Table X V I : Two-step B . D . F . , A t = l O M a m = Forcing function parameter, see (5.1.4) D — F lexural Rigidi ty Esoi = Relative Error in Solution Eder = Relative Error in Derivative Rp = Predicted degree of support Re = Exact 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. We expect the latter reason is the more significant since if we do not start wi th a good approximation to the solution at the current t ime 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 in both the solution and its derivative. Later we shall examine the effect, of reducing the time step. C . D . F . consistently provides better results than B . E . for the same amount of work making C . D . F . a very attractive method. C . D . F . also provides less spurious results in that as D increases, the errors in both the solution and the derivative are seen to increase as expected. This 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. They should both provide approximations which are second order in 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 in 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 addit ional work and storage involved in 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 wi 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 in less Chapter 5. The Viscoelastic Problem 67 m D Esol(%) Eder{%) RP(%) Re{%) 1 x 10 2 1 0.198 3.629 0.09 0.37 0.03 1 x 1 0 2 3 0.899 2.973 25.8 27.0 1 x 10 2 5 3.563 0.353 97.2 97.3 Table X V I I : Backward Euler , A t = 5 M a m D Esoi(%) Ed„{%) RP(%) Re(%) 1 x 10 2 1 0.120 0.074 0.63 0.37 0.03 1 x 1 0 2 3 0.298 0.249 27.4 27.0 1 x 10 2 5 0.278 2.257 97.4 97.3 Table X V I I I : C . D . F metb.od,..At = 5 M a m D Esol(%) Eder{%) £P(%) Re{%) 1 x 10 2 1 0.067 0.308 0.33 0.37 0.03 1 x 1 0 2 3 0.100 0.453 27.1 27.0 1 x 10 2 5 0.126 0.263 97.4 97.3 Table X I X : Two-step B . D . F . , A t =-5 M a m = Forcing function parameter, see (5.1.4) D = Flexura l Rigidi ty Esoi — Relative Error in Solution Eder — Relative Error in Derivative Rp = Predicted degree of support Re = Exact degree of support accurate solutions. As expected, reducing the time-step has led to more accurate results for al l three methods. Backward Euler, being a.'hrst order method, provides an O ( A t ) approxima-tion (in time) to both the solution and 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 in the solution and its derivative by 50%. The errors we have presented, however, account for errors arising from both the time and spatial discretizations and so we should not necessarily expect a 50% improvement in Chapter 5. The Viscoelastic Problem 68 our results. We expect that reducing the time step by a factor of 2 for the C . D . F . method should reduce the error in the solution by a factor of 4 whereas the error in the derivative should be reduced by only 50% since C . D . F . provides a second order approximation to the solution but only a first order approximation to its derivative. B . D . F . , on the other hand, provides a second order approximation to both the solution and its derivative and so we should expect the errors in both 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 in accuracy since we also have an error arising from the spatial discretization. Backward Euler in fact provides close to a 50% improvement i n the approximations to both the solution and its derivative. The C . D . F . method in some sense behaves better than expected in that a reduction of close to 75% is observed not only for the error in the solution, but also for the error in the derivative when only a 50% improvement in the latter is expected. B . D . F . provides a near 75% improvement for both the solution and its derivative, as expected, except for the case D = 1 x 1 0 2 5 N m where the error in the solution actually gets worse! A n explanation for this strange observation is extremely hard to find. On comparing the two second order methods, we see that B . D . F . once more pro-vides better estimates of the solution whereas C . D . F . provides better estimates of the derivative. In view of the added difficulty and additional expense in using B . D . F . we shall from hereon restrict ourselves to the C . D . F . method. We shall also discard Back-ward Euler since we clearly gain nothing in using this method from the point of view of either accuracy or efficiency. Chapter 5. The Viscoelastic Problem 69 5.3.2 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, in 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 > te we look for a solution of the equation: D V 4 u ) + 7 + w^j = P (5.3.1) where now, p is independent of time. The solution to this equation is: w=- (5.3.2) 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. Owing 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 te which we take to be lOOMa. For times greater than te, the forcing function has the form: . /2"7rx\ . / 27ry \ P = Pcgh0 sin [~^~J s m J (5.3.3) so that it maintains the value it had at lOOMa. We consider two cases having different values of L, D and m. The first takes L = 2.5 x 10 5 m, D = 1 x 1 0 2 3 N m and m = 0.03 while the second has L = 5.0 x 10 5 m, D = 1 x 1 0 2 5 N m and m = 0.01. We take At = l O M a for both and use the C . D . F . method. In Table X X we present computed values for the max imum deflection, wmax, at intervals of lOOMa for case (i). The maximum deflection expected when the lithosphere is in isostatic adjustment is 16800m. It is very clear that the numerical solution behaves as predicted and converges to the Chapter 5. The Viscoelastic Problem 70 t (Ma) 100 11995.21 200 16797.98 300 16798.74 400 16799.10 500 16799.76 Table X X : Steady state integration : case 1 t (Ma) 100 . 9408.81 200 14790.80 300 16253.75 400 16649.53 500 16756.77 600 16785.98 700 16794.06 800 16796.40 900 16797.17 1000 16797.49 Table X X I : Steady state integration : case 2 Chapter 5. The Viscoelastic Problem 71 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 unt i l it is in isostatic adjustment wi th 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 r igidity) wi l l deform wi th t ime unt i l it reaches isostatic adjustment. 5.3.3 Geophys ica l A p p l i c a t i o n Final ly , we turn our attention to a forcing function which is more meaningful geophys-ically but for which there is no analytic solution available. We wi l l assume a Gaussian distr ibution for the loading. So, allowing f(t) to represent the time dependence and taking hx to be as represented in equation (4.3.9) we have: p = P c g h T ( r ) f ( t ) (5.3.4) 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 ( — ) (5.3.5) 7 r \ r r i T J which has the required properties that / (0 ) = 0 and f(t) —> 1 as t —• oo. The application here is to the formation of sedimentary basins. A basin of shape KT is assumed to be formed by tectonic subsidence. The basin ini t ia l ly has zero depth but increases wi th time to reach a steady state asymptotically. The relaxation parameter, r , wi l l again take on the value 1 x 10 6 years so that the forcing grows by 50% in m M a . Chapter 5. The Viscoelastic Problem 72 m /(ioo) /(500) /(1000) 50 0.844 0.994 0.998 75 0.674 0.986 0.996 100 0.500 0.975 0.994 Table X X I I : T ime dependence in terms of m t(Ma) m = 50 m = 100 100 10102 5592 ,'.200 13783 11433 s* 300 14975 13752 •• 400 15537 14802 •• 5 0 0 : '15872 15381 600 16092 15743 700 16247 15987 800 16358 16158 900 16442 16284 1000 16505 16378 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 wi l l 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 in steady state towards the end of integration, (ii) the rate of increase at the beginning of integration is not too great. We wi l l consider the following two cases: (i) L = 2.5 x 10 5m, D = 1 x 10 2 3 Nm, (ii) L — 5.0 x 10 5m, D = 1 x 1025 N m and allow m to have the values 50 or 100 for both. Since, p = 0 at. t = 0 we have the in i t ia l condition w = 0 at t = 0. The C . D . F . method was used with one V ( l , l ) cycle per time step wi th At = l O M a . The results are presented in Tables X X I I I and X X I V . Chapter 5. The Viscoelastic Problem 73 t (Ma) m = 50 m = 100 100 6688 3565 200 9271 7585 300 10719 9579 400 11759 10904 500 12539 11873 600 13135 12607 700 13598 13173 800 13962 13617 900 14254 13969 1000 14491 14254 Table X X I V : M a x i m u m deflection (metres) for case (ii) The results show that for case (i) the deflection is close to being in isostatic adjust-ment towards the end of integration - this is as expected since the derivative, p, is close to zero at this stage and, as was seen in section 5.3.2, the steady state solution is one representing an adjustment to isostatic equil ibrium. We also see that increasing m re-duces the in i t i a l rate of increase in the loading but towards the end of integration there is l i t t le 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 in 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 wi l l be required before the plate reaches isostatic adjustment. This k ind of behaviour has already been observed in the previous section. The effect of increasing m is the same as was observed for case (i)-C h a p t e r 6 C o n c l u s i o n 6.1 R e s u l t s A derivation of the governing equation for deflection of a viscoelastic plate under an applied load has been presented. The 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. Boundary 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 prob-lems were presented. The discretization was applied after reformulating the fourth order equations in terms of two coupled second order equations. The mult igr id method has proved to be a very fast and efficient solver for the elastic, problem. It has dealt effectively wi th a wide range of plate lengths and flexural rigidities representing degrees of support ranging from 0.05% to 99.9%. The numerical solution has therefore proven to be stable wi th respect to varying degrees of importance of the biharmonic 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 in fact converges quickly to the solution of the reduced equation. As the degree of support increases and the biharmonic operator becomes more significant, the numerical solution takes longer to converge and becomes 74 Chapter 6. Conclusion 75 less accurate, however, this deterioration cannot be considered to be serious as the convergence rate and accuracy are st i l l wi th in acceptable l imits. A l l i n all we can be very pleased w i t h the performance of the mult igr id method as applied to the elastic problem. The 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 both the algorithmic level and on the level of choosing appropriate time discretizations. We examined two mult igrid 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)2 discretization error was obtained from the M F M G algorithm and compared to the exact change i n the derivative over the last t ime step. It was found that in those cases where the discretization error was higher than the change in 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 in the derivative. O n the whole, L M G provided better results at a lower cost and so the M F M G algorithm was discarded. The time dependence was dealt wi th i n three different ways: (i) Backward Euler - a first order method, (ii) A centred difference scheme which provides a second order approximation with the same amount of work as Backward Euler 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 Backward Euler scheme. A comparison of the two second order methods did 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. We also examined the effect of maintaining the forcing function at a steady state Chapter 6. Conclusion 76 after a given period of time and observing the subsequent deflection. Analy t ica l ly we expect the plate to relax, asymptotically, unt i l it is in isostatic adjustment regardless of plate size or elastic properties - this behaviour was also borne out in the numerical results. It appears that a very good degree of accuracy may be obtained wi th only one V ( l , 1) L M G cycle per time step. This represents a fast and highly efficient method. Perhaps the greatest advantage of the solution method we have developed lies in its adaptabili ty in that a min imum amount.of work is involved in solving the problem wi th different forcing functions once the program has been set up. This cannot be said of methods such as integral transform techniques and Fourier Series expansion methods which have been employed previously. The 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 in dealing wi th time-dependent mult igr id . 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. This appears to be due to the high discretization error on the next-to-finest level relative to the change in the solution from one time step to the next. However, further investigation into the applicabili ty 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 t ime 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 optimal choice. The 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 , in that integrating wi th a time step A t using, for example, one V ( 2 , 2) cycle per time step involves the same amount of work as integrating wi th a time step A t / 2 using one V ( l , 1) cycle per t ime step. This question of at taining opt imal efficiency at the algorithmic level is certainly worthy of further attention. In the geophysical context there is st i l l room for improvement in 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 pl is the density of the z t h layer having a thickness h,(x,y). A comparison of the elastic and viscoelastic models would also be of interest. The deflection of an elastic lithosphere wi th time-dependent loading is governed by the equation: DV4w(t) +jw(t) = p{t) (6.2.2) Having solved the viscoelastic problem, the solution of (6.2.2) should not present any further difficulty numerically since it represents a simple continuation problem. The L M G algorithm wi th 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 compar-ison 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 mod-elled, is the extension of the mult igr id method to deal with irregular domains. This involves much computational work since the discretization, interpolation and restric-t ion operators become non-standard at points near the boundary. The 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 unt i l only one interior grid point remains thus making 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 both the numerical and geophysical aspects of this topic are in order. We content ourselves wi th having provided what we hope to be a sound basis for future investigation in this field of study. Bibliography [l] Ascher, U . On numerical differential algebraic problems with application to semi-conductor device simulation S . I . A . M . Journal of Numerical Analysis (to appear). [2] Beaumont, C . The evolution of Sedimentary .Basins on a Viscoelastic Lithosphere Geophysical Journal of the Royal Astronomical Society 55 p.471, 1978. [3] Drucker, D . C . Introduction to the mechanics of deformable bodies M c G r a w - H i l l 1967. [4] Gendler, E . Multigrid methods for time-dependent parabolic equations Master 's Thesis, Weizmann Institute of Science, Rehovot, Israel 1986. [5] Hackbusch, W . Multi-Grid methods and Applications Springer Verlag 1985. [6] Hunter, S.C. Mechanics of Continuous Media E l l i s Horwood 1983. [7] Nunn , J . A . , and N . H . Sleep Thermal Contraction and Flexure of Iniracratonal ' Basins: a three dimensional study of the Michigan Basin Geophysical Journal of the Royal Astronomical Society 76 p.587, 1984. [8] Nunn , J . A . , N . H . Sleep and W . E . Moore Thermal Subsidence and generation of hydrocarbons in the Michigan Basin Amer ican Association of Petroleum Geologists Bul le t in 68 No. 3 p.296, 1984. [9] Turcot te, D . L . Flexure Advances in Geophysics, 21 1979. [10] Turcotte, D . L . and G . Schubert Geodynamics J . Wi ley 1982. 79 Bibliography 80 [11] Vinson , J .R . Structural Mechanics J .Wi ley 1974. [12] Watts , Karner and Steckler Lithospheric Flexure and the evolution of Sedimentary Basins Philosophical Transcripts of the Royal Society of London A 305 p.249, 1982. [13] Wesseling, P. in Multigrid Methods S.F. M c C o r m i c k ed., Frontiers in Appl i ed Mathematics , S . I . A . M . 1987.
- 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0079633/manifest