MSU LIBRARIES “ RETURNING MATERIALS: RTEce in book drop to remove this checkout from your record. FINES W111 be charged if book is returned after the date stamped be10w. INTEGRAL EQUATION METHODS FOR EIGENVALUE ESTIMATION IN SYSTEMS WITH DISCONTINUOUS COEFFICIENTS By‘ John Patrick Spence A DISSERTATION Submitted to Michigan State University in partial fulfiiiment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanica] Engineering 1981 ABSTRACT INTEGRAL EQUATION METHODS FOR EIGENVALUE ESTIMATION IN SYSTEMS NITH DISCONTINUOUS COEFFICIENTS By John Patrick Spence Systems with abruptly varying physical properties occur in many contexts, both natural and artificial. Of particular interest is the solution of eigenvalue problems which arise in the mathematical models of composite systems. The analysis of such problems hinges on the development of effective and appropriate approximate methods, since exact solution of eigenvalue equations is only possible in the simplest of discontinuous systems. It has been found that traditional techniques which are quite effective in determining the dynamic characteristics of systems with Slowly varying properties do not perform nearly as well in discontinuous systems without the application of considerable extra effort. In the present work, the advantages of taking an integral equation approach to problems in eigenvalue estimation in discontinuous systems are investigated. The inherent smoothing properties of integrals sug- gest that methods based on integral equations should be quite effective in handling systems with discontinuous material properties. It is found that the integral equation formulation of the Galerkin method leads to upper bounds for the eigenvalues which are superior to John Patrick Spence those obtained from the traditional differential equation formulation. Furthermore, the usually unavailable lower bounds to the eigenvalues may be readily computed by the use of trace identities involving the kernel of the integral equation and infinite sums of powers of the eigenvalues. These lower bounds complement any set of upper bounds and, by their convergence behavior, give insight into the distribution of the generally irregular spectra of discontinuous systems. ACKNOWLEDGEMENTS I wish to express my appreciation to the many people who have inspired, encouraged and supported me during the course of my research. First, I wish to thank my major advisor, Dr. Albert N. Andry, Jr., for convincing me to persevere in the doctoral program and for being an excellent sounding board for my ideas. I also wish to thank Dr. Cornelius 0. Horgan for bringing the problem addressed in this dissertation to my attention and for his continuing guidance and encouragement in bringing this work to fruition. I also appreciate the contribution of the other members of my guidance committee, DrS. David Yen and James Beck; for their careful reading and suggestions in clarifying the mathematical statement and areas of application of the methods used in this work. I am also grateful to all my instructors for the knowledge and intuition imparted to me during the course of my degree program and to Ms. Judy Duncan for her timely and accurate typing of this disserta- tion. But most of all, I wish to thank my wife, Jan, for her loving support and patience, without which this work would never have come to pass. ii TABLE OF CONTENTS List of Tables ................................................... v List of Figures .................................................. vi INTRODUCTION ..................................................... 1 CHAPTER l -- FUNDAMENTALS OF INTEGRAL EQUATIONS .................. 5 1.1 Basic Concepts of Linear Integral Equation Theory ................................ 5 1.2 Formulation of Integral Equation Eigenvalue Problems ............................ 10 CHAPTER 2 -- INTEGRAL EQUATION EIGENVALUE ESTIMATION PROCEDURES .......................................... 12 2.1 The Galerkin Method ............................ 12 2.2 Lower Bounds by Trace Identities ............... 15 CHAPTER 3 -- STURM-LIOUVILLE PROBLEMS WITH DISCONTINUOUS COEFFICIENTS ........................................ 18 3.1 Introduction and Literature Review ............. 18 3.2 Eigenvalue Problem Formulation ................. 20 3.3 Integral Equation Formulation .................. 22 3.4 Lower Bounds by Eigenvalue Extremization ....... 23 3.5 Example Problem ................................ 25 3.6 Numerical Results .............................. 29 3.7 Discussion of Results .......................... 30 CHAPTER 4 -— LONER BOUNDS FOR THE EIGENFREQUENCIES OF STEPPED ELASTIC BEAMS ............................... 39 4.1 Introduction and Literature Review ............. 39 4.2 Formulation of the Eigenvalue Problem .......... 41 4.3 Green's Functions and Integral Equation Formulation .................................... 41 4.4 Sample Problem #l--Cantilever Case ............. 43 4.5 Numerical Results and Discussion ............... 46 4.6 Sample Problem #2--Simp1y Supported Case ....... 49 4.7 Discussion of Results-~Prob1em #2 .............. 52 TABLE OF CONTENTS (continued) CHAPTER 5 -- CHAPTER 6 -- CHAPTER 7 -- HARMONIC WAVES IN LAYERED COMPOSITES ............... Introduction and Literature Review ............ Formulation of the Eigenvalue Problem ......... Integral Equation Formulation and Iterated Kernel Trace ......................... Sample Problem ................................ Numerical Results ............................. Discussion of Results ......................... 010101 010101 03014:: (JON—4 FREE VIBRATIONS OF A CIRCULAR MEMBRANE WITH RADIALLY STEPPED DENSITY .......................... Introduction and Literature Review ............ Formulation of the Eigenvalue Problem ......... Integral Equation Formulation ................. Sample Problem ................................ . Numerical Results ............................. 6. Discussion of Results ......................... momma mth—a SUMMARY AND SUGGESTIONS FOR FURTHER STUDY .......... 7.1 Summary ....................................... 7.2 Suggestions for Further Study ................. LIST OF REFERENCES .............................................. iv 55 55 56 60 62 64 65 90 94 LIST OF TABLES Comparison of Truncation and Extremization Lower Bounds for the Least Eigenvalue ...................... 35 Comparison of Six—Term Upper Bounds from Integral and Differential Formulations with Lower Bound Sets ..................................... 36 Bounds for the Spectrum with y = 10, e = 10 with Percent Relative Errors ......................... 37 Bounds for the Spectrum with y = 0.1, 6 = 0.1, with Percent Relative Errors ......................... 38 Bounds for the Least Two Frequencies of a Stepped Simply Supported Beam ........................ 48 Bounds for the Least Two Squared Frequencies of a Stepped Simply Supported Beam ................... 54 Lower Bounds to the Least Eigenvalue, Q = l .......... 68 Lower Bounds to the Least Eigenvalue, Q = 3 .......... 69 Lower Bounds to the Second Eigenvalue, Q = l ......... 70 Lower Bounds to the Second Eigenvalue, Q = 3 ......... 71 Bounds to the Spectrum with Q = 1, y = 10, e = 0.01 ............................................. 72 e = 0.01 ............................................. 73 Bounds to the Least Frequency of the Two- Piece Circular Membrane .............................. 86 Bounds for the Spectrum with a = l, e = .01 ......... 87 Bounds for the Spectrum with a = 7, e = 10,000 ...... 88 Bounds for the Spectrum with a = 9, e = 100 ......... 89 3.1 3.2 4.1 4.2 5.1 6.1 LIST OF FIGURES Layered Slab with Piecewise Constant Properties ....... 26 Fundamental Mode Temperature Profiles ................. 25 Stepped Cantilever Beam--Sample Problem #1 ........... 47 Stepped Simply Supported Beam--Samp1e Problem #2 ...... 47 Layered Composite with Periodic Structure ............. 57 Circular Membrane with Radially Stepped Density ....... 85 vi INTRODUCTION Systems with abrupt changes in physical properties occur in many contexts in nature and engineering. Among the naturally occurring systems are the layered structure of skin and teeth, the lipid bilayer biological membranes which separate regions of differing solute con- centrations, and rock layers in the earth. Problems of engineering interest include vibration and buckling of stepped shafts, wave pro- pagation through composite media, cooling of nuclear fuel rods and redistribution of impurities in semiconductors. Exact solution of the eigenvalue problems posed in the mathematical models of discontinuous systems such as these are only possible in the simplest of cases. Thus the analysis of such systems depends on the development of apprOpriate approximate methods. It has been found that traditional techniques which are quite effective in determining the eigenvalues of systems with slowly varying properties do not perform nearly as well in these discontinuous systems, unless con- siderable additional effort is applied. Further, accurate lower bounds to the eigenvalues are generally unavailable from these methods. In this dissertation, a series of eigenvalue problems arising from the application of the technique of separation of variables to linear partial differential equations with linear time invariant boundary conditions are analyzed. These problems are all essentially one dimensional, but the methods used are not restricted to cases which are or can be made to be such. As shall be seen in the following chapters, the crux of the integral equation formulations prerequisite to applying the methods of this dissertation lies in the availability or derivability of a static impulse response or Green's function. In one-dimensional cases over a finite region, the Green's functions may be found in a closed form which makes the eigenvalue problems in their integral equation form tractable by formally exact methods. However, in problems which cannot be reduced to one spatial dimension by techniques such as conformal mapping or separation of variables, the Green's function itself must be approximated. This, although by no means a trivial task, may be accomplished by methods such as finite differences or finite elements to any desired degree of accuracy. Once a suitably accurate impulse response has been determined, the methods of this dissertation may be applied to find eigenvalue estimates for the new model defined by the Green's function approximation. Past work in the area of eigenvalue estimation includes use of Laplace transform methods [1,2] finite difference methods [3-5], Sturm-Liouville theory [6-9,12] and many variational schemes. Among the variational approaches taken are Rayleigh-Ritz schemes using both smooth, non-smooth and piecewise defined test functions [6,7,10-17,19] and saddle point variational principles [12-16,19]. Lower bounds for eigenvalues have been obtained by Weinstein's method of intermediate problems [18] and application of results from Sturm-Liouville theory [6,7,12] and eigenvalue extremization [7,19,20]. Some of the lower bound results based on integral equation techniques from the third chapter of this dissertation also appear in [20]. More detailed expositions of the advantages and pitfalls of these methods are in the application chapters of this dissertation. In the present work, the advantages of taking an integral equation approach to problems in eigenvalue estimation are investigated. The inherent smoothing properties of integrals suggest that methods based on integral equations should be quite effective in handling systems with discontinuous material properties. For example integral equation methods have proven quite effective in a problem of beam vibration with random coefficients [21]. Variational methods, such as the Rayleigh- Ritz or Galerkin methods, developed in the context of differential formulations of the problems, may also be applied to the integral equation formulation to obtain upper bounds for the eigenvalues. Add- itionally, lower bounds to complement the upper bounds from variational methods may be readily computed, a distinct advantage over the dif- ferential formulations. The lower bounds in each case are presented in terms of basic parameters of the system in question and thus allow a parametric analysis of the dependence of the eigenvalues on system characteristics. In Chapters 1 and 2 of this dissertation, the fundamental concepts of integral equation theory and some of the results applicable to eigen- value estimation are discussed. Chapter 3 presents the results of applying both upper and lower bound integral equation methods to heat conduction in a layered slab, comparing the results to those obtained by eigenvalue extremization and differential equation variational methods. In Chapter 4, lower bounds are computed for the frequencies of vibration of stepped beams, and are compared to those from the method of intermediate problems. Chapter 5 addresses a problem of wave propagation, where lower bounds are found for the frequencies of non-dispersive wave forms travelling normal to the interfaces of laminated composites. In Chapter 6, upper and lower bounds are found by integral equation methods for the frequencies of vibration of a circular membrane with a stepped radial density. The results of this study, together with some suggestions for further study, are summarized in Chapter 7. The results of this study indicate that integral equation methods may be profitably employed in the analysis of the dynamics of composite structures. CHAPTER 1 FUNDAMENTALS OF INTEGRAL EQUATIONS This chapter presents the basic definitions and theorems of the theory of linear integral equations used in the subsequent chapters of this dissertation. The particulars of the integral equation based techniques used to bound eigenvalues are presented in Chapter 2. For a more detailed exposition of the theory of integral equations, see the books of Cochran [22], Tricomi [23] and Stakgold [24]. 1.1 Basic Concepts of Linear Integral Equation Theory Generally speaking, any equation which incorporates an unknown function under one or more signs of integration is an integral equation. Foregoing technical assumptions concerning smoothness, two particular classes of interest are the Fredholm integral equations of the first kind Mxl = fa k(x,y)¢(y)dy (l .1) and of the second kind MM=OU)+A§MMHMHW (La In these equations, o is the unknown function, A is a complex parameter, and w and k are known functions. The function k(x,y) is called the kernel of the integral equation. Notationally, the independent variables are suppressed and the integral is written as KO 2 gfk(x,y)¢ s 4Ff(x)g(x)dx where the superior bar denotes complex conjugation. Definition 1.3 Two L2(a,b) functions f and g are said to be orthogonal if = 0. Theorem 1.1 The complex inner product <-,-> satisfies the following properties for f, g, heL2(a,b) and a, B complex scalars: (1) = a + B (2) <9,f> = 3_O with equality if and only if f = 0 (almost everywhere). Theorem 1.2 Let f and g be elements of L2(a,b). Then the functional is an innerproductnorm; that is (l) llell = IAI ||f|| where A is a complex scalar (2) l|f+9|| :,||f|| + llgll (Minkowski inequality) (3) ||f|| 3_0 with equality if and only if f = 0 (almost everywhere). Further, the inner products <-, -> and norm - satsify (4) l|_: ||f . g|| (Schwartz inequality) Theorem 1.3 L2(a,b), under the inner product norm . , contains the limits of all Cauchy sequences of elements of L2(a,b). Thus L2(a,b) is a Hilbert space (complete normed linear vector Space with inner product norm). Definition 1.4 The adjoint k*(x,y) of the kernel k(x,y) is given by the complex conjugate of the transpose of k; that is k*(st) = ENYSX): and, in operator notation, the adjoint K* of K is defined by * - b 'k K ch = fa k (x,y)¢(y)dy. Remark: An alternative definition based on the L2 inner product is that K* is the operator with kernel k* satisfying, for all f. g e L2(a.b). = < f,K*g> Definition 1.5 The trace of the operator K with L2 kernel k(x,y) is given by _ b tr(K) : fak(x,x)dx Definition 1.6 The composition k of two kernels k1 and k2 is given by b k(x,y) = fal k1(X.2)k2(Z.y)dz and one writes in operator notation K = K1K2. Theorem 1.4 Let K1 and K2 be linear integral operators with L2 kernels k1 and k2. Then the composite kernel of K E K1K2 is L2, and (1) 11K11=11K1K211:11K] K211 9 Further, if o is L2(a,b) so is w = K¢ with (2) llwll i IIK ¢>l| Definition 1.7 A kernel k(x,y) (or operator K) is said to be Hermitian if k = k* (or K = K*). If k is real valued and Hermitian, it is said to be symmetric. If KK* = K*K, K is said to be normal; all Hermitian kernels have this property. Definition 1.8 A Fredholm integral equation of the second kind, 4 = w + AK¢, is said to be homogeneous if u = 0. Definition 1.9 A value of the complex parameter A in the homogeneous equation ¢ = AK¢ which admits only the trivial solution 6 E 0 is called a regular value of K. Values of A admitting a nontrivial solution o i O are called eigenvalues of K, and o is an eigenfunction of K belonging to A. Theorem 1.5 Let K be a linear integral operator with L2 kernel k(x,y). Then: (1) The eigenvalues of a Hermitian kernel form a non-void finite or countable sequence {An} of real numbers ordered by their magnitude 0 < |A1| < | 2| < ...; this sequence has no finite limit point. " " (2) The eigenfunctions of the Hermitian kernel k belonging to distinct eigenvalues are orthogonal and may be chosen to have norm 1; that is, they can be chosen orthonormal. (3) The eigenfunctions of a symmetric kernel k may be chosen to be real. (4) The finite number of eigenfunctions belonging to any distinct eigenvalue A of a Hermitian kernel k may be chosen to be orthonormal. (5) 1.2 10 If 61, oz, ... ¢N are distinct orthonormal eigenfunctions of a Hermitian kernel k belonging to the eigenvalues A], A2, ... A (not necessarily distinct), then Bessel's inequality ho ds: N |¢1(x)|2 b 2 .21 ———;?r—*’ j: .g lK(X,Y)I dy 1: Integration of the inequality (5) leads to m4: IIKII2 1 M2 —1 1: and if each eigenvalue of K, repeated according to its multiplicity, is included in the sum, equality holds. Formulation of Integral Equation Eigenvalue Problems One way of obtaining linear integral equation eigenvalue problems is to transform linear differential eigenvalue problems to integral form. This is accomplished by using an associated Green's function which captures both the internal and boundary behavior of the underlying differential equation in terms of the response to unit point excitation [25]. For example, examine the following linear differential eigen— value problem: L[u] B[u] AM[u] on domain 0 (1.4) 0 on boundary 30. The operators L, M, and B are linear and do not depend on the eigen- parameter A. If one then can solve the associated problem L[G(X.y)l = 6(x-y) on D B[G(x,y)] = O on 30, 11 where L and B are understood to operate on the X dependence of G(X,y) and 6 is the Dirac delta distribution, the solution to an equation L[u] = v, B[u] = O is given by u(x) = IDG(x.y)V(y)dy. In particular then, the eigenvalue problem 1.4 can be expressed as U(X) = AIDG(x,y)MLUJ(y)dy (1.5) In the case where L and 8 lead to a self-adjoint problem, that is, if for all functions u and v (for which the equation makes sense), = , the resulting Green's function is Hermitian. However the integral equa- tion 1.5 does not have a Hermitian kernel unless M[u](y) = m-u(y) for some scalar m f 0. In the case where M is algebraic, M[u](y) = m(y)u(y), with m(y) > 0, the equation 1.5 can be made to have a Hermitian kernel by the substitutions ¢(-)=fifil_-Tu(-) k(x,y) = memwmw (1'6) giving the eigenvalue problem with Hermitian kernel k, ¢(X) = Aka(x,y)¢(y)dy or (1.7) ¢ = AK¢. CHAPTER 2 INTEGRAL EQUATION EIGENVALUE ESTIMATION PROCEDURES This chapter outlines the basis and mechanics of the procedures used in this dissertation to provide upper and lower bounds for the eigenvalues in the examples of the subsequent chapters. Upper bounds to the eigenvalues are obtained by the Galerkin method, which for Hermitian kernels coincides with the traditional Rayleigh-Ritz method. Lower bounds are computed using these upper bounds and the trace of a related kernel. For a more detailed exposition of these methods and many others, see the excellent book of Baker [26]. 2.1 The Galerkin Method The essence of the Galerkin method is the use of finite expansions of the true eigenfunctions in a series of appropriately chosen known functions. That is, a true eigenfunction w of the eigenvalue problem 0 - AKw = O is approximated by {13 Ill 355% (2.1) where the aj are as yet undetermined constants and the ‘bj are linearly independent functions. If the approximant is substituted into the expression w - AKw one obtains an error term r(x), w - AKv = r(x). (2.2) 12 13 Various minimizations of the error term lead to the expansion approxima- tion methods. In the general Galerkin method, the constants aj and A are found by requiring that the error term be orthogonal to the test functions used under a weighted inner product with known positive L2 weight w(x) defined by w E £Ef(x)§(x)w(x)dx. One possible choice of w is w(x) E 1, leading to orthogonality under the usual L2 inner product. If Equation 2.2 is multiplied by 6} and integrated over the domain, one obtains §8j<¢j’¢1> ' A §3j = for each i = l, 2, ... n. If the matrices A and B are defined by Aij <¢j’¢1> (2.3) requiring that the error term be orthogonal to the test functions used leads to the general matrix eigenvalue problem for the a1 and A (A - AB)§ = 0. (2.4) The choice of linearly independent test functions guarantees that A is nonsingular. Further, it is clear that A is Hermitian, since by Theorem 1.1, <¢j’¢1> = <¢1’¢j>' 14 If the kernel k(x,y) is Hermitian, the matrix B is as well, since = <¢j’K*¢i> (Definition 1.4) = <¢j’K¢i> (Definition 1.7) = ' (Theorem 1.1) That the approximate eigenvalues obtained by the Galerkin pro- cedures are upper bounds in magnitude to the actual eigenvalues is a consequence of the Poincare characterization of the eigenvalues, namely max min 1 = An 8(n) ¢ES(n) /<¢a¢> (n) where S ranges over all n-dimensional subspaces of L2. This result is stated in the following theorem proven in Baker [26], page 321. Theorem 2.1 Let k(x,y) be a Hermitian L2 kernel, o], ... on linearly independent L2 functions, and A and B the Galerkin matrices from Equations 2.3a and 2.3b . Then if A:(A;) is the rth positive (negative) eigenvalue satisfying (A - AB)3 = 0, then . - + - - . prov1ded Ar, Ar(Ar,Ar) both eXTst. A choice of a complete sequence of test functions leads to the monotone convergence of the approximate eigenvalues to the actual and thus completeness in L2 is a desirable property for the test functions chosen to have. Further characteristics of the test functions leading to generally better results are that they share as many of the known properties of the true eigenfunctions as possible, such as boundary 15 conditions, number of zeros in the interval, and, of some consequence in discontinuous systems, smoothness conditions at internal points. It should be noted that the convergence of the approximate eigen- functions as more terms are added is in a mean-square sense to the actual eigenfunctions. In cases where there are jumps in the actual eigenfunctions (or their derivatives), the use of continuous (or smooth) test functions can lead to poor uniform approximations. These effects are ameliorated by the use of K¢j rather than the oj in the eigenfunction approximations; the iterated Galerkin method [27]. 2.2 Lower Bounds by Trace Identities Define the iterated kernels k(n)(x,y) by the recursion k"’ A.'n, j = l, 2, ... i=l1—J 16 giving the lower bound Aj > (tr(k(")))"/". (2.6) This bound is of course most accurate for the least eigenvalue. If a set of M upper bounds Xi is available, then the inequality (2.6) may be improved considerably: tr(K(n)) = :30).- | v >2 + M >J I v >J + M >2 giving the lower bounds N Am > (tr(K(n)) - z ‘A'i'nfl/n —' ifm (2.7) This bound differs from each m, and is improvable by obtaining either more or better upper bounds. In this dissertation, the bounds obtained from the particular case of n=2 are used and are denoted by _I§°) 2 (tr(K‘2)))"/2 (2.8) and M _;é“1 ; (tr(K(2)) z 'xi'2)"/2, (2.9) 17 and are called the truncation lower bound and the corrected truncation lower bound respectively. CHAPTER 3 STURM-LIOUVILLE PROBLEMS WITH DISCONTINUOUS COEFFICIENTS 3.1 Introduction and Literature Review This chapter addresses the problem of eigenvalue estimation for Sturm-Liouville problems with discontinuous coefficients in the context of diffusion in a laminated medium. Other problems of the Sturm- Liouville form that have been the subject of much recent attention include vibration problems in geophysics [8,9,28], buckling of stepped beams [12,18], and harmonic waves in layered composites [6,10-15,l9], among others. The latter problem is discussed in Chapter 5 of this dissertation. Considerable emphasis has been placed on the development of compu- tational schemes for estimating the eigenvalues and eigenfunctions for such problems. These efforts have met with serious difficulties due to the non-smoothness of the coefficients and the resulting spectral irregularities. Early attempts were focused mainly on variational tech- niques, with emphasis on obtaining upper bounds for the eigenvalues. These techniques include Rayleigh-Ritz approximation using smooth test functions and improved test functions with appropriate derivative dis- continuities [7.12.19]. Also in the variational mode, mixed variational principles where two field quantities are independently varied have been employed [7,12,19]. Alternative methods, such as finite difference 18 19 and finite element methods, which lead to matrix eigenvalue problems, have been investigated [3-3. A direct variational scheme employing piecewise polynomial eigenfunction approximation has also been used [17]. More recently, results from classical Sturm-Liouville theory and eigenvalue optimization techniques have been adapted for these problems to obtain upper and lower bounds for the eigenvalues [7,12,19]. The most valuable contribution of [7] was the reduction of the general Sturm-Liouville problem to that of a string with discontinuous density. In this simpler form, the various techniques are more easily applied and give better results for a comparable amount of effort. Lower bounds have also been found by a variation on Weinstein's method of intermediate problems [18]. A complete theoretical discussion of dif- fusion in laminated media under general linear interface conditions has been set forth in [29], with an extensive bibliography. In this chapter, integral equation methods are applied to a problem of heat conduction in a layered slab with fixed temperature boundary conditions and perfect thermal contact interface conditions. Lower bounds to the least eigenvalue obtained using the iterated kernel trace are compared to those based on eigenvalue extremization discussed in [19]. Upper bounds are generated using Galerkin's method applied to both the differential and integral equation formulations of the problem. These upper bounds are compared in terms of their accuracy and are also used to generate the sequences of lower bounds for higher eigenvalues obtained by correction of the truncation of the series summing to the iterated kernel trace. The actual eigenvalues,obtained by numerical solution of the actual transcendental eigenvalue equation, are used 20 to compute errors in the bounds. Numerical results are presented for a variety of material property combinations, and sources of error in the approximate methods are discussed. The results suggest that integral equation techniques, with their inherent smoothing properties, may be particularly appropriate for investigation of the problems of concern here. 3.2 Eigenvalue Problem Formulation The general partial differential equation modelling diffusion with no internal sources is v-(KvU) = c g3 (3.1) In the one-dimensional problem of heat conduction in a layered slab with fixed surface temperature and perfect thermal contact between the layers, separation of spatial and temporal variables leads to the Sturm-Liouville problem: (K(x)u'(x))' + A C(x)u(x) = 0 O < x < L u(O) = u(L) = O (3.2) u(x), K(x)u'(x) continuous for 0 < x < L The coefficient functions K and c, representing conductivity and capacity respectively, are piecewise continuous positive functions with step discontinuities at the interface locations x1, x2, ... xn‘ It has been recently demonstrated [7] that conversion of the problem (3.2) to Liouville normal form leads to computational advan- tages. Thus, let 21 K-](S)dS, t = T-1 A? K-](S)dS and (3.3) v(t) = u(x(t)). f(t) = T2K(x(t))c(x(t)). Then the eigenvalues of the system (3.2) are the same as those of V + Afv = 0 0.: t_: l v(O) = v(l) = 0 (3.4) v, v continuous for O < t < l where the superposed dot represents differentiation with respect to the new independent variable t. The coefficient function f(t) is positive and bounded for 0.: t.fi l and admits step discontinuities at the points ti = T.1 Cfi K-](S)dS. The effect of the transformation has been to remove the discontinuous coefficient K from its position subject to differentiation in the system (3.2) and move the locations of the discontinuities. When standard variational methods using smooth test functions are applied to the transformed system (3.4), the re— sults are equivalent to those obtained using improved test functions satisfying the condition of continuity in flux, Ku', with less computa- tional effort. A further advantage of this form is the availability of lower bounds to the eigenvaluesbased on eigenvalue extremization in vibrating strings. 22 3.3 Integral Equation Formulation The eigenvalue problem (3.4) is readily transformed to the integral equation v(t) = A 10' G(t,s)f(s)v(s)ds (3.8) where f (l-s)t 0 §_t < s < l G(t.S) = 4 (3.9) (l-t)s O §_s < t < l L is the Green's function, obtained by solving 2 - LG- (t,s) = 6(t-s); G(o.s) = G(I.s) = 0. (3-10) Btz Note that if we view the kernel of the integral equation as E(t,s) = G(t,s)f(s) (3.11) then the kernel is not symmetric in s and t. However, we may obtain an integral equation with symmetric kernel by multiplying (3.8) by /T(t) and thus obtaining the equivalent eigenvalue problem w(t) = A 10' k(t,s)w(s)ds = my (3.12) where w(t) = JTTET v(t) and k(t,s) = JTTET G(t,s) JFTST. Equation (3.12) is a Fredholm integral equation of the second kind with a kernel which, although discontinuous in s and t, is real, symmetric, square-integrable and continuous in the mean [30]. The trace of the second iterate of this kernel, which equals the L2 norm of the kernel, is given by the double integral 23 tr(K(2)) = [0' 10' f(t)GZ(t,s)f(s)dsdt (3.13) In the case of piecewise constant f(t), this equation involves easy integrations of polynomials of degree 3 3 in s and t over appropriate subregions of the unit square, and are easily carried out either by hand or numerically. 3.4 Lower Bounds by Eigenvalue Extremization Lower bounds for the eigenvalues of (3.2) have been obtained in [19] based on the normal form (3.4) and the results of Krein [31] on eigenvalue extremization. In [31], Krein'hsconcerned with the problem of maximizing and minimizing the eigenfrequencies w§(f) of a string of variable density f, as in (3.4), subject to the constraints on f of a fixed total mass 10 f(t)dt = M (3-14) and various bounds on the density, such as O < f §_H with H 3.M (3.15) 0 < h §_f with h < M Under conditions (3.14), (3.15b) the minimum A1 is attained by the singular function f(t) = h + (M-h)6(t-l/2) (3.16) which may be viewed as a uniform string with a bead at the center. Under the conditions (3.14), (3.15a) the minimum A1 is attained by 24 H |2t - 1| < M/H f(t) (3.17) 0 otherwise. The maximizing densities for the higher eigenvalues An have n equally spaced beads or dense intervals at the antinodal points of the nth eigenfunction of a uniform string. Thus we see that Krein's work leads to consideration of problems directly related to the investiga- tion of eigenvalue problems with discontinuous coefficients. This relationship was exploited in [19] and their development follows. Under the conditions (3.14), (3.15a), Krein [31] has shown that A (f) > 513% X(%) (3.18) An(f) ZTXHT) (3.19) where x(d) is the least positive root of the transcendental equation Ji— tan JET = d/(l-d) (3.20) The lower bounds may be made explicit on obtaining bounds for the lowest root. One such bound, obtained by Krein [31], is -1/2 x(d) > W - 5‘39 + (£- + iii-NZ] (3.21) This result, in conjunction with (3.18), was utilized in [19] for n = 1. When used with (3.19) more bounds are found; these results are pre- sented in the next two sections on the example problem. 25 3.5 Example Problem In order to evaluate the various bounding techniques, we consider the example of a three layer slab composed of two identical homogeneous outer layers enclosing an inner homogeneous layer as has been treated in [7,12,19]. The piecewise constant material property coefficienusK and c are given by K], c1 in l_: |2x-l| > b (outer layers) (3.22) K2, c2 in |2x-l| §_b (inner layer) The problem is then nondimensionalized and parametrized by the sub- stitutions K c y = E; , O = 52-, n1 = 1-b, n2 = b, 1 1 (3.23) K = n1K1 + n2K2, c = nlc1 + n2c2, and so, using the normalization :2= 1, E'= 1, we obtain c = (n + n G)'] c = G(n +n O)-1 1 l 2 ’ 2 l 2 _] 1 (3.24) K] (1'11 + nzY) 9 K2 = Y(n'| + nzY)- The corresponding dimensionless eigenvalue is then denoted by v and given by v = (IE/E)”2 = AUZ. (3.25) For given values of the geometric parameters 111 and n2, the effect of the material discontinuities on v is conveniently analyzed through consideration of the dependence of v on the dimensionless material K2 = 0.1‘K] C2 = 100::l XL. (KT')' ACT 0 T(O) = T(L) = 0 T,KT' ntinuo FIGURE 3.2 Fundamental mode temperature profiles 27 parameters y and 6. For continuous conductivities y = 1, while for continuous capacities 0 =1. Henceforth, we will set b = 1/2, giving 111 = n2 = 1/2 and interfaces at x1 = 1/4 and x2 = 3/4. For this case, the exact eigenvalue equations, obtained by solving (3.2) exactly in each material subinterval and matching temperature and flux at the interfaces, are: II 0 /A5 SinBO Sinav - cost cosav (3.26) JAB cost sinav + sinBv COSaV = 0 where 0': 1+0 ’ 8:4— .ppa 1+y 1/2 1 1+v" 1/2 (——) ( _,) . 1+9 The two separate equations for v1, v3, ... and v2, v4, ... result from the symmetry of the problem about the slab center. Under this parametrization, the coefficient f(t) in (3.4) is given by '11 It - 11> ' 1 2 217117' f(t) = l (3.27) 1 1 1 h2 It ' 21 f-ETTTTT where 3 h1= 1’” , 112 =Y0h]. 4v (1+9) For this f(t), the trace of the second kernel iterate is given by 28 2 tr(K(2)) ._. h${-L M Q(]-Q,)3(3+2,2) + 90 ' 144 2 £1%%5ll-2(15-1022+324)I (3.28) where C = (1W)-1 is the center interval length and ye is the coefficient discontinuity ratio in the Liouville normal form. The truncation lower bound (2.8) for the least eigenparameter v becomes 9].: (tr(k(2)))"/4, (3.29) and the lower bounds for higher eigenparameters (2.9) becomes v > mod”) - z (3.)'4)"/4 m_ f 1 (3.30) T1 111 An additional parameter to be used in the bounds from Krein [31] is the total "mass“ 2 M = (;f(t)dt = §1§%1—- (3.31) Thus from (3.18) and (3.19), using (3.20) and the above parametrization, we obtain the lower bounds 4ny1/2 4 11 1L. 2 -1/4 VnZ'TTFT—[1'3d+(24+fl4)d] (3.32) where d is given by one of the four expressions: d = 1+3"1 1+) ’ if v0.3 1 under f :.H = hz, 29 _ HY"1 d-W, 1fY931 under th=h], d=—'l9—, ifyO_h= 112. 1+6 3.6 Numerical Results In order to assess the accuracy and nature of the error involved in the lower bounds based on the iterated kernel trace, a number of bounds were computed for various combinations of material properties. The lower bounds for the least eigenvalue based on truncation of the series at the first term (3.29) are Shown in Table 3.1, with the best Krein bound (3.32) and actual eigenvalues from (3.26) for comparison. Upper bound sets from the Galerkin method using trigonometric test functions in the differential and integral formulations are presented for four cases, together with the corresponding lower bound sets from (3.30), in Table 3.2. An additional study of the progressively improv- ing upper and lower bounds obtained as the number of test functions used in the integral Galerkin method increases appears in Tables 3.3 and 3.4. Also in these tables, the best possible lower bounds from (3.30), where the actual eigenvalues are used as the upper bounds, are presented for comparison purposes. All of the numerical results in this chapter were computed in double precision on the Prime 750 computer of the Case Center for Computer Aided Design at Michigan State University. 30 3.7 Discussion of Results Table 3.1 provides a comparison between the integral equation lower bounds (3.29) and the eigenvalue extremization bounds (3.32). In all cases where the Liouville normal form (3.4) has a center weighted density f(t), that is, where the product of the discontinuity ratios, Y6, exceeds or equals one, the extremization bounds (3.32) are quite accurate. When both ratios Y and e exceed one, the only cases con- SldETEd'UlL7,12,l9],the integral equation lower bounds (3.29) compare favorably with the extremization bounds (3.32). In mixed cases, where either Y or a, but not both, and the product v6 is less than one, the coefficient function f(t) is no longer similar to the extremizing densities (3.16) and (3.17) of Krein, and the corresponding bounds (3.32) are quite poor. In these edge-weighted cases, the accuracy of the integral equation bounds (3.29) suffer as well, but for a quite different reason. The error induced by using (3.29) depends on the appropriateness of the truncation of the series that sums to the iterated kernel trace, and thus depends on the distribution of the spectrum. The closer the least eigenvalue is to the rest of the spectrum, the less accurate the results of the truncation can be. Table 3.1 shows that the accuracy of (3.29) is exceptional for highly center weighted cases with ya >> 1 and falls off for the edge weighted cases with ya < 1. Using the vibrations of a string as an analogy, one can interpret this phenomenon in terms of the effect of the heavier region on the modes of vibration. In the center weighted cases, the concentration of the mass in the center about the antinodal point of the first mode significantly lowers 31 the fundamental frequency compared to that obtained for a uniform string. However, this mass concentration occurs about the nodal point of the second mode, and the lesser amplitude of motion reduces the effect on the frequency. Thus the first two frequencies are more widely separated, and the series can be truncated at the first term and still give good results. In contrast, in the edge weighted cases, with ya < l, the motion of the concentrated masses near the ends is comparable in the first two modes, and the frequencies are correspondingly closer than those of a uniform string. Such cluster- ing of the lower eigenvalues increases the error due to truncation in (3.29). This defect is remedied by the use of upper bound information as in the bound (3.30), as can be seen in Table 3.2, case 4, where Y = e = 0.1. In Table 3.2, two sets of six upper bounds and the resulting lower bound sets are presented for four combinations of material properties. Both upper bound sets were obtained using GalerkinS method with test functions of the form 6 . t = .. ' ' v1( ) jilcustflt in the differential formulation, and 11m = «1117 v,(t) in the integral formulation. In this particular example problem, the symmetry of the system about the midpoint allows separation into odd and even matrix eigenvalue problems, since the test functions and true eigenfunctions are even and odd about the midpoint, and the 32 resulting inner products, which constitute the Galerkin matrices, vanish for pairs with different parity. As can be seen in all four cases in Table 3.2, the upper bounds from the integral equation formulation are superior to those from the differential equation formulation. This leads to a marked superiority of the lower bounds based on the integral formulation. The error in these lower bounds from Equation (3.30) comes from two sources: truncation of the series after six terms and the approxima- tion of the second, or correction, term. The error due to truncation depends on the number of terms used as well as the distribution of the spectrum. The error due to the approximation of the correction term depends primarily on the accuracy of the upper bound for the least eigenvalue used in the correction term. In Table 3.2, Case 1, the eigenvalues are nearly equally spaced and the upper bounds are quite accurate. In this case, the lower bound error is primarily due to series truncation, the accuracy of the lower bounds increasing only slightly when the better integral equation bounds are used. In Case 2, the spectrum is again nearly equally spaced, but the differential upper bound to v] is significantly less accurate than the integral, leading to far less accurate lower bounds for the rest of the spectrum. This difference is even more pronounced in Cases 3 and 4, but the effects of eigenvalue spacing enter into these bounds as well. Within each case, the accuracies of lower bounds within a cluster of closely spaced eigenvalues are comparable. In Case 3, the lower bounds to the least two isolated eigenvalues are accurate (particularly those from the integral upper bounds), and the lower bounds to the cluster v3 - us are of 27-36 33 percent relative error, jumping to 53 percent for the isolated 06. In Case 4, one cluster is v] - v3 and the corresponding lower bounds range from 2-8 percent relative error, v4 is isolated with 34 percent relative error, and 05 - O7 is another cluster with 51-56 percent relative error in the lower bounds. Tables 3.3 and 3.4 present another examination of Cases 3 and 4 designed to isolate the effects of the number of bounds used and of the spectral distribution from that of the accuracy of the upper bounds used. Thus in these tables, the best possible lower bounds obtainable from Equation (3.30), where the actual eigenvalues are used as upper bounds, are given along with the integral formulation upper bounds and their resulting lower bound sets. Again, within a cluster of closely spaced eigenvalues the error in these best possible lower bounds is seen to be comparable for a given number of upper bounds. One obvious effect of the number of upper bounds used is the increase in the accuracy of the lower bounds. A more subtle effect is the interaction with the effect of clustering. When upper bounds for all members of a cluster are used, there is a significant improvement in the accuracy of the lower bounds in that cluster which is greater than the improve- ment obtained upon adding a bound to an isolated eigenvalue. For example, examine the relative error changes down the columns of Table 3.3. When the third upper bound is added, the relative error in columns 1 and 2 drOps sharply; the addition of the fourth upper bound cuts the error by a more typical rate. A similar phenomenon occurs in Table 3.4, columns 4-6. 34 One other source of error in these series approximations which must be mentioned is that of the numerical truncation errors inherent in summing numbers of disparate magnitudes and in taking smaller differences of larger numbers. Both problems occur in the use of the corrected truncation lower bounds, especially in bounding the higher eigenvalues. Any use of this method must take account of these sources of error by summing from small to large and using high precision arith- metic. 35 TABLE 3.1 Comparison of Truncation and Extremization Lower Bounds for the Least Eigenvalue Y 6 0.1 1.0 10. 100. EV 5.058 2.421 1.8063 1.73177 0.1 LK 1.190(76.) 1.430(41.) 1.8014(.27) 1.72653(.30) LI 4.447(12.) 2.364(2.4) l.7709(2.0) 1.69822(1.9) EV 4.211 3.142 2.5292 2.44333 1.0 LK 2.488(41.) 3.133(.27) 2.4877(1.6) 2.43646(.28) LI 3.947(6.3) 3.080(2.0) 2.5043(.99) 2.42059(.93) EV 1.806 1.454 1.2251 1.19037 10. LK 1.801(.27) 1.430(1.6) 1.1901(2.9) 1.18666(.31) LI l.771(2.0) 1.448(.37) 1.2246(.04) 1.18994(.04) EV 0.596 0.484 0.4100 0.39868 100. LK 0.595(.30) 0.482(.28) 0.4087(.31) 0.39737(.33) LI 0.587(l.5) 0.483(.27) 0.4099(.01) 0.39867(<.001) EV = actual eigenvalue LK = best of Krein lower bounds, Equations 3.32 LI = iterated trace truncation lower bound, Equation 3.29 Figures in parentheses are percent relative errors. 36 TABLE 3.2 Comparison of Six-Term Upper Bounds from Integral and Differential Formulations with Lower Bound Sets Case 1: y = 0.1, e = 100 v 1.731772 3.469050 5.21637 6.9767 8.751 10.54 12.3 UD 1.731776 3.469163 5.21736 6.9804 8.767 10.58 LD 1.731284 3.453456 5.10183 6.5369 7.625 8.34 9.43 UI 1.731772 3.469052 5.21642 6.9770 8.754 10.55 LI 1.731301 3.454027 5.10578 6.5510 7.655 8.39 9.54 Case 2 v = 100, e = 0.1 v 0.59643 1.298409 1.79655 2.5966 3.0140 3.8943 4.25 U0 0.59697 1.298409 1.80962 2.5966 3.0676 3.8943 LD 0.59622 1.263957 1.65526 2.0044 2.1017 2.1796 2.24 UI 0.59644 1.298409 1.79844 2.5966 3.0364 3.8943 LI 0.59627 1.290145 1.75801 2.3813 2.6311 2.9364 3.24 Case 3 1'= 10.0, e = 10.0 v 1.22511 6.283 11.34 12.57 13.79 18.85 23.9 UD 1.23884 8.165 12.52 14.79 24.63 27.76 LD 1.22477 2.656 2.66 2.66 2.66 2.66 2.66 UI 1.22515 6.422 12.06 12.96 22.52 26.12 .LI 1.22503 6.038 8.29 8.41 8.78 8.80 8.83 Case 4 'Y= 0.1, e = 0.1 v 5.058 6.283 7.508 12.57 17.62 18.85 20.1 UD 5.316 8.165 10.827 14.79 21.84 27.76 LD 4.601 5.368 5.552 5.62 5.65 5.65 5.65 UI 5.088 6.422 7.879 12.96 19.87 26.12 LI 4.950 6.022 6.936 8.33 8.64 8.70 8.72 v = actual eigenvalues UD upper bounds, differential formulation UI. upper bounds, integral formulation corresponding lower bounds corresponding lower bounds HO 37 moapn>zompm pmauuu h: :o woman mucaoa Euro. u an «canon Luna: #2 :0 women mucaon Luzo— n ma covuapascow —ucumu=_ .vOgums :vxco_ao soc» mvcaoa Loan: u m: mucaon Luna: so» vow: mcovuucaw ummu Co guess: u z _ w:—m>:mmvm pmauus u > Am_-vm.o~ A¢.K-V¢.K_ Ae.~-ve.m_ AK.F-V3.~_ A~._-V_~.__ __.-vm~.c Amoco.-VKo.m-.. an mm-vm._. Ace-vn.._ A.~-.m.o. .c~-Vp.c_ .¢_-VmK.a K..-vfi_.e once.-vmsom-.p as o_ e.+v~.K~ A_m+vn.e~ Ae~+ve.5p Aea.+vs.~_ Ac.e+vom.P. mm.+vpn.o mooo.+vn._m-._ a: AmN-.~.K_ AK_-VK.m_ Am.m-va.~. .K.¢-vc.~_ .~.m-vKa.o_ Am”.-.o~.m Amocc.-vec_m-._ am po-ve.m cm- «.8 mm-vm.a Aa~-va.m AMN- 8K.» o.~-.~p.o A8888.-vomom-._ as K m~+vo.om am+ ..o~ oe+v..o~ A..m+.c.n_ A~.m+ mo... ~.~+.~¢.o Ao_oo.+v-_m-.. a: ow- N.K. om- _.m_ Am.m- e.~_ ...e- m.F_ Am.¢- om.op Ace- w~.o MKooc.- No.m-._ an me- m.m mm- m.m Ann- m.m Ann- e.w ANN- ow.» 3.”- 48.9 moco.- ~mow-._ as m Amm+v_.o~ Amo+vm.- A..m+vc.n. Ae.o+vec.~_ ~.~+v~v.o Amuoo.+vm¢_m-._ m: o~-v_.mp A~_- _.~_ Ae.m- e... AK.G-me.c_ AeK.-V.~.o ..oo.-voaom-._ an mm-.a.~ me- m.K Aam- K.K Amm-vom.K A~.e-v~c.o o__o.-vuoo¢-._ mg m mm+vm.- Am.m+vs.m_ Ae.m+voo.~. Am.m+vmo.c Am~oc.+vmepm-._ a: A~_-VF.~_ Am_-vm.o_ Am.-v~m.a .m._-.K_.o AKNco.-VKKom-.. am Ame-v_.K Ae¢-vc.K .mm- mm.o c.m- mK.m M_~_o.- Neoe-._ a; e Am.m+vfi.m_ Am.w+ .m.~_ m.o+ no.9 oc_c.+ ~m~m-._ a: m.- m.c_ cw- No.3 ..m- ao.o Macao.- Seem-._ am av- 3.8 me- m~.e K.m- 3K.m oewo.- aow¢-.p a; m Am.m+vpm.~_ Ao~+v¢m.K Aoo_8.+v~m~m-.. a: Low-vKo.o A..m- Ko.m ”moo.-vwoom-.. an Amm-v_m.. Aom- no.4 oKNO.-vaKKe-._ m4 N Ao~+vem.K A_mho.+.Kooo-._ a: A..m-.Km.m A4338.-Voome-._ am Aow-vmo.o «438.- meme-._ m4 _ _mhc.+ Kmoo-._ a: o - h: 3 mm m a, m.m_ e.~_ em.., m~.o mc_m-._ IF mLoLLm m>_um_wm ucmucoa cum: op u ¢ .o— u > saw: Eaguumnm on» Low mvcsom m.m m4m:mmpo paauuu b: no woman mvcaoa gaze. an mvcsoa song: #2 co woman mucsoa Lozop m4 cowuu—assow —ngamu:P .co;uos cwggm—mc soc» mucson song: a: mvcaon Lona: Low tum: mac—aucsm and» we guess: #2 * wapu>cmovm —o=uum *9 AK.a-v_.m. Am.K-ch.K_ Am.o-VFm.o_ .w.p-vqm.wp Aem.-vae.K AN..-VGKN.Q Ame.-vmmmo.m an ~e- K.__ mm- 38... “mm- Km... Am_- o~.o_ o.m- m~.K m._-v_a_.o NK.-VK_No.m mg a. mm+ ...m pm+ es.¢~ Am.m+ o~.m_ Aea.+ ae.~p No.+ mm.K ov.+ve_m.m ...+vmmoc.m a: Am.-vc.o_ Ac.-.Kw.m_ .m_-vpm.m_ Ae.¢-v~o.~_ A_o.-vm¢.K Acm.-veo~.e Am..-vo_mo.m am Nm-vo.m Ame-vom.m Ac¢-Vme.m Am~-voa.m 3.8-Vmo.K Am.~-vmnp.e. Ae._-vmoKa.e a; K pm+v~.~n Aom+v~_.e~ Ap.m+vmo.ap A_.m+vom.~_ a._+vmm.K A~.~+v-¢.o Ae~.+v~cuc.m m: op- e.o_ N~-v~e.e_ Aa_- e~.e_ A~.K- om... ...-vm¢.K Aem.- ao~.o m~.- coco.m an em. K.m em-voK.m “.m- «8.x Aem- mm.m e.K-Vem.o A~.e- -o.o ._~- «cma.e m; o Amm+vw..e~ Am_+.5m.a_ Am_+vom.~_ Ao.e+vmm.K A~.~+v-¢.m Aam.+vKKmo.m a: A--v~o.¢_ .mN-VK~.m_ .o_-vK~.__ A~._-me.s Lem.-v_m~.o Ana.-v_o¢c.m am “am-vmx.s Mom- mK.K Lam- mo.K A~_- om.e e.¢- eoo.o m.m- ammo.e m3 m m_+ Km.a_ Am.m+ mo.m_ Aa.e+ mm.~ m.o+ wKo.o om.4 KKmo.m a: Amm-vo~.m. A¢_-Vem.o_ A..~-vmm.K AN..-vKo~.o “Nm.-.K_mo.m an floo-vma.o Ame-va.o m_- Ne.e K.m- omK.m m.e-.-.m.e m; e Ao.m+vmo.m_ K_+ mK.m m.e+ m~o.o o._+voum_.m m: e_-vem.o_ Ao.m-.m_.K o.~- mp_.m ~._-Vmaom.e am om-vm~.o A_~-vmm.m 4.3- Nae.m ¢.K-v~mmo.¢ as m AK_+VwK.m Ac~+vamm.K Ao.p+voKm_.m m: Ao.m-vm_.K A__-qum.m Am.m-vmomk.e am Am~-VKK.m m_- oem.m ~.a- Pmam.o m4 N o~+ omm.K ~.m+ m-~.m a: A__-vewm.m A~_-VMK¢e.e am Amp-voom.m A~_-VMK¢¢.Q m4 _ A~.m+vm-~.m a: #2 ..ON mw.w_ No.K_ Km.~_ _m.K mmm.m pmmc.m , 9 mcoccu m>wum—mm acmugwm new: ~.o u o .p.c u > new: Eacaumam as» Low muczom v.m m4m m a an» area .om.¢ cowuascm =_ > can J» .z m1 . ANV» .m.¢ :oPuumm 5., an 9 .> cw muvmpu acuxo uwsammm msu co vmmmn men mocamPC omw;p« .uoma mew: mucaoa cmaa: > :o comma mucaon Lozop us» ogm zvl .m.~ covenacm scum ucsoa cwzop on» m, ona .mopu soc» :oxmu mvcaon Lona: ~upm-gm_mp»am Egmuuwmsgu new .mcorumspxocaam acupuoza so: .mmzpm>=mm_m —u:auu one «so 4» use .z> .> "p.¢ m—amh on ncumug Km_mK.K NemepK.c Kmme.K -- mmo~.mm me~m.m mamm.op «mesmo.o mmesmo.o -- amehmo.c Ko~_.P msmc.o mee.o cc, oo— mommh.e 33¢8~._ c~5o8.m -- Npmo.m NKm_.m mmam.m «mammo.o mmmmme.o -- mqmmmo.o F¢em.o oemo.o ammo.o m cop 4~¢~Km.K 4O¢m88.¢ «_ofimm.o -- mmmP.m mmmm..m Kch.K m_ammm.o Kmemm.o spammm.o mpmmmm.o Amvamwm.o Amvmmmm.o Am_vammm.o . co. mpro.om ~¢~om.m weem.m~ -- mmom.Pm omme.am eoom.om mmmmm.m mmmmm.m -- Kemmm.m meo¢.m Kmmm.m mmmm.m oo_ m oKom.mN em_m.~N ¢m_m.K~ -- mm__.m~ ¢m__.m~ meNm.m~ msm_e.¢ Mkmpe.e msm_e.e mpm_¢.e mm_¢.¢ amp¢.¢ empe.e cop _ mm qm 2N .m m9 29 9 r E E E 2: .. 8 Emom Lm>mppgcmu umaamum a $0 mmvocoacmgm o3» “mam; any Low muczom p.¢ m4m) 2 < 0] (4.29) M l M 2 -— -2 --— 9; ) = (tr(K( )) - Z ¢n ) 2 < ¢m (4.30) nfim where the $6 are any upper bounds. The value of tr(K(2)) in terms of system parameters and for x1 = 0.5 is, in matrix product form, F791 710 1757—1 7 (2) _ 08 2 tr(K ) - 2 [1 m0 m0] 710 1372 710 bO , (4.31) 6144-9450-bO 2 L175 710 791‘Lb0d Bickford's results,and the lower bounds for the least two eigenvalues based on the lower bounds of equations 4.29 and 4.30 using his upper bounds,are summarized in Table 4.2. 4.7 Discussion of Results--Problem #2 As can be seen in Table 4.2, the lower bounds based on the iterated kernel trace either with or without using upper bound information give comparable results in bounding the least eigenvalue to the far more expensive bounds from the method of truncation. In fact, the bounds computed from the trace and only the four low precision bounds of the twenty-four Rayleigh bounds are superior in all cases considered. These lower bounds would be even better if all digits and all of the upper bounds were available. In the second case, the bound for the second eigenvalue was based on a pessimistic extra digit in $$’24, 53 since the reported figure from [18] is not an upper bound. In all likelihood, the true upper bound was much better, and would have yielded a much better lower bound, though even the reported bound is quite good. Also of note is the contrast between the results of the second and fourth cases, where the eigenvalues are the same. The coincidence of the spectra, which Bickford found remarkable, is due to the fact that the cases are duals of each other. That the results are the same when b6] replaces m0 and ma] replaces b0 is apparent from an examina- tion of the equation for tr(K(2)). As in all variational methods, discontinuities in coefficients associated with differentiated terms lead to poorer results. Thus when a dual problem is available, one should analyze the one which has a lesser discontinuity in the coef- ficient of the most highly differentiated term. In the fourth case, with discontinuous stiffness, Bickford used the beam functions r} to improve the results obtained from using the 1 p.. 1 The most accurate figures from each of these are those reported in Table 4.2. These bounds are much more expensive than those which were obtained for the dual second case, and far less accurate as well. The question remains, however, whether a comparable computational investment in a Rayleigh-Ritz method, either through the use of more terms, improved test functions, or integral equation formulation, rather than in the method of truncation, would pay off by giving far better upper bounds and comparable lower bounds based on the iterated kernel trace. 54 TABLE 4.2 Bounds for the Least Two Squared Frequencies of a Stepped Simply Supported Beam b0 m0 $32,16 $64,32 $R,24 9(0) 9(4) 4 1 1.5198 1.5452 1.5945 1.57118 1.57298 32.907 33.237 33.789 -- 9.24096 [[1 $8,4 ¢l6,8 $—R,24 9(0) 9(4) L, 1 0.5 1.3262 1.3262 1.3262 1.32394 1.32622 23.111 23.111 23.111 -- 22.9544* $32,16,4 $64,32,4 $-R,24 i(O) 2(4) 2 0.5 1.7701 1.8212 1.8479 1.83768 1.84064 33.132 33.199 33.215 -- 17.6167 4, 2 64, 2 —- , 4 $3 3 ¢r 3 ¢ R 6 9(0) 9(1) 2 1.0 1.3129 1.3149 1.3281 1.32394 1.32394 23.077 23.082 -- -- 16.7328 Legend to Table 4.2 ¢n,k,£ lower bounds from method of truncation E’R’N N term Rayleigh-Ritz upper bounds 9(0) Lower bound based on Equation 4.29 $(J) Lower bounds based on Equation4.30 and J upper bounds *This bound computed using EiR’Z4 = 1.32625 CHAPTER 5 HARMONIC WAVES IN LAYERED COMPOSITES 5.1 Introduction and Literature Review This chapter discusses the application of the integral equation lower bound technique to the problem of elastic waves propagating normal to the interfaces of layered composites with periodic structure. This problem was proposed by Lee [10] as a test case for variational methods to be applied to composites with periodic structures, and has been the subject of a number of papers [6,10-15]. Composite materials are highly dispersive in their wave propagation behavior, but Floquet theory shows that they admit certain stable wave systems, called Floquet waves, which do not disperse and retain their form relative to the periodic structure of the composite. In two- and three-dimensional lattices, the wave forms may only be determined approximately. However, for the one-dimensional composite formed by parallel plates, the wave forms for waves travelling normal to the interfaces can be determined exactly, and thus provides an objective measure of the accuracy of approximate methods. In this chapter, numerical results are obtained for composites of two materials for a variety of material property combinations and wave numbers. Upper bounds are obtained by application of a Rayleigh-Ritz technique to the Sturm-Liouville eigenvalue problem with discontinuous 55 56 coefficients and quasi-periodic boundary conditions that arises from the application of Floquet theory to the governing partial differential equations. These are the upper bounds as discussed in [6]. Correspond- ing sets of lower bounds are obtained from the L2 norm of the complex— valued kernel of the analogous integral equation. The actual eigen- frequencies were also obtained from the transcendental eigenvalue equation, and errors in the approximations are discussed. 5.2 Formulation of the Eigenvalue Problem The governing partial differential equation for plane waves is the wave equation in each material subregion 2 2 820 x t v U(x,t) = c (X) 47—1 (5.1) at where U may represent displacement, stress or strain and c is the spatially varying wave speed. Floquet waves are represented by a solu- tion of the form T(QTX-wt) U(x,t) = u(X)e (5.2) where u(X) is a periodic function with the period of the composite, Q is a vector wave number giving the direction and wave length of the wave, and w is the frequency. For displacement waves travelling normal to the interfaces in a layered composite, Q and X are scalars, and insertion of (5.2) into (5.1) yields the eigenvalue problem for a cell of one period, the length of which is assumed to be 1, W\ '/////2\\K// ////% fl//////// &\\\ \SS 0 - C u a . 0. w 58 (Hour-HENOu=o (5.3) u(1) = e'°u(0) (nu')(1) = e'°(nu')(o) The unit cell is chosen to start and end at the midpoint of two layers of the same material, and thus the second material layer is centrally located in the unit cell. Further, we require that certain inter- face conditions hold, namely continuity of displacement and stress, u(x ) = u(x+) and n(x-)u'(x—) = n(x+)u'(x+) (5.4) at the location x of any interface. It has been shown that transformation of this problem to Liouville normal form by a substitution for the spatial variable x gives better results in the application of the Rayleigh-Ritz method for upper bounds [6]. Accordingly we set t = I- 70" n"(s)ds, T = 70' n"(s)ds (5.5) f(t) = 126(x(t))6(x(t)) (5.6) V(t) = U(X(t)) (5 7) With these substitutions, the eigenvalue problem becomes V + wzfv = 0 - V(0)eiQ (5.8) (7(1) = (1(0)e"Q < A —I v I v, 9 continuous at interfaces The problem (5.3) may also be cast in terms of the streSSCI = nu', to be called the dual formulation: 1-- 59 (6"6')' + 426"o = 0 6(1) = e'°o(0); (6"6')(1) = (6"6')(o)-e'Q (5.9) o, 0’10' continuous at interfaces The system (5.9) also admits transformation to Liouville normal form by the substitutions s = %-IX o(t)dt; S = 4} p(t)dt (5.10) .n 1 . | 0(s) = s26“'(x(s))p"(x(s)) (5.11) h. w(S) = o(x(s)) (5.12) With these substitutions, the dual problem becomes W + wzgw = O w(l) = w(O)e'Q (5.13) 0(1) W(0)eiQ w, w continuous at interfaces Both the primal formulation (5.8) and the dual formulation (5.13) have the same eigenvalues, but their eigenfunctions differ. It has been shown in [6] that, depending on the relative magnitude of the discontinuities in the material constants,these formulations give differing eigenvalue bounds when a Rayleigh-Ritz procedure is applied to the corresponding variational formulations. This will be discussed further in section 5.5 on numerical results. 6D 5.3 Integral Equation Formulation and Iterated Kernel Trace The complex valued Green's function for either primal or dual formulation is the solution of the boundary value problem for each wave number 0 7‘ 3111 2 3—3 (s,t;Q) = 6(s-t) as G(l,t;Q) = e'QG(0.t;O) (5.14) 3% (l.t;Q) = e'Q 32-(o.t;0) giving the solution 1'0 '(S_-§'% - ——§—1—d—_2— 0 _<_ S i t _<_1 l-e (l-e ) G(s,t;Q) = 1 (5.15) 10 1'0 e (SIB) - 6 i0 2 0 §_t 5_s 5_1 L 1'9 (1'9 ) Note that G is hermitian, that is, G(t,sul) = G(S.t;Q) (5.16) Using this Green's function, one transforms the equations of the form of (5.8) and (5.13), -U = wzhu u(1) = e'Qu(0) 0(1) = e'Qa(0) (5.17) 61 to the integral equation u(s) = 02 70' G(s,t;Q)h(t)u(t)dt. (5.18) The kernel of this equation, Gh, is not Hermitian, but an integral equation with a Hermitian kernel is obtained by the substitutions kQ(s.t) = mama) W (5.19) giving 2(5) =62 70' kQ(s.t)c(t)dt 2 02(02. (5.20) Further, the kernel of the second iterate K6 is (2) - 1 _ 1 kQ (s,r) — A) kQ(s,t)kQ(t,r)dt - {)kQ(s,t)Eb(r,t)dt (5.21) giving tr(KZ) = f1k(2)(s s)ds Q (1 Q ’ = 70' folle(s,t)|2dtds. (5.22) For the problem at hand, we define €10 1 C = ' (l-eiQ)2 = 2(1-cosQ) (5'23) and thus IkQ(s,t)|2 = h(s)h(t)C(C-|s-t| + (s-t)2). (5.24) The value of tr(Kg), which is the same as the L2 norm of the kernel kQ, will be obtained for the sample problem in the next section. of equ 511) f N119 the i0) 1e pr 91 62 5.4 Sample Problem We consider the sample problem of a layered composite consisting of equal amounts of two materials. The unit cell stiffness and den- sity functions, then, are, for 0 §_x §_l 1 1 9 IX - -l > n(x) = ' f 2f (5.25) 1 1 p1 1" ‘ 2' > 4 00‘) = 1 1 (5-26) 02 'X ‘ 2" < 4 We then non-dimensionalize the problem by introducing the parameters _ _ 2 _ 2—-— v - n2/n1 6 - 02/01 v - w o/n (5-27) where E and 5 are the average material properties. We then analyze the problem for the case fi'= EI= l, the frequencies for other cases following from the definition of v. Noting that the Liouville transformation changes the relative lengths of the two material regions, in either case we are led to a problem of the form (5.17), here repeated for convenience, and using . 2 eigenparameter v , U +vzhu = 0 u(1) = e'Qu(0) 0(1) = e'Qa(0) (5.28) ‘where the coefficient function h has the form lt-1§> 1 h2 |t - §1 < For the primal formulation, 1 b = (l + Y)- h - (1 + 3 2 1 T Y) /(4Y (1 + 6)) h2 = yeh]. For the dual formulation, b = 8/(1 + 0) h1 = (l + Y)(1 + 6)/4Y9 h2 = h1/v0. In either case, one obtains the same value for tr(KQZ), namely 2 _ 2 tr(KQ ) — h1 oflra or in terms of the parameters 0, y, 0: tr(KQz) = (111)2[(2+C°$Q)(1+v)2 + ( 1-COSQ Either of these are seen to reduce in the continuous case to “éFJ nfl;d [(6C-l)(l + ( h 2 1 h )2]/l92(l-cosQ) tr(KQZ) = %~(6C-1) = 2+°°SQF 12(l-cosQ)72 (5.29) (5.30) (5.31) h - 1)6)2 + b2(l-b)2(l - E$2] (5.32) (5.33) (5.34) The transcendental equation for the eigenparameter v in terms of 0, y and 0 is 4776 cosQ = (l + 775)2cos€]v - (l - JT§)ZCOSEZv (5.35) 64 l «177' /___:T where E = —-( __X. + .LEL_T (5 36) 'I o 2 1+9 1+6- andg =l(’—'+Y - “HY-1) (5.37) 2 2 1+0 1+6_1 5.5 Numerical Results Numerical results have been obtained for the sample problem on the Cyber 750 computer of the Computer Laboratory at Michigan State University. Upper bounds for the eigenparameter u were obtained by applying a Rayleigh-Ritz procedure to both the primal and dual formula- tions of the problem. For this procedure, test functions of the form u(t) = g C NW“)t = n with M ranging from 0 to 6 provide sets of l to 13 upper bounds from each formulation. Preliminary lower bounds to the least eigenparameter u] were found from (5.33) and the lower bound (2.8). Increasing sets of lower bounds were obtained from the primal and dual upper bounds and the lower bounds (2.9). Additionally, lower bounds based on (2.9) and a merged set of upper bounds chosen from the primal and dual upper bound sets were also computed. The actual eigenvalues were obtained from the eigenvalue equation (5.35) by using the best bounds obtained, searching these intervals for a sign change, and then using the secant method to converge on the roots. These roots were then used as upper bounds to generate the best possible lower bounds obtainable from Equation (2.9). 65 The lower bounds for the least eigenvalue for various combinations of y and e are in Tables 5.1 and 5.2 for Q = l and 3 respectively. These figures contain the actual eigenvalues, the preliminary lower bounds based on (2.8), the better of the primal and dual upper bounds based on test functions with M = 1 (three term approximations), and the lower bound based on the best bound set obtained from the primal and dual upper bounds with M = 1. Lower bounds for the second eigenvalue for the same values of y and 0 as Tables 5.1 and 5.2 are presented in Tables 5.3 and 5.4 for Q = 1 and 3 respectively. These figures Show the effect of using more and better upper bounds in the lower bound (2.9). The figures contain the actual second eigenvalues, lower bounds based on the better of the primal and dual upper bounds for M = l and M = 6, and the best possible lower bounds for the corresponding number of upper bounds. A detailed examination of the upper and lower bounds for y = 10, 0 = 0.01 and Q = l and 3 are presented in Tables 5.5 and 5.6. The tabulated numbers, with percent relative errors, are: the actual eigenvalues; primal and dual upper bound sets; corresponding lower bound sets; lower bound sets based on choosing the best upper bounds from the primal and dual formulations; and the best possible lower bounds for the corresponding number of upper bounds. 5.6 Discussion of Results An examination of Tables 5.1 and 5.2 shows that the preliminary lower bounds are very accurate for Q = l and less so for Q = 3. This is due to the relatively wide separation of the first two eigenvalues 66 in the former case and the lesser separation in the latter. In fact, the continuous cases, where the product ya = 1, produce the least accurate lower bounds. This defect is remedied by use of some upper bound information, and the lower bounds obtained using this informa- tion have less than 0.4% error for all combinations of y and 0 studied. An examination of Tables 5.3 and 5.4 on the second eigenvalue shows the marked increase in the accuracy of these lower bounds as more upper bound information is included. Of particular note are the cases where ya >> 1, where the bounds based on the variational upper bounds are very poor, indicative of the poor performance of the upper bound techniques in these highly discontinuous cases. Clearly, the use of improved test functions with appropriate interface behavior may be worth the additional computational cost in such cases. The use of integral equation variational methods would also be of aid here, as is seen in Chapter 3 and Chapter 6. The best possible lower bounds shown in Tables 5.3 and 5.4 Show how much better the lower bound technique could perform with good upper bounds, and thus provide a measure of how poor the upper bounds are. Tables 5.5 and 5.6 allow comparison of the primal and dual upper bounds for the first seven eigenvalues. In Table 5.5, it is seen that the primal bounds are better than the dual, with the exception of the highest frequency estimated in each line. However in Table 5.6, it is seen that the primal is better for the odd numbered eigenvalues, the dual for the even. Similar analysis was applied to many combina- tions of material properties and wave numbers with the result that neither formulation was ever superior for the bounding of the whole 67 spectrum, further there were cases where the performance of primal and dual formulations were the opposite of that noted in Tables 5.5 and 5.6. Also of note in Tables 5.5 and 5.6 is the performance of the lower bound technique, particularly when the better upper bounds were selected from the primal and dual bounds to obtain the lower bound sets. When these lower bounds are compared in accuracy to the best possible lower bounds, it is seen that, for at least the first half of the number of frequencies estimated, the performance of the lower bound technique is excellent. This reflects the accuracy of the correspOnding upper bounds in this case with a relatively mild dis- continuity, Y0 = 0.1. The results of more highly discontinuous cases which are not here tabulated are quite poor due to the inability of the Rayleigh- Ritz method based on the differential equation form with smooth test functions to cope with the important interface conditions and con- comitant curvature changes in the true eigenfunctions. The mixed variational scheme of Nemat-Nasser [13] deals more successfully with these cases, and upper bounds he obtained, used in conjunction with the lower bound technique developed in Chapter 2, give better results. It is quite possible however that the more standard upper bound methods using either the integral equation formulation or improved test func- tions could perform as well as has been seen in Chapter 3 and Chapter 6. These methods were not tested on the case of harmonic waves. 68 .mgogcm m>¢umpmg «cwucwa mum mmmmgucogma cw mwgamwu cowuwswxogaam snug : co sumac ucaon Luna: 0 Acyp mucaon Long: a co cmmmn canon waop u Acv> mapm>cwmvm —o=uuu u 9 AK-m_-V~mfim_om_. A_oo.-v~o~_m.. Auoo.-vemmmop. Amos.-vmmmnm_. Ammo.-ve_oma_. Anma Am¢._+vo_m~ap. Ao~8.+vemomm_. Repo.+vnmmmm.. Awoo.+yshahmp. Ao.ovo~owmp. Amy». oo— AK-u¢-VmKKm_oa_. Awoo.-vom~,m_. A_mo.-vamwmm_. A588.-vawKop. Aeeo.-vmmmsa_. Roy» mmhm_om_. ¢a~_a_. coama_. _oanm_. omcmap. a ((x_oo.-vuwvmmm. Awoo.-v.mmmmm. Aooo.-v_omoKm. Amos.-ve¢ackm. Amoo.-vpKKcKm. Ammfl Aowm.+VPemoom. AoKK.+VKom~8m. Ao_o.+va8oKm. Ao.ovoemehm. Amoc.+vmmweKm. Amy». op Awoo.-vowemmm. Aeoo.-v~_mmmm. Ammo.-vwm¢oKm. Aces.-vmosesm. Aeeo.-vmmmesm. ona omemmm. Nmmmmm. mwooKm. oomeKm. oaseum. p .Koo.-vm-mmm. Aooo.-v.oewmm. Amoo.-V~Kmamm. Aooo.-v_oe~am. Anoo.-vmmwmwa. an!» Aepo.+vmmemm8. Ao_o.+vemm~ma. Ao.ovoooooc.F Aopo.+vemmwam. Ae_o.+v~m¢awm. Amvm. _ Apmo.-vummwm¢. Ammo.-vmmp~mm. Aces.-vmommmm. Ammo.-vmm.mma. A.mo.-vumommm. Acv> mmmamm. chwmm. oooooo._ Kmemam. oawmmm. > oo. o_ P _. _o. o » _ u o .m:~w>:mmpm ammo; us» ca mccsom —.m m4m=mawm Fungus E» A=V> fl? Amoco.-v_88Kmm. Amos.-vomwmoe. Ac..-vommpme. Amm.-vmmm~mm. Lem.-vmmo~mm. Amy» A8.K+VNKKK~5. A_.m+vm~wm~e. Amm.vommmme. Amo.+vopmemm. Ao.o+vamoemm. Amy». oo_ A_oo.-vmmoKam. Ace.-vmmcwoe. Am.~-vaKoKe. Am.p_-vommopm. AK.NP-vm_om_m. Rev» mooKom. Kommoe. 8NKpme. Kmmemm. amoemm. a ”moo.-vm_fim_._ Ago.-vo_.mm._ Am..-vuomee.n Lem.-vmoapK._ Amm.-vmwomm.. Amy» A..m+v~om¢~._ AK.¢+V~mmK~._ ANN.+VPSOm¢.P Ao.ovmm¢~K.P Amo.+vaKmo._ Amy» Ace.-vmeom_.P A__.-vKNONN._ A..m-vmmmoe.r AK.NP-memom._ Am.__-v_o_Om._ ADV». op 8_Km_._ mmFNN.. 88858., wmch._ _K888._ 5 Ac..-vpmome.~ Amp.-vmmm_m.m Aem.-vmmmmm.~ Amp.-vmmm_m.~ Ac..-v_mom¢.~ Amwm Amm.+v__Fec.N ANN.+vommmm.~ Ao.ovooooo.m ANN.+Vomm~m.~ Amm.+v__Fe¢.~ “Myw Am.~-vmmKKm.N Ap.m-vmawme.m AK.~.-vmoapo.~ A_.m-vmmmm¢.~ Am.~-vm~KKm.N onm. _ ~KNM¢.N moopm.~ ooooo.m m88_m.~ ~KNM¢.N > > cop op P P. .o. a m u o .mzpm>:mmmu ummmA mg» 0» mccaom Logo; N.m m4m:mmpm unoumm _oauum u 9 salmon: :o...§~..~ 29-338; 5:325; 39-58.. 3:» Aoa-vsmmm.o A8.-.namp._ Ase.-chKo._ Am_.-v¢_eo.. Ana.-.amgo._ An_9m Am.m-Vmew.m .cm.-.muop.~ Ao._-vmmKo._ Ao.~-v_-o._ A_.~-voe~o._ any». oc— A.e¢-vo.am.o A.mK-voome.o A.._-V¢mmo.o Aa.o-.~_um.c A..~-.m¢~o._ Amy» ammo.8 cKm~9.~ «Kemo._ momeo._ m_oeo._ > Ape.-vm_Kp.e Rec.-vm-K.m AND.-VK58_.m Ame.-vuemo.m Ame.-.KK~o.m Ampmw A.me-vomoM.m A._¢-Vweov.m Amm.-vooe_.m Amc.-vsmmo.m Ao._-vmm~o.m Amp.» Aew.-vmo~_.o AK.~-VmcKm.m Am._-vms__.m A..N-v_msa.~ Ac.~-vaem.~ Amvm. o. A.¢K-CKQK~.P A.Kh-vemmm.9 Am.m-vmmmm.~ AP.~-VFme.~ A¢.m-vmo_m.m Amy» KGNLP.8 moeNK.m ~8m8_.n NcKmo.m mmmmo.n > “Po.-vmshe.m Ame.-VNQOm.m Ame.-vo_m~.m Awe.-vmeom.m Ape.-vmuse.m Anpvn. AK.o-Vooee.m Amm.-voKKe.m Awe.-.o_m~.m Amm.-voKKe.m Ak.o-vaooe.m Am_9m Ao.p-ve_~e.m Am._-vomme.m A_.~-veek_.m Am.9-vo-e.m Ao._-ve_~e.m Amy». _ “...-veoam.¢ Aw.w-v_o~o.m A_.~-veeK_.m Am.m-v_o~c.m A__-V¢oom.e Ammfl mmKKe.m oemom.m a_mm~.m oemgm.m mmKK¢.m 9 oo_ op _ p. —o. > F u a .ma—e>cma.m ucouom may cu mvcaom szoA m.m u4m

:mm.m ucoumm .mauuo « 9 35:88.... 25:55.85 855.-....55 355.8285 .85.-..585 .55» 35-5.5555 .8522... 85:58.55 38:58.85 .85.-..585 .2.» 8.32:5. .5..-:85.N 85.55255 25:58.85 35:55.85 .2» 55. 85.8585 8.5285 3.5-5885 25:55.85 85:55.85 .2» N885 .285 5.255 52555 5.5.85 .9 25:58.85 35:588.... .555.-.5.N5.N 38:885.. 38:558.. 5% $252.: £1.58... .55.;2N5N 355.588.. 355.588.. .2.» 22.55255 8.7.55.5... 35:58.55 85:588.. ..5.-55..5.. EN 5. 838.5. 85.5355. 3.3885 85:588.. £2588. .2» 5255.5 55 .85 N .5N5.N 528.. 852... .9 38:82.5 855.5825 355.5825 88:522.. 355.588... .5...» $5.558: 35:55.25 $8.558: 35:55.85 35.555255 .2.» .8585: 35:585.... 85:55.85 35:558.... .8585: EN . .55-:585 8.2.5.5.... 85:55.85 3.2.5.5.... 9.5-5.5525 EN :58... :83 5.825 ..NNN... :58... 9 8. 5. . .. .5. 5 9 e.m m4m:mm.m vacuum mg» o» mucaom 85:54 '72 mas—59555.5 paauua u .9 m5czoa .555: .555 :5 55.55 m5caon .535. u a. mu:.5>cmm.o .mauuu :5 55555 «55:55 .535. 5.5.mmon ammo «45m 555555 .555: .55..5 :0 55555 m5caoa .535. u 5. mcopuapgom cot—23:50» :25 .5.—5:33 Luna: n a: .5=5\.5s.55 so.» .55505 .555: umaa no 55555 «5:555 .535. u a. 55.95.:Ego. .55—La .555555 .555: u 5: .mN-V...5. M.N-Wm...5 M.- 555.5 .5.5- m.... ...n- 55... 55.- .5.... M5 5N.N5.. N555.- m55...5.5 ..5 .5.-5.5... 5.- .55.. N.- 555. 5 ..N- ....5 .5.- 5.... 5..- .5.5.. .. .- .m.55.N .555.- ..5...m.5 .. ...-5555.. .5.-55.5.. .5- .555. . 5.- 55... ...-5N.5.. ...5.-.55..n Mn. 5- ..5.N5.N .m555.-5555...m.5 5. . .5.-.55... .5.-5N55.5 N.-..55.5 5N- 555.. .5.-5.5... ....-..555.. .. .- .55.55.N ..555.-.5.5...m.5 .5 .... .5... N... 58.5. ..... ..5.. 5... ..5.. .... N..N.. Mm. . 5mNm5. m .55.. 5..5..m.5 55 ..5. 5..N. .... ...5. .5... .N5.. mN.. ..5.. .... N..N.. .5. . ..5N5. 5 N55.. 555.....5 .5 .N- 5... 5 m.-.NN..5 .5.5-55.5.5 ...- N..... .5. 55..5.m .55.- N55...5.5 ... .m- ..5. m 5.-...5.. .5N-55.... ..5- 5.N5.n .5. N- 5...5. N ..5.... 5 5. ..5- 555.. 5.-.555.. M5m-w..N.. ...- .5.5.m .5. .- 5Nm55. N .55: 555...m. 5 5. 5 ..m- ..5.. 5m-...5.. 5N 55... .N..- ..N5.m .5. N- .55.9 N .55.- m.5...m. 5 .. .N....5.5.. 55.. .55.. 5N.. ...N.. 5. . .5.55.5 559 . 5mN5..5.5 55 .5.N..5.5.5 N... 5.5.5 5N.. N..N.. 55.. m..m5.n .55.. 5555....5 .5 .5.-...... .5.5-.N5N5.m .5 N-.N5.55.N ..55.- 5.....m.5 ... .5.-5555.. .5.-5N5...m 5.5-55.5.5.N .55.- ......5.5 5. ...-...... ..N-vm..m.m ..5-.N.m...N .55.- 55......5 5. m .5.-...... .5.-.N5.... .5 . 55..5.N m55.-.55....m.5 .5 .5....5.NN.. ..m.. ..5.5.. ..5.. 5.85....5 55 ..5...NN.N.. .mN....mmm5.5 555.. ..55..m.5 .5 .5..-.m5.55.N ...5.-...mm..5.5 ... .5.-...55..N ...5.-...mm..m.5 5. .5.-5......N ...5.-...mm..m.5 5. . .5.-...55..N ..5.-.......m.5 .5 5m5...5555..5.5 55 .5.5...5555..5.5 .5 .5.5. 5.... 5N.5. 5.5.. NM... 555N.. 5m5N5.m 555.....5 .9 5 . 5 m . m N . . . .5.5 u 5 .5. - 5 ...: 55.....5 5;. 5. .5555. m.m m4m 0 m(S-ITI The L2 norms of the corresponding kernels are given by ||km||2 = [0‘ 4} rd(r)Gm(r,s)sd(s)dsdr (6.19) (6.20) (6.21) 80 For the density d(r) given above, these norms are 2 d ||k0||2 = 5%?-{1 +-©2 - 1)[(1 - a4) - 8a4(lna)2 + 464ina] + 86(1 - e)[a2(1 - a2) - 2a4(ina)2 + 2a4inaj} (6.22) 2 d 2 ||k1l| = T%—-{1 + (62 — 1)[6(1 - a4) - 8(1 - a6) + 3(1 - a8)] 4 2 4 + 66(1 - e)a [24lna - 24(1 - a ) + 6(1 - a )1} (6.23) and for m > 1 d 2 2 4 2m+4 llk I'2 = I { m + (62 _ ])[I-a - 1' 'a m 4m2(m+l) 4(m+1)(m+2) 4 m+2 4m+4 2- 2m + '2Ifi%ri) 3 T 9(1'9)a 2m 2[] 3-2m ‘ (1‘32) 2m+2 l- -a 2m,2 1} (6.24) Note that in these formulae as a approaches 0 or 1 or as e approaches 1 these formulae reduce to those of the continuous case, namely ||km ((2 /(16 (m+1)2(m+2)); m 0, l, 2, ... (6.25) vzrdu have as their solutions The differential equations Lm[u] linear combinations of Bessel functions of order m in each material subinterval, namely 81 C J (or) + C Y (or) 0 < r < a u (r) = 11 m 12 m (6.26) C21Jm(8r) + C22Ym(6r) a < r < l where a = J3; °v and B = J3; ~v. Using the four conditions, two boundary and two interface lim rum(r) = O, um(1) = O r+0 . (6.27) um, rum continuous at r = a a nontrivial solution is obtained when u satisfies the equation ( ) J$(aav) Jm(a/§av)Ym(a/6p) - Ym(a/5av)dm(a/§u) J aav = - - , . , m ,6- Jm(a/§au)Ym(a/e—v) — vm(a/6a6)am(a/66) (5 28) In the continuous case, these reduce to the familiar equations J (u) = O, m = 0, 1, 2, ... (6.29) Upper bounds to the eigenvalues may be obtained by applying Galerkin's method to either the integral or differential form of the problem. For the problem at hand, both methods were applied, using approximate eigenfunctions of a form appropriate for the axisymmetric modes: n u = Z aisininr/inr in the differential form (6.30) l or w = er(r) u in the integral form. (6.31) 82 The matrices for the method, in terms of the usual L2 inner product are (H).. (M).. (G) II A 2 J- C V 13 1 J H A ‘6 an... €— V II < rd(r)ui, u. > (6.32) U J 'ij <¢1,K¢j>, The approximating eigenvalue problems then are: Ha and Ma 32Ma in differential form (6.33) 326a in integral form. 6.5 Numerical Results The numerical results in this section include the actual eigenvalues and both upper and lower bounds obtained by implementing the following procedure on the Prime 750 computer of the Case Center for Computer Aided Design at Michigan State University: 1. Generate increasing sets of upper bounds and corresponding lower bound sets from the Galerkin methods and equation (2.9). Use the bounds generated to solve the actual eigenvalue equation by a. finding an interval where the equation has a sign change b. applying the secant method until the root is found Use the actual eigenvalues to generate the best possible lower bounds obtainable from equation (2.9) and to compute relative errors. 83 Table 6.l contains lower bounds to the least frequency for a number of internal radii and density ratios. These bounds were obtained from the kernel norm alone, using none of the upper bound information. Tables 6.2 through 6.4 contain bounds to the spectrum of 3 cases of internal radius and density ratios, as well as the actual eigenvalues and percent relative errors of the bounds. The numbers tabulated are the upper bounds from both the integral and differential equation based Galerkin method, the lower bounds based on the integral upper bounds. and the best possible lower bounds, based on using the actual eigenvalues as bounds. 6.6 Discussion of Results Examination of Table 6.2 shows that the upper bounds from the integral equation formulation are clearly superior to those from the differential equation formulation. However, both methods have dif- ficulty coping with highly discontinuous cases. The superior per- formance of the integral formulation stems from the iteration of the test functions with the kernel, yielding essentially improved test functions having more of the properties of the actual eigenfunctions at the interface. Examination of Table 6.l shows that the preliminary lower bound based on knowledge only of the kernel norm are quite accurate, with relative errors ranging from 0.4 percent to 2 percent in the cases studied. The accuracy of this bound depends on the spectral structure, and for the cases studied the least frequency is well separated from the higher frequencies, yielding good lower bounds. 84 Examination of Tables 6.2, 6.3 and 6.4 shows that accurate lower bounds can be obtained provided that good upper bounds are used. The accuracy of the lower bounds for higher frequencies depends primarily on the accuracy of the upper bounds to the lesser frequencies. Given a good upper bound to the lower frequencies, the effects of eigenvalue spacing can be seen, especially in Table 6.2, where the actual fre- quencies come in closely spaced pairs after the isolated initial frequency. The lower bounds to each of the close pairs are of com- parable accuracy,though the first lower bound obtained for the higher of the frequencies is better than the comparable bound for the lower of a pair; this is the effect of clustering on lower bound accuracy. 85 F ‘U (P F P2\ § pl ,. i W“! A ,4 R' l J QJJJJQJJ JJJJJJJ JJ J JJ \\\\\\\\“\\-. (ru')' - QE-u + w2 gF-u = 0 u, ru' continuous u(O) finite, u(R) = 0 m = # of nodal diameters FIGURE 6.l Circular membrane with radially stepped density 86 TABLE 6.1 Bounds to the Least Frequency of the Two-Piece Circular Membrane 10 3 1/3 1/10 v 2.43282 2 42571 2 33789 2 09196 UD 2.44599(+.54) 2.43715(+.47) 2.34503(+.30) 2.16144(+3.3) '1 UI 2.43496(+.09) 2.42765(+.08) 2.33848(+.03) 2.10638(+.69) L0 2.40815(-1.0) 2.40067(-1.0) 2.30659(-1.3) 2.05683(-1.7) v 3.02808 2.83469 1.91760 1.61114 U0 3.39475(+12) 2.97186(+4.8) 1.92658(+.47) 1.63398(+1.4) .5 UI 3.06703(+1.3) 2.85978(+.89) 1.91852(+.05) 1.61259(+.09) L0 3.00730(-.7) 2.804397(-1.1) 1.90554(-.63) 1.60399(-.44) v 3.72529 3.09024 2.01280 1.84025 U0 4.36498(+17) 3.18748(+31.) 2.01352(+.04) 1.84090(+.04) .7 UI 3.88142(+4.2) 3.11528(+.81) 2.01288(+.004) 1.84027(+.001) L0 3.67186(-1.4) 3.03337(-1.8) 1.99646(-.81) 1.82659(-.74) v 3.87949 2.81390 2.25019 2.19328 UD 3.92151(+1.1) 2.82612(+.43) 2.25681(+.29) 2.19952(+.28) .9 UI 3.89011(+.27) 2.81613(+.08) 2.25132(+.05) 2.19434(+.05) L0 3.80068(-2.0) 2.77939(-1.2) 2.22614(-1.1) 2.17003(-1.1) v 3.28768 2.62759 2.32573 2.29739 UD 3.30066(+.39) 2.63647(+.34) 2.33318(+.32) 2.30471(+.32) .95 UI 3.29001(+.07) 2.62913(+.06) 2.32701(+.06) 2.29865(+.05) L0 3.24764(-1.2) 2.59826(-1.1) 2.30030(-1.1) 2.27230(-1.1) = actual least frequency V U0 = upper bound, differential formulation 01 = upper bound, integral formulation L0 = truncation lower bound, Equation 2.8 £37 E.Ho—.= + u E + u: + usage» povucocoaxw c. ago —oo. v mgoggu .mgoggu u>.uupus vacuum: was mmmogucogua c. mwgauvu a. N co.uo:om new > sag» mezzo; gmzop u a; o. N ccvucaam can _3 son» avenue 5936. u .4 co.uo_:ELcw pusamucv .mvcaon Long: a _a =o.uap=sgou —o_a:oguu~—e .mccaoa Luna: u a: .o_-vm.mp «.m-vm.~. Am.~-.c.o. Aa._- m_.m AGN.-.¢~.m A... -vce. c-mo-.moo~._ as Ame-vc~.~ me-vso.~ A--vwm.~ Aeu- 9..“ M_.enw~m.m m. N- on. mp°.- ooo~._ _s op A_m+v_.¢~ .am+ve.o_ A_.+Vm.._ A...+ _~.a e.m+ ma.m ...+ 5.. e-m¢+ ooo~._ ~= A~m+vp.m~ Ane+vp.o~ A~¢+v~.m_ A..+V~.o_ .o.+v.o.o .m._+..m. ....+vomc~.. a: Am~-Vm.c_ ..- mm.» ..- em.m Aa._-wmo.m An~.- «4. ..cc.- aco~.p as Aeo-vmc.m _m- ~o.m me. am.. Am_- no.4 Ao_- co. Aoo.- amo~._ .4 m Po+ m.o_ m~+ N... oe+ N..m a..+wm . No.+ ooa~._ _= ca+ o.m_ _m+ ..v. vo+ v..o ..m+ . . ~.~+ “an~._ a: e e e e e c e . .._-.mm.m m.. "m." Ac..- mm.m .m.—-vo e _o.-vmoo~._ as Aom-vem.. NN- a... A¢_-w. m Kc.- mmo~.. _4 e M? $_. n. cm+ mo.m m.~+ mm.. co.+ m~o~.. _= mm+ m. «p mc+ o~.a n.~+. .. o.m+ mmo~.. a: m_- eo.~ .m.o- Kn.m .~.~- e um- oo.e Apm- “a.” .a_- n e e o m c m o N No.1 vmo~.p ma No.1 umc~.— ~4 wmcom O m n m m «8+ we. a A~.n+w o. m—.+ ~ao~._ _= m~+ aa. a A__+ a. m.m+ _Ncn.. a: Am.o- kn.m Am.a- me. we.. mmo~.p as Ame. o_.n AN”- mo. mo.- mmo~.. _s N Am. 8+ 05. ~¢.+ o—_~._ _a Akp+ mm. .m_+ m~am.. a: Am.a-vmc. Ao~.-V~eo~._ as Apm-v_~. o~.- ~e¢~._ .3 ..~+ m~m~._ _a P .-+vmonm.. a: on: ...: 3: 2.... o; 3... 38.. r, 223. po. u a .p. n a guvz Ezguumam may go» mvcaom ~.c m4m cub. : + e + u + a u — .u . poo m.~ copuuaau vcm .us—mg acmugmn ago mummgucmgaa c. mmgzm.u > so.» mucaoa Lose. a a. a.~ =o.um=au use .3 so.» mucacn Lyso— u .4 :o.uo.:Egcm .ogmmu:. .mvcaoa swan: u .3 co.uu.=sgom .m.u=ogmum.u .mccaco Lona: n a: ...-.N.N. ....-.m.mm .N.m-.m.Nm .N..-.o.mN .mn.-.me.¢. .eo.-...m... ..-NN-..N.Q.¢ N. .N.-.¢.NN ...-.e.mN ..- e..N .....-w..MN .N..- am... .Nm.-.mN... ..o.- ON.o.e .. o. ..om+..oe .mm¢+.ecN mm+ a... ....+ ...N .ooo.+ m..m. ..oo.+.o.m... .m-NN+ «N.o.e .= .mm..+.oo. .mm..+..om .Nmm+.a..N .NN.+.N.Nm .m.N+.on.a. .en.+.¢mn... .mc.+.mm.o.¢ a: ON- ..mn ..-...NN .o.m-...¢N «.N-.on.m. .cm.- NN... .moo.-.NN.o.e N. mm- ..m. N.-.e.a. ..N- .... m..- ea... o.N- No... .mo.- oo.o.e .. m .omm+.emm ..mm+ o..m. .m.+ N...N No.+ NN... ..-..+ .N.o.. .= ..oa.+.o.e ..eo.+...amN .Nom+.N.eN. ..N+..o.e. ....+.No.c.e a: ..- ..NN ..- N.NN ...e- Na... m..- «N... .c.- ON.o.. N. .m- N... 8.. ..c. ..- mo.o. ..m- cm.o. No.- .aoo.e .. e .Noa+ N.eoN N¢+Vnm.mN ....+ me... .....+ .N.o.e .a .mmm.+ ..NN. ..ao.+.m.eNN ...+ N”... .om.+ NNN°.. a: ...-.N.NN .o.m-.wo... .m..-.m.... .Nc.- ...e.e N. .me-.om.e. .mN-.ae.¢. ....-.No.o. .m..- occc.e .. m N.e+.eo.oa a.e+ .m... ..oo.+ .N.o.e .= mee.+.o.mmN .cN+ mc.¢m .N.N+ aNN... a: .o.a-.co... ....-.mm.o. .Nc.-.m.oo.e a. .4.-.mm.N. .8..- .m.c. .N.- .moo.. .. N .m.+ «N.N. mc.+ «..9.4 .= .eme+.me.mm .m.+..m¢m.. a: c.4-.m¢.o. o..- Nam... N. 8N-.amm.m a..- .mmm.m .. . ..N.+ ow.o.e .= .mm+ .cc... a: o.me N... .... N.NN m..m. ...... ..N.c.. .9 ”sea. . N c m e m N . . ooo.o. u m ... u a cup: Eaguumam on» so; mccsom m.m m4m.uu—mc acmucoa m N Numwgucugaa :. Nugaopu > ..5.... 3.53 .28. o.~ :o.uo:au was .3 so.» chacn gaze. a. ~. co.uopasgo» .mcmwuz. .Nucaoa song: a .2 co.ua.=sgom .N.u=mguw».u .chaoa song: a a: ..- N.N. .N..- N.NN .N.N-...NN .N..- N..N NN.- N..NN MNN.- NN..N. .N.- NNNN..N N. NN- N.NN .NN- «.NN .NN-...NN .NN- N..N N.N- NN.N. N..- «.N... N..- NNNN..N .. N. NN+ N... NN+ N.NN NN+ N.NN .N..+ N.NN ..+ N..NN NN.+ ....N. .N.+ NNNN..N .: NN+ N.N. NN+ ..NN NN+ N.NN .N.+ ..NN NN+ NN.NN NN.+ N.N N. N..+ NNNNN.N N: ..N-.N.NN ...- N.NN .N.N- N.NN .N.N- N..NN ...- NNN.N. NN.- N..N..N N. .NN-.N.N. NN- N.N. NN- N.N. .NN- NN.N. «.N- NNN... NN.- NNN...N .. N NN+ ..NN NN+ N... .NN+ NN.NN .N..+ NNN.N. N¢.+ NNN...N .: ...+.N.NN .NN+.N.NN .NN+.N..NN ...N+.N.N.N. .e.+.NNNNN.. N: ...-.N.NN ..- N.NN .N.N- NN.N. N..- .NN.N. .N.- NNNN..N N. .NN-..... N.. ..N. NN- NN.N. «.N- N.N.N. N..- NNN...N .. . NN+ ..Ne .N+ .N..N .N..+ NNN.N. .NN.+ N..NN.N .: NN+ N.N. Nn+ NN.NN .N.+ N.N.N. .N.+ NNNNN.N N: ..- N.NN ...- NN.N. ....- N.N... ....- .NNN..N N. .N- N... ...- N.... .N.- NNN.N. .NN.- «N.N..N .. N .N+ .N.NN .N.N+ NNN.N. ...N+ .NNNN.N .: NN+ N..NN .NN+ N...N. .NN+ N.NNN.N N: ...-.NN.N. ....-.N.N... .NN.-.N.N...N N. .NN-.NN.N .N- N.N.N N..- ..NNN.N .. N .... NNN.N. ...N+ 8N... .: .NN+.N.N... .NN+.NNNNN.N N: ....-.N.N... N.N- NNNNN.N N. .NN-.NNN.. N.N- NNNNN.N .. . .NN+ NNNNN.N .: ..e+ N.N.¢.N N: ..N. ..N. N.NN N..N NN.NN N...N. N¢NN..N .9 Nata. . N N N N N N . . oo. u N .m. u a cup: Escuuoam mg» no» Nccsom e.o m4m